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
}).
combine
({
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
=
intervals
.
map
(
function
(
i
)
{
// Calculate temporal composite.
var
startRangeL
=
ee
.
Date
(
startDate
).
advance
(
i
,
intervalUnit
);
var
endRangeL
=
startRangeL
.
advance
(
interval
,
intervalUnit
);
var
temporalStat
=
dataset
.
filterDate
(
startRangeL
,
endRangeL
)
.
reduce
(
temporalReducer
);
// 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
statsL
.
map
(
function
(
feature
)
{
return
feature
.
set
({
composite_start
:
startRangeL
.
format
(
'YYYY'
),
// or 'YYYY-MM-dd'
});
});
});
zonalStatsL
=
ee
.
FeatureCollection
(
zonalStatsL
).
flatten
();
print
(
'Spatiotemporal statistics (long)'
,
zonalStatsL
);
Export
.
table
.
toDrive
({
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
=
keys
.
map
(
function
(
key
)
{
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
()
.
reduce
(
temporalReducer
);
zonalStatsW
=
zonalStatsW
.
map
(
reduce
);
}
print
(
'Spatiotemporal statistics (wide)'
,
zonalStatsW
);
Export
.
table
.
toDrive
({
collection
:
zonalStatsW
,
description
:
'zonal_stats_wide'
,
});
