Multi-temporal NDVI covariate stacking for species distribution models
Satellite-derived NDVI phenology, land-surface temperature and seasonal vegetation statistics form the environmental covariate layers that feed species distribution models. Choosing the wrong temporal window or spatial grain quietly degrades transferability and inflates apparent accuracy.
Sensors
- MODIS MOD13Q1 (Terra/Aqua): 250 m NDVI and EVI composited at 16-day intervals. The long archive from 2000 onwards makes it the standard source for multi-year phenology statistics: growing-season length, peak NDVI, green-up timing. Coarse resolution limits use to species with broad habitat associations or large range extents.
- MODIS MOD11A2 Land Surface Temperature: 1 km, 8-day daytime and night-time LST composites. Thermal seasonality, mean cold-quarter temperature and frost-day frequency are among the most ecologically informative covariates for ectotherm and plant distribution models. Cloud-contaminated pixels are flagged but gaps remain a persistent issue in humid tropics.
- Landsat OLI (Landsat 8 and 9): 30 m multispectral imagery with a 16-day revisit per satellite, halved to roughly 8 days when both are combined. Enables NDVI, NDWI and bare-soil indices at a spatial grain appropriate for habitat patches of a few hectares. Archive depth to 1984 supports long-term phenological baselines.
- Sentinel-2 MSI (ESA): 10 m (visible/NIR) and 20 m (red-edge, SWIR) bands, 5-day revisit at the equator with both satellites. Red-edge bands (B5, B6, B7) add sensitivity to canopy chlorophyll content beyond what broadband NDVI captures. Shorter archive (from 2015) compared to Landsat, but finer grain suits narrow habitat features.
What the stack actually contains
A covariate stack for a species distribution model is not a single image. It is a set of summary statistics computed across time: mean NDVI for the wettest quarter, standard deviation of NDVI across a calendar year, the date of maximum greenness, amplitude of the seasonal NDVI cycle, mean and minimum land-surface temperature for the coldest 8-day composite period. Each statistic is a separate raster layer, and the model treats each pixel as a point in that multi-dimensional environmental space.
The practical consequence is that the quality of the model depends less on any single image than on the consistency of the temporal aggregation. A 16-year mean of MOD13Q1 peak NDVI is a stable, reproducible covariate. A single-year Sentinel-2 NDVI from an anomalously dry season is not, and projecting a model trained on it to a different year or region will produce quietly wrong predictions.
Temporal window choice is a modelling decision, not a data-processing detail
The aggregation window should match the organism's ecology. A migratory bird using a site only during the boreal spring requires spring-phenology covariates, not annual means. An amphibian constrained by breeding-pond availability needs wet-season surface-moisture indices. Averaging across the full year dilutes exactly the signal that separates suitable from unsuitable habitat.
Standard practice borrows from the BIOCLIM variable set: 19 bioclimatic variables derived from temperature and precipitation, several of which can be approximated from LST and NDVI-based moisture proxies. But satellite-derived phenology metrics can go further. Green-up day-of-year, calculated from MODIS or Landsat time series by fitting a double-logistic curve to the annual NDVI trajectory, captures the timing of resource availability in a way that quarterly means cannot. The published literature on this approach is extensive; the method is not experimental.
One honest caveat: phenology fitting requires at least 10 to 15 cloud-free observations per year at each pixel. In persistently cloudy landscapes, such as montane cloud forest or the humid tropics, annual NDVI time series have large gaps. Gap-filling by harmonic regression or Savitzky-Golay smoothing introduces assumptions that propagate into phenology metrics and, downstream, into model predictions. Analysts should report the proportion of gap-filled observations in each covariate layer.
Spatial resolution and the mismatch problem
Species distribution models are sensitive to the grain of their covariates. A 1 km LST pixel covers a wide range of micro-habitats in hilly terrain. A 10 m Sentinel-2 pixel may resolve a single tree crown. Neither is universally better; the right choice depends on the species' home-range size and the spatial scale of the environmental gradient being modelled.
Mixing resolutions within a single stack, for example combining 250 m MODIS NDVI with 30 m Landsat LST, is common practice but introduces spatial autocorrelation artefacts at the coarser pixel boundaries. Resampling everything to a common grain before modelling is necessary; bilinear resampling for continuous covariates and nearest-neighbour for categorical layers. The choice of output grain also sets a hard floor on the spatial resolution of predictions: a MaxEnt model trained on 1 km covariates produces 1 km probability surfaces, regardless of how fine the occurrence point data are.
For small-island or narrow-habitat contexts, such as riparian corridors or coastal ecotones, the 250 m MODIS grain is simply too coarse to resolve meaningful habitat variation. Landsat or Sentinel-2 derived phenology is necessary, which means accepting a shorter archive or greater cloud-gap risk depending on the site.
What the model predicts, and what it does not
This point is not a disclaimer buried at the end. It belongs at the beginning of any briefing to a decision-maker. A species distribution model trained on NDVI covariates predicts the environmental envelope within which recorded occurrences have been found. It does not confirm that the species is present in any predicted pixel. It does not account for dispersal barriers, interspecific competition, hunting pressure or land tenure. A pixel with a high habitat-suitability score may be physically inaccessible to the species, or may have been colonised and then abandoned.
Model validation requires independent field occurrence records, meaning records not used in model training, ideally collected using a systematic survey design rather than opportunistic observations. Metrics such as AUC (area under the receiver operating characteristic curve) are frequently over-interpreted. An AUC above 0.9 on training data can coexist with poor spatial transferability when the model is projected to a new region or a future climate scenario. Reporting partial AUC, continuous Boyce index or true skill statistic alongside spatial cross-validation results gives a more honest picture of predictive performance.
Building the stack in practice
For a continental or national-scale model, MOD13Q1 NDVI and MOD11A2 LST provide the most consistent multi-decadal baseline. Google Earth Engine's public catalogue holds the full MODIS archive and allows per-pixel temporal statistics to be computed at scale without local data storage. For finer-grain national assessments, Sentinel-2 Level-2A surface reflectance products, available through the Copernicus Data Space Ecosystem, support NDVI, red-edge chlorophyll index and SWIR-based moisture covariates at 10 to 20 m.
Satellize runs covariate-stack production on open constellations, including Sentinel-2 and Landsat, as part of its analytics services. The Tonga crop-estimation programme demonstrated that phenology-based time-series analysis on small island territories is feasible even with limited local cloud-computing infrastructure, which is a relevant constraint for many biodiversity programmes in island or low-income-country contexts.
Covariate layers should be delivered with accompanying uncertainty metadata: number of valid (non-gap-filled) observations per pixel, coefficient of variation across the aggregation period, and, where applicable, the root-mean-square error of the phenology curve fit. Without these, downstream modellers have no way to weight pixels by data quality or flag predictions built on thin temporal evidence.
Where the method reaches its limits
NDVI-based covariates are blind to below-canopy habitat structure. A closed-canopy forest pixel and a structurally degraded forest with the same leaf-area index may return near-identical NDVI values while supporting entirely different understorey communities. Spaceborne lidar (GEDI, ICESat-2) adds canopy height and vertical structure metrics that partially address this, but GEDI's 25 m footprint and non-contiguous sampling mean it supplements rather than replaces spectral covariates.
Aquatic and subterranean species are largely invisible to optical vegetation indices. Soil-moisture proxies from SWIR bands and LST can indicate wetland suitability, but the method's primary domain is terrestrial vegetation-associated fauna and flora. For species whose distributions are structured by geology, soil chemistry or subsurface hydrology, NDVI stacking adds limited explanatory power without complementary geophysical layers.
Typical figures
| Typical spatial resolution (NDVI covariates) | 250 m (MODIS MOD13Q1); 30 m (Landsat OLI); 10–20 m (Sentinel-2 MSI) |
| Typical spatial resolution (LST covariates) | 1 km (MODIS MOD11A2); no direct Sentinel-2 or Landsat thermal at comparable revisit for LST time series |
| Temporal compositing interval | 16-day (MODIS NDVI); 8-day (MODIS LST); 16-day per Landsat satellite, ~8-day combined; 5-day (Sentinel-2 pair) |
| Archive depth | MODIS: 2000–present; Landsat: 1984–present; Sentinel-2: 2015–present |
| Cloud sensitivity | All optical sensors. Persistent cloud cover in humid tropics or montane zones can reduce valid annual observations to fewer than 10 per pixel, requiring gap-filling |
| Phenology metric precision | Green-up day-of-year uncertainty typically ±5–15 days depending on gap-fill fraction and curve-fitting method; higher uncertainty in cloud-affected pixels |
| Minimum habitat patch resolvable | ~1 ha at 30 m Landsat; ~0.25 ha at 10 m Sentinel-2; ~6 ha at 250 m MODIS (rule of thumb: patch must span several pixels) |
| Delivery formats | GeoTIFF raster stack per covariate; NetCDF for multi-temporal arrays; CSV occurrence-covariate extract for direct model input; QGIS/ArcGIS compatible |
Analytics Satellize can run
| Annual phenology covariate stack | Double-logistic or Savitzky-Golay curve fitting to MODIS or Landsat NDVI time series; per-pixel extraction of green-up date, peak NDVI, growing-season length and NDVI amplitude | Multi-band GeoTIFF with one layer per phenology metric, plus per-pixel observation-count and gap-fill-fraction metadata layers |
| Seasonal NDVI summary statistics | Temporal aggregation of Sentinel-2 or Landsat surface-reflectance NDVI by ecologically defined quarter (wettest, driest, warmest, coldest) following BIOCLIM conventions | GeoTIFF stack of quarterly mean, maximum and standard-deviation NDVI layers, ready for direct import into MaxEnt or BRT modelling workflows |
| Land-surface temperature covariate layers | MODIS MOD11A2 8-day LST composites aggregated to mean cold-quarter, mean warm-quarter and annual temperature range; resampled to match NDVI stack grain | GeoTIFF LST covariate layers with valid-pixel count per aggregation period; flagged where fewer than four valid composites contribute to a quarterly mean |
| Occurrence-covariate extract table | Spatial join of species occurrence records (GBIF or client-supplied) to covariate stack at each occurrence point; duplicates thinned to one record per pixel to reduce spatial autocorrelation | CSV table of occurrence coordinates with extracted covariate values per point, formatted for direct input to MaxEnt, R-dismo or biomod2 |
| Model transferability uncertainty assessment | Spatial block cross-validation (checkerboard or geographic k-fold) to estimate predictive performance under geographic extrapolation; MESS (Multivariate Environmental Similarity Surface) mapping to flag extrapolation zones | PDF report with spatial cross-validation AUC and TSS distributions, MESS raster identifying pixels outside the training environmental envelope |
| Multi-year covariate change detection | Comparison of phenology statistics between two user-defined epochs (e.g. 2000–2010 vs 2013–2023) to identify pixels where the environmental envelope has shifted, flagging areas where historical model predictions may no longer be valid | Change-magnitude GeoTIFF per covariate layer; summary table of mean shift per habitat class |
Who does the work
We can get this done for you. Satellize runs its own analyst desk and a strong science team. You do not buy a data feed and work out what it means; our people source the imagery, run the analysis described on this page, and hand you the answer with its confidence limits stated. Discuss this requirement.