1

I want to look at NDVI and other spectral indices over time and compare these time series across different land cover types. I believe I have the code to do this but I'm running into the error that I have too many elements: Error generating chart: Collection query aborted after accumulating over 5000 elements. Would combining all those vectors into a single multi-part vector help? I would really like to do this over a much larger area so I assume I'll need a different approach?

// Function to cloud mask from the pixel_qa band of Landsat 8 SR data.
function maskL8sr(image) {
  // Bits 3 and 5 are cloud shadow and cloud, respectively.
  var cloudShadowBitMask = 1 << 3;
  var cloudsBitMask = 1 << 5;

  // Get the pixel QA band.
  var qa = image.select('pixel_qa');

  // Both flags should be set to zero, indicating clear conditions.
  var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0)
      .and(qa.bitwiseAnd(cloudsBitMask).eq(0));

  // Return the masked image, scaled to reflectance, without the QA bands.
  return image.updateMask(mask).divide(10000)
      .select("B[0-9]*")
      .copyProperties(image, ["system:time_start"]);
}

// fcn to add NDVI band to each image
var addNDVI = function(image) {
  var ndvi = image.normalizedDifference(['B5', 'B4']).rename('NDVI');
  return image.addBands(ndvi);
};

// Use the NLCD to filter areas by cover type and see how NDVI changed in those areas
// pull in NLCD
var NLCDdataset = ee.ImageCollection('USGS/NLCD');
var landcover = NLCDdataset.select('landcover');
var landcoverVis = {
  min: 11.0,
  max: 95.0,
  palette: [
    '5475A8', 'ffffff', 'E8D1D1', 'E29E8C', 'ff0000', 'B50000', 'D2CDC0',
    '85C77E', '38814E', 'D4E7B0', 'AF963C', 'DCCA8F', 'FDE9AA', 'D1D182',
    'A3CC51', '82BA9E', 'FBF65D', 'CA9146', 'C8E6F8', '64B3D5'
  ],
};
Map.addLayer(landcover, landcoverVis, 'Landcover');

// mask to your land cover class of interest
var landcover2011 = landcover.filter(ee.Filter.eq('system:index', 'NLCD2011')).first()

// so first we need a vector with all the wetland area

// Create a raster that is just the wetland class.
var wetland_raster = landcover2011.eq(95);
// mask the wetland raster to areas that are only wetland
wetland_raster = wetland_raster.updateMask(wetland_raster.neq(0));

// Convert the zones of the wetland areas to vectors.
var wetland_vectors = wetland_raster.addBands(landcover2011).reduceToVectors({
  geometry: roi, // this doesn't need to be just the ROI, could be much larger
  crs: wetland_raster.projection(),
  scale: 30,
  geometryType: 'polygon',
  eightConnected: false,
  labelProperty: 'zone',
  reducer: ee.Reducer.mean()
});

print(wetland_vectors);
//print(ee.Geometry.MultiPolygon(wetland_vectors))

// Make a display image for the vectors, add it to the map.
var display = ee.Image(0).updateMask(0).paint(wetland_vectors, '000000', 3);
Map.addLayer(display, {palette: '000000'}, 'vectors');


// now let's look at NDVI at this section over time

// Create image collection of l8 imagery
var l8 = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
//filter start and end date
.filterDate('2019-03-01', '2019-04-01')
//filter according to drawn boundary
.filterBounds(roi)
.map(maskL8sr);


// Add NDVI band to image collection
var l8 = l8.map(addNDVI);

// Create a time series chart.
var plotNDVI = ui.Chart.image.seriesByRegion(l8, wetland_vectors, ee.Reducer.mean(),
'NDVI',30,'system:time_start', 'system:index')
              .setChartType('LineChart').setOptions({
                title: 'NDVI time series',
                hAxis: {title: 'Date'},
                vAxis: {title: 'NDVI'}
});

print(plotNDVI);

Here is a link to this code: https://code.earthengine.google.com/64daced9e1c453f1fcd98180196e410d

1

To speed up this process, you will need to use grouped reducers. You can then skip the reduceToVectors which is computationally intensive and keep the reducers in the ee.Image space, which the Earth Engine is good at.

However, the output of a grouped reducer is a bit unhandy to rewrite to usefull properties. It is therefore that it takes some lines of code.

// get the number of times you want to use the reducer
var bands = l8.first().bandNames().length();
// get the values of the image landcover
var values = ee.List(landcover2011.get('landcover_class_values'));
var keysList = values.map(function(value){
  var groupName = ee.String('group_').cat(ee.String(value)).cat(ee.String('_'));
  var withBandNames = l8.first().bandNames().map(function(bandName){
    return groupName.cat(ee.String(bandName));
  });
  return withBandNames;
});

Now map over the collection and use the grouped reducer.

// Get statistics for every landcover type using a grouped reducer
var information = l8.map(function(image){
  image = image.addBands(landcover2011);
  // per landcover type statistics for all bands
  var reduced = ee.List(image.reduceRegion({
    reducer: ee.Reducer.mean().repeat(bands).group(bands), 
    scale: 30,
    geometry: roi,
    maxPixels: 10e12
  }).get('groups'));

  // values actually present inside the aoi of the landcover 
  //   (masked landcover types will be omitted in the reducer output)
  var actualVals = reduced.map(function(val){
    return ee.Dictionary(val).get('group');
  });

  // rewrite the grouped reducer output
  var rewritten = keysList.map(function(keys){
    keys = ee.List(keys);
    var group = ee.String(keys.get(0)).slice(6,8);
    var indexGroup = actualVals.indexOf(ee.Number.parse(group));
    var vals = ee.List(ee.Dictionary(reduced.get(indexGroup)).get('mean'));
    return vals;
  }).flatten();

  // write to a dictionary
  var dictionary = ee.Dictionary.fromLists(keysList.flatten(), rewritten);
  return ee.Image(image.setMulti(dictionary));
});

Now use regular expressions to get the properties you are interested in. As an example, I used NDVI.

// for example, get statistics for the NDVI
var NDVIstats = ee.FeatureCollection(information.map(function(feat){
  var dict = ee.Feature(feat).toDictionary().select([".*_NDVI.*"], true);
  return ee.Feature(null, dict.combine(
                 ee.Dictionary({'system:time_start': feat.get('system:time_start')})));
}));
print(NDVIstats)
print(ui.Chart.feature.byFeature(NDVIstats, 'system:time_start').setSeriesNames(legend).setOptions({
colors: palette, interpolateNulls: true
}));

In the link to the code you can also find the palette and names ('legend') variables I defined for making a nice looking graph: link

  • Wow this is incredibly helpful! Thanks so much. The noisiness and jumpiness of the data seems like it is probably not reflecting true on the ground changes. I realize this is a big question in and of itself but do you have any advice for filtering that data, ruling out unrealistic values, or just better understanding some jumps (e.g. the huge decline in NDVI in December of 2016). Thanks! – clifgray Apr 14 at 21:46
  • Yes that is inherent on time series analysis. Especially for NDVI, you should do some harmonic modelling such as in this example: docs.google.com/document/d/… – Kuik Apr 15 at 8:28

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.