Earth Engine
for Dummies
Part III
Michael A Menarguez
mamenarguez@ou.edu
Contents
- Paddy rice mapping algorithm solution (1st part)
- Filtering bad observations (where + mask + bit operations)
- Powerful ee.Image.iterate for algorithms!
- Guided exercise: LST > T for k days
-
Download Results
- Image
- ImageCollection
- FeatureCollection
Paddy rice mapping (I)
1 - Obtain MODIS collection
var year = 2014 //Easy way to change year in the script
var modisCollection = ee.ImageCollection('MODIS/MOD09A1') //Load MOD09A1 Collecion
.filterDate(year+'-01-01',(year+1)+'-01-01') //Filter by year
print('Modis collection');
print(modisCollection);
Paddy rice mapping (II)
2 - Calculate vegetation indices for every image
//This function receives a MOD09A1 image and returns ndvi, evi, and lswi
function getVegetationIndices(img){
var ndvi = img.normalizedDifference(['sur_refl_b02','sur_refl_b01'])
var lswi = img.normalizedDifference(['sur_refl_b02','sur_refl_b06'])
var evi = img.expression(
'2.5 * (nir - red) / (nir + 6 * red - 7.5 * blue + 1000)',
{
red: img.select('sur_refl_b01'), // 620-670nm, RED
nir: img.select('sur_refl_b02'), // 841-876nm, NIR
blue: img.select('sur_refl_b03') // 459-479nm, BLUE
});
return img.select([]) // Implementation improvement. Need to return same image
.addBands([ndvi,evi,lswi]) // Add bands to image
.select([0,1,2],['ndvi','evi','lswi']) //Rename bands
.copyProperties(img,img.propertyNames())
}
// Previous CODE from slide 1 ...
var modisViCollection = modisCollection.map(getVegetationIndices)
print('Modis VI collection');print(modisViCollection);
Paddy rice mapping (III)
3 - Obtain flooding condition for every image
function getFlooding(img){
// Extract VI bands to variables for legibility
var ndvi = img.select('ndvi'), evi = img.select('evi'), lswi =img.select('lswi');
// Calculate flooding mask [(lswi>ndvi or lswi>evi)] = [(lswi>min(ndvi,evi))]
var floodingMaskNdvi = img.select([0]).multiply(0).where(lswi.gt(ndvi),1);
var floodingMaskEvi = img.select([0]).multiply(0).where(lswi.gt(evi),1);
var floodingMask = img.select([0]).multiply(0).where(floodingMaskNdvi.or(floodingMaskEvi),1);
return img.select([])
.addBands([floodingMask,floodingMaskNdvi,floodingMaskEvi])
.select([0,1,2],['floodingMask','floodingMaskNdvi','floodingMaskEvi'])
.copyProperties(img,img.propertyNames())
}
// Previous code here ...
var floodingCollection = modisViCollection.map(getFlooding)
print('Flooding collection');
print(floodingCollection);
Paddy rice mapping (IV)
4 - Fast visualization of results:
// Visual representation using median values for each vegetation index
var palette = ['FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718',
'74A901', '66A000', '529400', '3E8601', '207401', '056201',
'004C00', '023B01', '012E01', '011D01', '011301'];
var studyRegion = /* color: 98ff00 */ee.Geometry.Polygon(
[[[117.333984375, 54.97761367069625],
[109.3359375, 47.398349200359256],
[96.240234375, 31.353636941500987],
[101.337890625, 18.729501999072138],
[119.00390625, 17.811456088564483],
[144.931640625, 32.32427558887655],
[153.017578125, 52.3755991766591],
[139.482421875, 54.059387886623576]]]);
Map.addLayer(modisViCollection.select('ndvi').median().clip(studyRegion),{palette:palette,min:-0.1,max:1},'ndvi')
Map.addLayer(modisViCollection.select('evi').median().clip(studyRegion),{palette:palette,min:-0.1,max:1},'evi')
Map.addLayer(modisViCollection.select('lswi').median().clip(studyRegion),{min:-0.1,max:1},'lswi')
Map.addLayer(floodingCollection.select(['floodingMaskNdvi']).sum().clip(studyRegion),{max:46,min:0},'Flooding ndvi')
Map.addLayer(floodingCollection.select(['floodingMaskEvi']).sum().clip(studyRegion),{max:46,min:0},'Flooding evi')
Map.addLayer(floodingCollection.select(['floodingMask']).sum().clip(studyRegion),{max:46,min:0},'Flooding signal')
Paddy rice mapping (IV)
5 - All together:
**Note that we did not use day of the year filtering to identify the transplanting season yet,we will use LST for this soon... BE PATIENT!
Filtering bad observations (I)
Use MOD09A1 StateQA band and remove cloud, cloud shadow, and snow pixels
- Test algorithm on a single image
- Test algorithm on Image collection
Filtering bad observations (II)
Create a function to extract bits from a position i to j
//Extracted from GEE examples
var getQABits = function(image, start, end, newName) {
// Compute the bits we need to extract.
var pattern = 0;
for (var i = start; i <= end; i++) {
pattern += Math.pow(2, i);
}
// Return a single band image of the extracted QA bits, giving the band
// a new name.
return image.select([0], [newName])
.bitwiseAnd(pattern)
.rightShift(start);
};
Filtering bad observations (III)
Filter bad observations based on StateQA Band (More info)
var filterBadObs = function(image){
// [bit 0-1] MOD35 cloud 0:clear 1:cloudy 2:mixed 3:not set, assumed clear
// [bit 2] Cloud shadow 0: no 1: yes
// [bit 15] Snow 0:no 1:yes
var cloudQA = getQABits(image.select('StateQA'),0,1,'cloud')
var cloudShadowQA = getQABits(image.select('StateQA'),2,2,'cloud_shadow')
var snow = getQABits(image.select('StateQA'),15,15,'snow')
var maskedImage = ee.Image(0).where(
cloudShadowQA.eq(0) // No cloud shadow
.and(cloudQA.eq(0).or(cloudQA.eq(3))) // Clear or unset
.and(snow.eq(0))// No snow
,1); //Put a 1 on good pixels
return image.mask(maskedImage) //Return only good pixels
}
Filtering bad observations (IV)
Test a single image
//---------------------------------------------------
// 1 Test a single image for bad observation filtering
//---------------------------------------------------
var sampleImg = ee.Image('MODIS/MOD09A1/MOD09A1_005_2010_01_09')
var sampleFilteredImage = filterBadObs(sampleImg)
Map.addLayer(sampleImg,{min:100,max:8000},'Original')
Map.addLayer(sampleFilteredImage,{min:100,max:8000},'Filtered')
Filtering bad observations (V)
Test image collection
//---------------------------------------------------
// 2 Test a complete collection for bad observation filtering
//---------------------------------------------------
var year = 2010
var col = ee.ImageCollection('MODIS/MOD09A1')
.filterDate(year+'-01-01',(year+1)+'-01-01')
var filteredCol = col.map(filterBadObs)
Map.addLayer(col.select(['sur_refl_b01','sur_refl_b04','sur_refl_b03']).mean(),
{min:100,max:8000},year+' Original collection')
Map.addLayer(filteredCol.select(['sur_refl_b01','sur_refl_b04','sur_refl_b03']).mean(),
{min:100,max:8000},year+ ' Filtered collection')
See the complete script results here
ee.Collection.iterate()
Sometimes, we just need to apply an algorithm to a collection IN ORDER
Arguments:
- algorithm (Function):The function to apply to each element. Must take two arguments: an element of the collection and the value from the previous iteration.
- opt_first (Object, optional):The initial state.
Examples:
- Gap filling algorithm
- LST > T for k consecutive days
ee.Collection.iterate()
function getMax(image,dict){
var currentMaxImage = ee.Image(ee.Dictionary(dict).get('currentMax'))
return ee.Dictionary({
currentMax: currentMaxImage.max(image)
})
}
var collection = ee.ImageCollection('MODIS/MOD11A2')
.filterDate('2010-01-01','2011-01-01')
.select(['LST_Night_1km'])
.sort('system:time_start')
var iniDict = ee.Dictionary({
currentMax : ee.Image(Math.NEGATIVE_INFINITY)
})
var result = ee.Dictionary(collection.iterate(getMax,iniDict))
.get('currentMax')
print(result)
How to get the maximum value of a collection at the pixel level withouth using ee.ImageCollection.max()?
LST > T for K consecutive days (I)
1 - Create a single pixel solution
Make groups of 2 or 3.
You have 5 minutes to think about it!
LST > T for K consecutive days (II)
-
// Initialization:
- t=0,k=3,contDaysGT0=0,firstDay=-1,fillValue=0;
-
// Iterate on every image in order
-
for (image in sorted(imageCollection)) do:
-
if contDaysGT0<k and LST>T and LST!=fillValue then:
- if (contDaysGT0==0) then firstDay=doy(image);
- contDaysGT0+=1;
-
else if contDaysGT0<k and (LST<=T or LST=fillValue )then:
- contDaysGT0 = 0;
-
if contDaysGT0<k and LST>T and LST!=fillValue then:
-
for (image in sorted(imageCollection)) do:
-
//Return results
- return firstDay
1 - Create a single pixel solution
LST > T for K consecutive days (III)
2 - Adapt solution to raster (2d)
It is exactly the same approach but we need to change variable types to be in 2d too!
A float / integer variable becomes float / integer Image
LST > T for K consecutive days (IV)
function getFirstDay(k,temp){
return function getFirstDay(image,dict){
//Retrieve input parameters
dict = ee.Dictionary(dict)
var cont = ee.Image(dict.get('cont'));
var firstDoY = ee.Image(dict.get('firstDoY'));
var fillValue = 0
var doy = ee.Number(ee.Date(image.get('system:time_start')).getRelative('day','year')).add(1)
//1: Get pixels that meet [contDaysGT0<k] and [LST>T] and [LST!=fillValue]
var positiveCondition = image.select([0]).multiply(0).where(
image.gte(temp)
.and(image.neq(fillValue))
.and(cont.neq(k))
,1)
//1.a: Set to current doy where positiveCondition==1 and cont==0
firstDoY = firstDoY.where(positiveCondition.and(cont.eq(0)),doy)
//1.b Increment cont for obsGtY pixels
cont = cont.add(positiveCondition)
//2: Restart counter on pixels where [contDaysGT0<k] and ([LST<=T] or [LST=fillValue] )
// It is the same as multiplying positiveCondition to counter, making sure we do not
// touch pixels that were found, this is, cont==k
var negativeCondition = positiveCondition.or(cont.eq(k))
cont = cont.multiply(negativeCondition)
//Return values for next iteration
var newDict = ee.Dictionary({k:k,cont:cont,firstDoY:firstDoY})
return newDict;
};
}
var initialParams = ee.Dictionary({
cont:ee.Image(0),
firstDoY:ee.Image(-1),
})
LST > T for K consecutive days (V)
Click here to see to entire script!
Google Earth Engine for Dummies (III)
By Michael a Menarguez
Google Earth Engine for Dummies (III)
- 6,630