Spatiotemporal Statistics of Vegetation Indices

Edit on GitHub
Report issue
Page history
Author(s): saumyatas
This tutorial demonstrates how to create a comma delimited table of zonal statistics of vegetation indices (NDVI or EVI) over a study area, for a given range of years.

Import datasets

Vegetation index

Google Earth Engine has a range data products that provide time series of vegetation indices. Here, we use the MODIS Terra Vegetation Indices for 16-days Global 250m product (also available at 500m and 1km resolution). After importing, we select the 'EVI' band.

var dataset = ee.ImageCollection('MODIS/006/MOD13Q1')
                  .select('EVI');  // or 'NDVI'

Region of interest

A FeatureCollection (or Geometry) is needed to define regions to summarize vegetation index data over. For example, you can use the Global Administrative Unit Layers (GAUL) dataset, to extract zonal statistics for administrative regions. The GAUL provides administrative unit boundaries for all countries in the world at three levels. Here, we'll use the districts of Maharashtra, India; we'll get zonal statistics for each district (35 districts).

var regions = ee.FeatureCollection('FAO/GAUL/2015/level2')
                  .filter(ee.Filter.inList('ADM1_NAME', ['Maharashtra']));

Alternatively, you can select a different FeatureCollection from the Earth Engine Data Catalog or upload your own ShapeFile by following the instructions in the Importing Table Data section of the Earth Engine Developer's Guide.

Defining spatiotemporal reduction parameters

The statistics we're after include spatial and temporal components. Above, we've defined the bounds of the spatial component, now we define the temporal component, i.e. the series of time windows to generate representative composite vegetation index images for. The following variables set the length a time window and the duration of the series.

The variables are set to generate zonal statistics (spatialReducers) for image composites constructed from three time periods (intervalCount) that start on 2018-01-01 (startDate) with each time window being 1 (interval) year (intervalUnit) reduced by mean (temporalReducer). In other words, we'll be generating annual mean vegetation index zonal statistics for years 2018, 2019, and 2020.

var startDate = '2018-01-01';  // time period of interest beginning date
var interval = 1;  // time window length
var intervalUnit = 'year';  // unit of time e.g. 'year', 'month', 'day'
var intervalCount = 3;  // number of time windows in the series
var temporalReducer = ee.Reducer.mean();  // how to reduce images in time window

// Defines mean, standard deviation, and variance as the zonal statistics.
var spatialReducers = ee.Reducer.mean().combine({
        reducer2: ee.Reducer.stdDev(),
        sharedInputs: true
        reducer2: ee.Reducer.variance(),
        sharedInputs: true

Calculating spatiotemporal statistics

Now, we'll calculate spatiotemporal vegetation index statistics. Two ways to arrange the statistics table are provided: a long table and a wide table. Depending on your application, one arrangement might be more suitable than the other. The resulting tables are exported to Google Drive as a CSV file, but there are multiple other ways to export, including downloading the CSV file directly using getDownloadURL.

Long table

The long table format will have one row per unique combination of region and time window. So in this example, there will be 105 rows (35 districts * 3 time windows). First, a list of 0-based time window indices is constructed, then we map over the time window index list to generate composite images for each time window, and then apply reduceRegions to calculate the zonal statistics per region. Finally, the start date of the time window is added as a feature property so the statistics can be tied to a given image composite. The result is a FeatureCollection of FeatureCollections, which must be flattened. The flattened FeatureCollection is then exported to Google Drive as a CSV file.

// Get time window index sequence.
var intervals = ee.List.sequence(0, intervalCount-1, interval);

// Map reductions over index sequence to calculate statistics for each interval.
var zonalStatsL = {
  // Calculate temporal composite.
  var startRangeL = ee.Date(startDate).advance(i, intervalUnit);
  var endRangeL = startRangeL.advance(interval, intervalUnit);
  var temporalStat = dataset.filterDate(startRangeL, endRangeL)

  // Calculate zonal statistics.
  var statsL = temporalStat.reduceRegions({
    collection: regions,
    reducer: spatialReducers,
    scale: dataset.first().projection().nominalScale(),
    crs: dataset.first().projection()

  // Set start date as a feature property.
  return {
    return feature.set({
      composite_start: startRangeL.format('YYYY'),  // or 'YYYY-MM-dd'

zonalStatsL = ee.FeatureCollection(zonalStatsL).flatten();

print('Spatiotemporal statistics (long)', zonalStatsL);

  collection: zonalStatsL,
  description: 'zonal_stats_long',

Wide table

The wide table will have one row for each region (35 rows in this case) with a column per unique combination of statistic and time window. The new column names are the concatenation of the statistic and the time window separated by an underscore. The process uses a client-side for loop on the number of time windows, and for each one, zonal statistics are calculated and appended as new properties to the FeatureCollection. An alternative approach is to use Earth Engine's iterate function, but a comparison of the approaches' performance was equal, so we've chosen the for loop for improved readability.

var temporalStatW;  // defined dynamically for each temporal image composite
var startRangeW;  // defined dynamically for each temporal image composite
var reduce = function(feature) {
  // Calculate zonal statistics.
  var statsW = temporalStatW.reduceRegion({
    reducer: spatialReducers,
    geometry: feature.geometry(),
    scale: dataset.first().projection().nominalScale(),
    crs: dataset.first().projection()

  // Append date to the statistic label.
  var keys = ee.Dictionary(statsW).keys();
  var newKeys = {
    return ee.String(key).cat('_')
               .cat(startRangeW.format('YYYY'));  // or 'YYYY-MM-dd'

  // Add the statistic properties to the feature.
  return feature.set(statsW.rename(keys, newKeys));

var zonalStatsW = regions;  // make a copy of the regions FeatureCollection

// Loop through sequence of intervals to calculate statistics for each.
for (var i = 0; i < intervalCount; i++) {
  var startRangeW = ee.Date(startDate).advance(i, intervalUnit);
  var endRangeW = startRangeW.advance(interval, intervalUnit);
  temporalStatW = dataset.filterDate(startRangeW, endRangeW).mean()
  zonalStatsW =;

print('Spatiotemporal statistics (wide)', zonalStatsW);

  collection: zonalStatsW,