Earth Engine

for Dummies

Part III

Michael A Menarguez

mamenarguez@ou.edu

Contents

  1. Paddy rice mapping algorithm solution (1st part)
  2. Filtering bad observations (where + mask + bit operations)
  3. Powerful ee.Image.iterate for algorithms!
  4. Guided exercise: LST > T for k days
  5. Download Results
    1. ​Image
    2. ImageCollection
    3. 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

 

  1. Test algorithm on a single image
  2. 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:

  1. 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.
  2. opt_first (Object, optional):The initial state.

Examples:

  1. Gap filling algorithm
  2. 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)

  1. // Initialization:
    1. t=0,k=3,contDaysGT0=0,firstDay=-1,fillValue=0;
  2. // Iterate on every image in order
    1. for (image in sorted(imageCollection)) do:
      1. if contDaysGT0<k and LST>T and LST!=fillValue then:
        1. if (contDaysGT0==0) then firstDay=doy(image);
        2. contDaysGT0+=1;
      2. else if contDaysGT0<k and (LST<=T or LST=fillValue )then:
        1. contDaysGT0 = 0;
  3. //Return results
    1. 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)

  • 3,286
Loading comments...

More from Michael a Menarguez