2

This question refers to Calculating NDVI per region, month & year with Google Earth Engine

I have modified the code posted by @Kel Markert at https://code.earthengine.google.com/349615d7802d59f677181bef0badad9f to attempt to get a maximum monthly NDVI value from a small polygon over a number of years from Landsat 8 in Google Earth Engine and export to CSV.

The output .csv table only gives -999 values for all months when using the max reducer reducer: ee.Reducer.max(). When using the median reducer reducer: ee.Reducer.median(), the output provides meaningful values.

Is Max the incorrect reducer to use?

Code Link to Google Earth Engine Script

Table asset link - region polygon

    var region = table,
L8 = ee.ImageCollection("LANDSAT/LC08/C01/T1_TOA");

var cloudlessNDVI = L8.map(function(image) {
  // Get a cloud score in [0, 100].
  var cloud = ee.Algorithms.Landsat.simpleCloudScore(image).select('cloud');

  // Create a mask of cloudy pixels from an arbitrary threshold.
  var mask = cloud.lte(20);

  // Compute NDVI.
  var ndvi = image.normalizedDifference(['B5', 'B4']).rename('NDVI');

  // Return the masked image with an NDVI band.
  return image.addBands(ndvi).updateMask(mask);
});

var startDate = ee.Date('2013-05-01'); // set analysis start time
var endDate = ee.Date('2019-12-31'); // set analysis end time

// calculate the number of months to process
var nMonths = ee.Number(endDate.difference(startDate,'month')).round();

// get a list of time strings to pass into a dictionary later on
var monList = ee.List.sequence(0, nMonths).map(function (n) {
  return startDate.advance(n, 'month').format('YYYMMdd');
})
print(monList)

var result = region.map(function(feature){
  // map over each month
  var timeSeries = ee.List.sequence(0,nMonths).map(function (n){
    // calculate the offset from startDate
    var ini = startDate.advance(n,'month');
    // advance just one month
    var end = ini.advance(1,'month');
    // filter and reduce
    var data = cloudlessNDVI.filterDate(ini,end).max().reduceRegion({
      reducer: ee.Reducer.max(),
      geometry: feature.geometry(),
      scale: 30
    });
    // get the value and check that it has data
    var val = ee.Number(data.get('NDVI'));
    val = ee.Number(ee.Algorithms.If(val,val,-999));
    // return max
    return val;
  });
  // create new dictionary with date strings and values
  var timeDict = ee.Dictionary.fromLists(monList,timeSeries);
  // return feature with a timeseries property and results
  return feature.set(timeDict);
});

// print to see if it is doing what we expect...
print(result);

// Export the data to a table for further analysis
Export.table.toDrive({
  collection:result,
  description:"MCM1_NDVI",
  fileFormat:"CSV",
  //selectors:["HRpcode","timeseries"]
})
2

You can obtain the answer you want with the small geometry you use using a different approach. Note that a Landsat pixel and the scale you are currently using is 30m. Your geometry is just about 10x10m.

Add the moment you are using an unweighted reducer (ee.Reducer.max), which will only take pixels into account if the centroid is within the geometry (see here). Neither the geometry nor the scale you are using apparently include the centroid of the pixels, thus you return no valid values (which you set to -999 afterwards).

You have two options:

1) Use a scale which will include the centroid and use the unweighted reducer. Set the scale for example at 10m or calculate approximately the scale of a square geometry using:

var scale = table.geometry().area(1).sqrt();

2) Use a weighted reducer which will include pixels which approximately include at least 0.5% of the pixel. An example would be using a mean reducer. In such a small region you probably aggregate one pixel, so that is fine for your max outcome.

var data = cloudlessNDVI.filterDate(ini,end).max().reduceRegion({
  reducer: ee.Reducer.mean(),
  geometry: feature.geometry(),
  scale: 30
});

Here an example returning both values (and you will see they are similar, confirming there is just one pixel in the region).

| improve this answer | |
  • Fantastic. Thanks @Kuik. Very clearly answered. Much appreciated. – Jock Jan 30 at 0:26
1

I think your region is too small. I've tried running this code with a full county in Hawaii, and it worked fine:

var region = ee.FeatureCollection('TIGER/2016/Counties')
  .filter(ee.Filter.eq('NAME', 'Kalawao'));

Additionally, I created a polygon on the map centered on your table object but about 3x the size, and the code ran fine with -999 results occurring only rarely:

var geometry = ee.Geometry.Polygon(
  [[149.18765874159647,-21.07471411401736],
   [149.1883158828051,-21.07471411401736],
   [149.1883158828051,-21.074181016177615],
   [149.18765874159647,-21.074181016177615],
   [149.18765874159647,-21.07471411401736]]);

var region = ee.FeatureCollection(geometry);

Substituting either of these regions for your original table work as expected.

| improve this answer | |
  • Many thanks @JepsonNomad. Weird that it works with Median but not Max? I would've though that 100m2 would be suitable for Landsat eg. contain enough pixels? – Jock Jan 29 at 23:22
  • 1
    A landsat pixel should be about 900m2 if it's 30m on a side - so if you're just sampling one pixel at each timepoint, maybe EarthEngine will give you the median (which is just that pixel's value) but doesn't like the idea of calling max? That's just a guess though... – JepsonNomad Jan 29 at 23:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.