Google Earth Engine Script for Extracting MODIS EVI Seasonal and Raw Time-Series Zonal Statistics

This repository provides a set of Google Earth Engine (GEE) scripts developed to extract Enhanced Vegetation Index (EVI) statistics from the MODIS MOD13Q1 (Version 6.1) dataset over Malawi. The scripts generate either (1) seasonal composites for the November–April vegetation period or (2) raw 16-day composite EVI values without temporal aggregation. They also implement zonal statistics for buffered regions of interest (ROIs) at 5 km and 10 km distances, with automated export of results as CSV tables to Google Drive.

The code is designed for reproducible vegetation and climate studies, including ecosystem productivity monitoring, land-cover phenology analysis, and climate–vegetation interaction assessments. It demonstrates best practices for temporal filtering, scale factor correction, buffer-based zonal reduction, and metadata tagging (year, month, buffer distance, image timestamp).

Data source

• MODIS Terra Vegetation Indices 16-Day (L3, 250 m): MODIS/061/MOD13Q1
  Provided via Google Earth Engine (https://earthengine.google.…).

Programming environment

Google Earth Engine JavaScript API.

How to cite

If you use this code or adapt it in your research, please cite as follows:

Ramdani, F. (2025). Google Earth Engine Script for Extracting MODIS EVI Seasonal and Raw Time-Series Zonal Statistics [Code]. Zenodo. [DOI: 10.5281/zenodo.17528257]

Keywords

MODIS, EVI, vegetation index, remote sensing, Google Earth Engine, Malawi, phenology, climate–vegetation interaction, reproducible research

License

Creative Commons Attribution 4.0 International (CC BY 4.0)

Critical points you should scrutinize if you want the results to stand up under peer review or publication scrutiny:

1. Radiometric scaling of MODIS EVI values

The MODIS EVI band in MOD13Q1 is stored as a 16-bit integer scaled by 0.0001.
If you use the raw values directly (without multiplying by 0.0001), your EVI range will appear as hundreds to thousands instead of −1 to +1, which renders all derived statistics physically meaningless.
Always rescale before any spatial or temporal aggregation.
Failing to do so is a fatal methodological error, not a minor oversight.

2. Definition of “raw” imagery

The code current logic treats each MODIS scene as “raw.”
However, MOD13Q1 is already a 16-day composite product, not a daily observation.
So, “raw” here means “unaggregated 16-day composites,” not truly raw radiances.
You should explicitly state this distinction, or reviewers will accuse you of misusing the term “raw.”
If you truly need uncomposited daily or 8-day values, you’d have to switch to another MODIS or VIIRS product.

3. Temporal sampling within the Nov–Apr season

Because each MODIS EVI scene represents a 16-day period, in this code Nov–Apr filter actually includes roughly nine composites per season (not six monthly points).
If you intend to model seasonal vegetation dynamics, make sure your temporal spacing is correctly interpreted — e.g., whether you treat each composite as independent or interpolate them to monthly medians.
Reviewers will immediately flag temporal aliasing if your analysis assumes monthly frequency but uses 16-day composites.

4. Spatial resolution and reduction scale

Using .reduceRegions() with a 5000 m scale on a 250 m MODIS grid averages roughly 400 pixels per sample.
That is acceptable for broad regional signals but can obscure spatial heterogeneity, especially for 5 km buffers.
To ensure your results are not artefacts of oversmoothing, you should test sensitivity to the scale parameter (e.g., 1000 m vs. 5000 m).
If your study area is where vegetation gradients are steep, 5 km resolution might erase meaningful contrasts.

5. Metadata traceability and reproducibility

You are exporting time-series data but without including the MODIS scene ID or acquisition date (only system_time_start).
That timestamp is sufficient for internal consistency but not for unambiguous cross-reference to MODIS tiles or quality flags.
Add at least the system:index property so each row can be linked back to the original image.
Without that, reviewers cannot verify which composites contributed to each extracted point.

6. Interpretation of EVI magnitudes

Once scaled, EVI values rarely exceed 0.8 even in dense vegetation.
If your final CSV still shows values >1.0 after scaling, it indicates inclusion of fill values, snow, or cloud contamination.
You should either apply the SummaryQA mask from the MOD13Q1 product or filter out anomalous pixels.
Unmasked artifacts will undermine any temporal trend analysis.