Seeing Green from Space: Vegetation Health with Satellite Indices

Vegetation indices compress canopy optics into stable, comparable numbers. Healthy leaves absorb red (chlorophyll) and strongly scatter NIR (leaf structure). Stress reduces red absorption and/or NIR scattering. The goal isn’t a prettier map; it’s a reliable, repeatable signal you can threshold, alert on, and verify in the field.

Index formulas, variants, and when to use them

NDVI  = (NIR - Red) / (NIR + Red)
EVI   = G * (NIR - Red) / (NIR + C1*Red - C2*Blue + L)   // common: G=2.5, C1=6, C2=7.5, L=1
EVI2  = 2.5 * (NIR - Red) / (NIR + 2.4*Red + 1)          // when Blue is noisy/unavailable
SAVI  = (1 + L) * (NIR - Red) / (NIR + Red + L)          // L≈0.5; soil-adjusted, early season
MSAVI2= (2*NIR + 1 - sqrt((2*NIR + 1)^2 - 8*(NIR - Red))) / 2
CI_re = (NIR / RedEdge) - 1                              // Sentinel-2: RedEdge ∈ {B5,B6,B7}

Band mapping — explicit wavelengths

Quality control — what to mask, what to drop

• Use surface reflectance. S2 L2A, Landsat SR, HLS. Avoid TOA for indices if SR is available.
• Clouds & shadows. S2 cloud probability < 30% (s2cloudless) + shadow geometry; Landsat QA_PIXEL bits for cloud/shadow/snow; MODIS/VIIRS use VI quality layers.
• View geometry. Exclude high view-zenith (> 15–20°) and extreme solar angles; reduce BRDF artifacts.
• Terrain. On slopes, apply topographic correction (e.g., C-correction with DEM-derived illumination).
• Edges. Buffer parcel polygons inward 1–2 pixels to avoid mixed pixels and georegistration jitter.

Minimal, robust compositing

• Temporal: 10–16-day rolling window; median for stability; max-NDVI to counter residual haze.
• Spatial (per parcel): median or P75 across valid pixels for robustness; track pixel count (coverage).
• Phenology-aware: compare DOY-aligned windows year-over-year (e.g., DOY ± 15).

Anomaly math that scales

1) Baseline: For each parcel and day-of-year (DOY), compute multi-year median NDVI (rolling window DOY±15). 
2) Smooth current season (e.g., 3x compositing windows or Savitzky–Golay).
3) Anomaly_t    = NDVI_t - Baseline_DOY_median
   Zscore_t     = (NDVI_t - μ_baseline) / σ_baseline
4) Alert rule: anomaly < -0.15 OR z < -1.0 for 2 consecutive windows → field check.
   Dense canopies: also watch EVI drop ≥ 0.08 vs 4-week median.
   Red-edge early warning: CI_re drop ≥ 0.10 vs parcel DOY baseline.

Cross-sensor harmonization (S2 ↔ Landsat)

If you need cadence over clouds, use HLS (Harmonized Landsat–Sentinel). It applies bandpass adjustments and BRDF normalization so indices are comparable at 30 m. Keep your thresholds sensor-agnostic by anchoring to each parcel’s historical baseline rather than absolute values.

Google Earth Engine — Sentinel-2 NDVI/EVI/EVI2/CIre + compositing + parcel stats

var s2sr = ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED')
  .filterBounds(parcels.geometry())
  .filterDate('2025-03-01', '2025-10-31')
  .filter(ee.Filter.lte('CLOUDY_PIXEL_PERCENTAGE', 40))
  .map(function(img){
    var qa = img.select('QA60');
    var cloud = qa.bitwiseAnd(1<<10).neq(0);
    var cirrus = qa.bitwiseAnd(1<<11).neq(0);
    var mask = cloud.or(cirrus).not();

    var red = img.select('B4').multiply(1e-4);
    var nir = img.select('B8').multiply(1e-4);
    var blue= img.select('B2').multiply(1e-4);
    var re5 = img.select('B5').multiply(1e-4);
    var re6 = img.select('B6').multiply(1e-4);
    var re7 = img.select('B7').multiply(1e-4);

    var ndvi = nir.subtract(red).divide(nir.add(red)).rename('NDVI');
    var evi  = nir.subtract(red)
                 .multiply(2.5)
                 .divide(nir.add(red.multiply(6)).subtract(blue.multiply(7.5)).add(1))
                 .rename('EVI');
    var evi2 = nir.subtract(red)
                  .multiply(2.5)
                  .divide(nir.add(red.multiply(2.4)).add(1))
                  .rename('EVI2');
    var cire = nir.divide(re5).subtract(1).rename('CIre');

    return img.updateMask(mask).addBands([ndvi, evi, evi2, cire])
              .set('date', img.date().format('YYYY-MM-dd'));
  });

var step = 12;
var start = ee.Date('2025-03-01');
var end   = ee.Date('2025-10-31');

var windows = ee.List.sequence(0, end.difference(start,'day').subtract(step), step)
  .map(function(d){
    var winStart = start.advance(ee.Number(d), 'day');
    var winEnd   = winStart.advance(step, 'day');
    var comp = s2sr.filterDate(winStart, winEnd)
      .select(['NDVI','EVI','EVI2','CIre'])
      .median()
      .set('window_start', winStart.format('YYYY-MM-dd'))
      .set('window_end',   winEnd.format('YYYY-MM-dd'));
    return comp;
  });

var compCol = ee.ImageCollection.fromImages(windows);

var reducer = ee.Reducer.median().combine({
  reducer2: ee.Reducer.percentile([75]),
  sharedInputs: true
});

var zonal = compCol.map(function(img){
  var stats = img.reduceRegions({
    collection: parcels,
    reducer: reducer,
    scale: 10,
    tileScale: 4
  }).map(function(f){
    return f.set({
      window_start: img.get('window_start'),
      window_end:   img.get('window_end')
    });
  });
  return stats;
}).flatten();

Export.table.toDrive({
  collection: zonal,
  description: 'parcel_vi_summaries_2025',
  fileFormat: 'CSV'
});

Python — rioxarray + GeoPandas, parcel-level indices & 14-day medians

import xarray as xr, rioxarray as rxr, geopandas as gpd
import numpy as np, pandas as pd
from pathlib import Path

cogs_dir = Path("/data/S2_SR")
dates = sorted({p.name.split("_")[1] for p in cogs_dir.glob("S2_*_B4.tif")})
parcels = gpd.read_file("/data/parcels.gpkg").to_crs(32633)

def stats(a):
    d = a.data.compressed() if np.ma.isMaskedArray(a.data) else a.data
    d = d[np.isfinite(d)]
    if d.size == 0: return np.nan, np.nan, 0
    return np.median(d), np.percentile(d, 75), d.size

rows = []
for d in dates:
    red  = rxr.open_rasterio(cogs_dir/f"S2_{d}_B4.tif", masked=True).squeeze()
    nir  = rxr.open_rasterio(cogs_dir/f"S2_{d}_B8.tif", masked=True).squeeze()
    blue = rxr.open_rasterio(cogs_dir/f"S2_{d}_B2.tif", masked=True).squeeze()
    re5  = rxr.open_rasterio(cogs_dir/f"S2_{d}_B5.tif", masked=True).squeeze()

    ndvi = (nir - red) / (nir + red + 1e-6)
    evi  = 2.5 * (nir - red) / (nir + 6*red - 7.5*blue + 1)
    evi2 = 2.5 * (nir - red) / (nir + 2.4*red + 1)
    cire = (nir / (re5 + 1e-6)) - 1

    for _, poly in parcels.iterrows():
        geom = [poly.geometry]
        ndvi_med, ndvi_p75, n = stats(ndvi.rio.clip(geom, ndvi.rio.crs, drop=False))
        evi_med,  evi_p75,  _ = stats(evi.rio.clip(geom, evi.rio.crs, drop=False))
        evi2_med, evi2_p75, _ = stats(evi2.rio.clip(geom, evi2.rio.crs, drop=False))
        cire_med, cire_p75, _ = stats(cire.rio.clip(geom, cire.rio.crs, drop=False))

        rows.append({
          "date": pd.to_datetime(d),
          "parcel_id": poly["parcel_id"],
          "ndvi_med": ndvi_med, "ndvi_p75": ndvi_p75, "px_valid": int(n),
          "evi_med": evi_med, "evi_p75": evi_p75,
          "evi2_med": evi2_med, "evi2_p75": evi2_p75,
          "cire_med": cire_med, "cire_p75": cire_p75
        })

df = pd.DataFrame(rows).sort_values(["parcel_id","date"])

out = []
for pid, g in df.groupby("parcel_id"):
    g = g.set_index("date").sort_index()
    g_roll = g.rolling("14D").median()
    g_roll["parcel_id"] = pid
    out.append(g_roll.reset_index())

df_roll = pd.concat(out).sort_values(["parcel_id","date"])
df_roll.to_csv("/data/parcel_vi_14d.csv", index=False)

COGs fast path — GDAL to create web-ready tiles

gdal_translate input.tif output_cog.tif -of COG -co COMPRESS=DEFLATE -co NUM_THREADS=ALL_CPUS
gdalbuildvrt -separate stack.vrt NIR.tif RED.tif
gdal_calc.py -A NIR.tif -B RED.tif --calc="(A-B)/(A+B+1e-6)" --outfile NDVI_cog.tif --of COG --co COMPRESS=DEFLATE

PostgreSQL — parcel baselines, anomalies, and alert flags

WITH hist AS (
  SELECT parcel_id,
         EXTRACT(doy FROM date)::int AS doy,
         ndvi_med
  FROM parcel_vi_14d
  WHERE date < DATE '2025-01-01'
),
baseline AS (
  SELECT h1.parcel_id, h1.doy AS doy_ref,
         PERCENTILE_CONT(0.5) WITHIN GROUP (ORDER BY h2.ndvi_med) AS ndvi_baseline
  FROM hist h1
  JOIN hist h2
    ON h1.parcel_id = h2.parcel_id
   AND h2.doy BETWEEN h1.doy - 15 AND h1.doy + 15
  GROUP BY h1.parcel_id, h1.doy
)
SELECT cur.parcel_id,
       cur.date,
       cur.ndvi_med,
       b.ndvi_baseline,
       (cur.ndvi_med - b.ndvi_baseline) AS ndvi_anom,
       (cur.ndvi_med - b.ndvi_baseline) < -0.15 AS ndvi_alert
FROM parcel_vi_14d cur
JOIN baseline b
  ON b.parcel_id = cur.parcel_id
 AND b.doy_ref   = EXTRACT(doy FROM cur.date)::int
WHERE cur.date BETWEEN DATE '2025-03-01' AND DATE '2025-10-31'
ORDER BY cur.parcel_id, cur.date;

Heuristics (production defaults; tune locally)

• NDVI anomaly < −0.15 for two consecutive 10–16-day windows ⇒ task a field check.
• EVI drop ≥ 0.08 vs 4-week median in dense canopy ⇒ canopy health change likely.
• CI_re drop ≥ 0.10 vs parcel DOY baseline ⇒ early stress flag (pre-NDVI move).
• Only alert when coverage ≥ 60% valid pixels inside parcel (avoid false alarms on partial scenes).
• Buffer polygons inward 1–2 pixels when summarizing to avoid mixed edges and registration jitter.

Ops notes

• Harmonize across sensors: prefer HLS (S30/L30) for cadence. Keep thresholds relative to each parcel’s baseline.
• Phenology: always compare DOY-aligned windows (±15) to avoid seasonal phase bias.
• Terrain & BRDF: steep terrain → apply topographic correction; exclude high view-zenith angles.
• APIs: serve COGs via TiTiler/WMTS/XYZ; parcel metrics via REST/GraphQL; push alerts via webhook.

References

• NASA Earth Observatory — “Measuring Vegetation (NDVI & EVI)”. Link.
• MODIS Vegetation Indices (MOD13) overview. Link · Product page: MOD13Q1 (V6.1).
• ESA — Sentinel-2 mission & Plant Health. Link.
• Copernicus Climate Change Service — ESOTC 2022 Drought. Link.