Welcome all!
Hope you are staying safe and following the COVID-19 safety guidelines.
So recently, like everyone else, I had spent a lot of time indoors, and with my background in GIS and Remote Sensing, I was looking for reading materials on Harmful Algal Blooms all over the internet, and found a handful of ‘Open Access’ scientific papers on this matter.
I believe in Open Science and Open Source codes and knowledge sharing, which is why I decided to convert all that I have read over the past few days, into an executable code using Google Earth Engine platform. I hope this is useful.
Before going into the science and the coding blocks I would like to direct GEE first-timers/new-comers to the following links -
- https://developers.google.com/earth-engine/guides — For ‘How-to’ guides on GEE Python installation, and also lots of useful tips on using GEE’s web browser based ‘Code Editor’
- https://github.com/giswqs — For guides and tutorials on Google Earth Engine Python API usage and Qiusheng Wu sir’s own ‘geemap’ python library, which I have used extensively. Thanks for that Qiusheng sir!
THE SCIENCE -
To understand how to map/detect Harmful Algal Blooms, we have to first read what are the physical and biological properties of these blooms, and then proceed to deriving them from the satellite images.
With the risk of turning this medium article into a scientific paper, here is a direct quote from “Fahad Alawadi’s ‘Detection of surface algal blooms using the newly developed algorithm surface algal bloom index (SABI)’ 2010 paper. Linked here as (1)
Algae is a non-planktonic species that performs photosynthesis. It is commonly classified into two groups: macro-algae like seaweed and micro-algae like filamentous cyanobacteria. They are also divided into three groups based on their photosynthetic pigments: Phaeophyta (the browns), Rhodophyta (the reds) and Chlorophyta (the greens), but can also appear in various other colour combinations like yellow or silvery grey. They are capable of regulating their buoyancy to modulate their photosynthetic activity to the changing light levels. This mechanism of buoyancy is either to protect themselves from excessive radiation or to access more light from the photic region located in the top layer of the water surface. Surface emerging algae can experience, similar to land-based vegetation, the “red-edge” effect, where the greatest reflectance slope occurs between the maximum absorption in the red and the maximum reflection in the infra-red (0.7 µm-1.0 µm). The IR reflectance is caused by the combination of decreasing absorption by the chlorophyll pigments in the spongy mesophyll cells and the increasing of absorption by water. Texturally, surface blooms are often characterized by their thick extensive mats or scum that looks like oil slicks or slightly foamy pollution. (1)
This paper by “Singh, Rakesh & Shanmugam, Palanisamy. An optical system for detecting and describing major algal blooms in coastal and oceanic waters around India” linked as (2), is probably one of the most detailed studies conducted on a wide geographical scale. Dr.Shanmugam has some of the highest cited articles and papers on Algal Blooms.
All the graphs and explanation in this paper (2) boils down to this logic -
As the Algal Bloom gets more dense and begins to form a mat-like feature on the ocean surface, the reflectance of the blooms increases in the 700+nm Infrared range.
THE CODE -
With the theory and science established, let’s have a look at how MODIS Aqua Imagery can help in converting these theory into maps.
Although Aqua satellite has been malfunctioning intermittently since August first week 2020 (3), the amazing people at NASA’s Aqua Flight Operations team are working round the clock to keep the 18 year old satellite up and running. This 18 year pedigree is why I chose MODIS Aqua with its extremely large collection of daily imagery at 250m resolution, now at our fingertips via Earth Engine.
To view the MODIS Imagery in False color -
In the Web Based Code Editor — first define your Area of Interest (AOI), this is not a compulsory step since GEE can operate on global basis, but let’s be nice on Google and reduce our computation requirements.
You can simply define your AOI by using one of the shapes from the panel seen on the browser (attached screenshot below) , or import your own Shapefile and use that. My AOI is off the coast of Mumbai City, India, and for the date — 06th August 2020.
Once you have your AOI defined we can use it to clip all MODIS imagery to that area.
Below I have defined Aqua Imagery collection as the variable ‘AquaBloom’, you can give it any name you desire.
// On Web based JS editor -
var AquaBloom = ee.ImageCollection('MODIS/006/MYDOCGA').filter(ee.Filter.date('2020-04-06','2020-04-07')).first().clip(AOI).aside(print);
// 'aside(print)' prints out the details of the imagery
//// including time of acquisition, bands and more.// '.first()' takes the first image of the Image collection for your //// given dates, in my case the image of the 6th of August 2020.
var falseColorVis = {
min: -50,
max: 100,
bands:['sur_refl_b15','sur_refl_b14','sur_refl_b13']
};
// 'falseColorVis' is defined as the visualization parameters
//// so bands 15,14,13 as R G B respectively.Map.addLayer(AquaBloom, falseColorVis, 'FasleColor4Blooms');// a layer named 'FasleColor4Blooms' will show up on the map.
The reason for choosing bands 15,14,13 as RGB combination is due to the apparent ‘red-edge’ of Harmful Blooms in the Infra-red region of the EM spectrum(1,2). This should give us a nice reddish color for the Blooms, but sediment rich areas, will also shows up in red. This will be removed later.
Bands 15,14,13 of MODIS Aqua are 10nm windows of 743–753nm,673–683nm and 662–672nm respectively, hence the pixels in each band holds mixed values of surface reflectance across these nm windows. A point to note is that, on Google Earth Engine the pixel values are scaled with a multiplication by 10000.
So lets look at the first visualization.
Yep! That’s a massive Harmful Algal Bloom right in the middle of the image. Notice the scale on the bottom-right, ‘50km’ for ‘that much’ length, so this formation, visually, could be around 150kmx250km.
There are other features here too, the light blue/cyan whisks (G&B of the RGB ie., Band 14 and 13 of MODIS) in the water are possibly ‘sediments’ when they are right next to the coast, and may even by other ‘algae/non-algae organic matter’ which are reflecting in shorter wavelengths ie., 650–680nm range, compared to the NIR band, Band 15 of MODIS which is 740–750nm (nanometer).
The whiteout on the left side of the image is sunglint! Surely we have to avoid that too.
Anyway, for some more intuitive understanding of the visualization let’s check the actual band values of the bloom.
By using the “Inspector Tool” and clicking on the Algal Bloom affected pixel, we can see that the band values will form a two spike line if we were to graph them. Bands 10, 15, 16 are having higher values compared to the other bands. Taking advantage of this relationship between the bands we can perform a simple math calculation to bring out Algal Blooms from the image.
// We define the variables 'a' and 'b' as bands 15 and 13, which as we can see from the values, have almost an inverse relationship.var a = AquaBloom.select('sur_refl_b15')var b = AquaBloom.select('sur_refl_b13')var c = (b).subtract(a);var Blooms = ((a.subtract(c)));
By declaring variable 'Blooms' as the above equation, we can see that when the difference between Band 15 and Band 13 reflectance values is large, then the Bloom value also increases, and when its low, Bloom value decreases. The emphasis of this relationship is more on Band 15, since Band 15 is in the Infra-red region, which is indicative of our initial logic derived from reading the papers, that higher the Infra-red reflectance, more thick and harmful the bloom.
When we view the ‘Blooms’ layer with the following code, we get -
var viz = {
min: 0,
max: 30,
palette: ['cyan','green','yellow','orange','red']
};Map.addLayer(Blooms, viz, 'Initial Bloom Layer');
And after some careful ‘inspecting’ with the Inspector tool, we can add ‘mask’ onto this layer, so that we ignore all values above 30 and below 8, since it seems to me, that the ‘Blooms’ variable equation produces this particular range of useful pixels.
The following code results in a masked layer -
var maskedBloom = Blooms.updateMask(Blooms.gt(8));var maskedBloom = maskedBloom.updateMask(maskedBloom.lt(30));var viz2 = {
min: 8,
max: 30,
palette: ['cyan','green','yellow','orange','red']
};Map.addLayer(maskedBloom, viz2, 'Harmful Blooms');
Now that we have our required layer, can we remove some of the noise? The noise is due to the edges of clouds and other atmospheric artifacts.
Thankfully there is a Quality Control QC band, as part of MODIS products, in which all ‘Bits layers’ which are numbered as ‘0’ represents the best quality pixels, so we now mask our ‘maskedBloom’ layer with all the QC layers other than 0, and simply overwrite it, so that we only view the best pixels.
The code for that is the following, and we get -
var qc = AquaBloom.select('QC_b8_15_1km');var maskedBloom = maskedBloom.updateMask(qc.neq(0));Map.addLayer(maskedBloom, viz2, 'QC Harmful Blooms');
We that, we now have a working code that derives Harmful Algal Blooms from MODIS Aqua Level-1 imagery.
This final layer can be attributed as a simple ‘Harmful Algal Bloom Intensity’ layer, with value 8 as ‘High’ and 30 as ‘Extreme’, and the focus is on surface floating Algal Blooms which are already in a pretty mature state of growth.
The ‘logic’ works at any location at any time, but do keep in mind that this is quite an abstracted version of the algorithms implemented by scientists.
So, there is a lot more work that can be done, to remove noise from clouds, sunglints, and more.
I have also tried this ‘logical’ approach on Sentinel-3 Level 2 imagery, and it seems to work well there too.
