Module 5: Drought & Vegetation Monitoring

Satellite-derived vegetation indices and land surface temperature form the backbone of operational drought monitoring. This module covers the physics behind spectral indices, thermal retrieval from Planck inversion, and end-to-end Google Earth Engine pipelines for drought early warning systems combining NDVI anomalies, the Standardized Precipitation Index (SPI), and LST.

1. Spectral Vegetation Indices

Healthy vegetation absorbs strongly in the red (~660 nm) due to chlorophyll pigments and reflects strongly in the near-infrared (~860 nm) due to mesophyll cell structure. This contrast is the foundation for all vegetation indices. By combining reflectances from specific spectral bands, we can quantify vegetation health, water content, and burn severity from space.

The following four indices are the workhorses of operational drought and vegetation monitoring. Each exploits a different aspect of the plant spectral response, and together they provide a comprehensive picture of ecosystem status.

Normalized Difference Vegetation Index (NDVI)

The most widely used vegetation index. NDVI ranges from −1 to +1, where values above 0.3 indicate healthy vegetation and values below 0.1 correspond to bare soil or water.

$$NDVI = \frac{\rho_{NIR} - \rho_{RED}}{\rho_{NIR} + \rho_{RED}}$$

Sentinel-2: B8 (842 nm) and B4 (665 nm). Landsat 8/9: B5 and B4.

Normalized Difference Water Index (NDWI)

Sensitive to leaf water content. Negative NDWI values indicate water stress. Used operationally in drought early warning systems worldwide.

$$NDWI = \frac{\rho_{NIR} - \rho_{SWIR}}{\rho_{NIR} + \rho_{SWIR}}$$

Sentinel-2: B8A (865 nm) and B11 (1610 nm). Landsat 8/9: B5 and B6.

Enhanced Vegetation Index (EVI)

Reduces atmospheric and soil background effects using the blue band. Does not saturate in dense tropical canopies like NDVI, making it preferred for forests.

$$EVI = 2.5 \cdot \frac{\rho_{NIR} - \rho_{RED}}{\rho_{NIR} + 6\rho_{RED} - 7.5\rho_{BLUE} + 1}$$

Sentinel-2: B8, B4, B2. MODIS: standard EVI product (MOD13Q1).

Normalized Burn Ratio (NBR)

Highlights burned areas by exploiting the increase in SWIR reflectance and decrease in NIR reflectance after fire. dNBR (pre-fire minus post-fire) quantifies burn severity.

$$NBR = \frac{\rho_{NIR} - \rho_{SWIR2}}{\rho_{NIR} + \rho_{SWIR2}}$$

Sentinel-2: B8A and B12 (2190 nm). Landsat 8/9: B5 and B7.

Practical Note: Atmospheric Correction

All indices require surface reflectance (bottom-of-atmosphere, BOA) rather than top-of-atmosphere (TOA) values. Sentinel-2 Level-2A products from ESA's Sen2Cor processor provide atmospherically corrected reflectances. For Landsat, use the Collection 2 Level-2 Science Products. Failure to use BOA reflectance introduces systematic biases, especially in the blue and SWIR bands.

2. Land Surface Temperature from Thermal IR

Satellite thermal infrared sensors (Landsat 8/9 TIRS, MODIS bands 31–32, Sentinel-3 SLSTR) measure the radiance emitted by Earth's surface at wavelengths around 10–12 μm. Converting this radiance to physical temperature requires inverting the Planck function and correcting for surface emissivity.

The brightness temperature $T_{bright}$ is obtained by inverting the Planck function assuming a blackbody (emissivity = 1). The true land surface temperature accounts for the actual emissivity $\varepsilon$ of the surface material:

$$T_{LST} = \frac{T_{bright}}{1 + \frac{\lambda \cdot T_{bright}}{\rho_c} \ln\varepsilon}$$

where $\lambda$ is the effective wavelength of the thermal band,$\rho_c = hc/k_B = 1.438 \times 10^{-2}$ m·K is the Planck radiation constant, and $\varepsilon$ is the surface emissivity. Typical emissivity values: vegetation ~0.98, bare soil ~0.93, water ~0.99, urban surfaces ~0.91.

Emissivity Estimation from NDVI

When emissivity measurements are unavailable, a common approach uses the NDVI-based Proportion of Vegetation ($P_v$) method. The fractional vegetation cover is:

$$P_v = \left(\frac{NDVI - NDVI_{soil}}{NDVI_{veg} - NDVI_{soil}}\right)^2$$

Emissivity is then estimated as $\varepsilon = 0.004 \cdot P_v + 0.986$, valid for mixed pixels. Pure soil pixels ($NDVI < 0.2$) use a fixed emissivity of 0.97, and pure vegetation pixels ($NDVI > 0.5$) use 0.99.

Split-Window Technique

For higher accuracy, the split-window method uses two adjacent thermal bands (e.g., MODIS bands 31 and 32 at 11 and 12 μm) to simultaneously solve for atmospheric water vapor absorption and surface temperature. This eliminates the need for independent atmospheric profiles and achieves LST accuracy of ~1 K over most surfaces.

3. Spectral Signatures & Index Computation

The following simulation generates synthetic reflectance spectra for five common land cover types (dense vegetation, sparse vegetation, bare soil, water, and urban) across Sentinel-2's key bands. We then compute NDVI, EVI, and NDWI for each cover type and visualize the results.

Spectral Signatures & Vegetation Index Computation

Python

script.py78 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

# Sentinel-2 band centers (nm) and names
bands = {
    'B2 (Blue)': 490, 'B3 (Green)': 560, 'B4 (Red)': 665,
    'B5 (RE1)': 705, 'B6 (RE2)': 740, 'B7 (RE3)': 783,
    'B8 (NIR)': 842, 'B8A (NIR2)': 865,
    'B11 (SWIR1)': 1610, 'B12 (SWIR2)': 2190
}
band_names = list(bands.keys())
wavelengths = np.array(list(bands.values()))

# Synthetic reflectance profiles (realistic values)
covers = {
    'Dense Vegetation': [0.03, 0.06, 0.04, 0.08, 0.25, 0.42, 0.48, 0.47, 0.18, 0.08],
    'Sparse Vegetation': [0.06, 0.09, 0.08, 0.12, 0.18, 0.25, 0.30, 0.29, 0.22, 0.14],
    'Bare Soil':         [0.10, 0.14, 0.18, 0.20, 0.22, 0.24, 0.26, 0.26, 0.30, 0.28],
    'Water':             [0.06, 0.04, 0.03, 0.02, 0.02, 0.01, 0.01, 0.01, 0.00, 0.00],
    'Urban':             [0.08, 0.10, 0.12, 0.13, 0.14, 0.15, 0.16, 0.16, 0.20, 0.22],
}

# Extract band indices for index computation
iB2, iB4, iB8, iB8A, iB11, iB12 = 0, 2, 6, 7, 8, 9

fig, axes = plt.subplots(1, 2, figsize=(14, 5.5))
colors = ['#22c55e', '#a3e635', '#d97706', '#3b82f6', '#94a3b8']

# Plot 1: Spectral Signatures
ax = axes[0]
for i, (name, refl) in enumerate(covers.items()):
    ax.plot(wavelengths, refl, 'o-', color=colors[i], label=name, linewidth=2, markersize=5)
ax.set_xlabel('Wavelength (nm)', fontsize=12)
ax.set_ylabel('Surface Reflectance', fontsize=12)
ax.set_title('Sentinel-2 Spectral Signatures', fontsize=13, fontweight='bold')
ax.legend(fontsize=9, loc='upper left')
ax.set_xlim(400, 2300)
ax.set_ylim(0, 0.55)
ax.grid(True, alpha=0.3)
ax.axvspan(620, 700, alpha=0.08, color='red', label='Red')
ax.axvspan(770, 900, alpha=0.08, color='green')

# Compute indices for each cover
results = {}
for name, refl in covers.items():
    r = np.array(refl)
    ndvi = (r[iB8] - r[iB4]) / (r[iB8] + r[iB4])
    evi = 2.5 * (r[iB8] - r[iB4]) / (r[iB8] + 6*r[iB4] - 7.5*r[iB2] + 1)
    ndwi = (r[iB8A] - r[iB11]) / (r[iB8A] + r[iB11])
    nbr = (r[iB8A] - r[iB12]) / (r[iB8A] + r[iB12])
    results[name] = {'NDVI': ndvi, 'EVI': evi, 'NDWI': ndwi, 'NBR': nbr}
    print(f"{name:22s}  NDVI={ndvi:+.3f}  EVI={evi:+.3f}  NDWI={ndwi:+.3f}  NBR={nbr:+.3f}")

# Plot 2: Index comparison bar chart
ax2 = axes[1]
x = np.arange(len(covers))
width = 0.2
indices = ['NDVI', 'EVI', 'NDWI', 'NBR']
idx_colors = ['#22c55e', '#06b6d4', '#f59e0b', '#ef4444']
for j, idx_name in enumerate(indices):
    vals = [results[c][idx_name] for c in covers]
    ax2.bar(x + j*width, vals, width, label=idx_name, color=idx_colors[j], alpha=0.85)
ax2.set_xticks(x + 1.5*width)
ax2.set_xticklabels([n.replace(' ', '\n') for n in covers.keys()], fontsize=8)
ax2.set_ylabel('Index Value', fontsize=12)
ax2.set_title('Vegetation Indices by Land Cover', fontsize=13, fontweight='bold')
ax2.legend(fontsize=9)
ax2.grid(True, alpha=0.3, axis='y')
ax2.axhline(y=0, color='white', linewidth=0.5, alpha=0.5)
ax2.set_ylim(-1.1, 1.1)

plt.tight_layout()
plt.savefig('output.png', dpi=130, bbox_inches='tight', facecolor='#0a0a0a')
plt.close()
print("\nSpectral signature and index comparison plots saved.")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

4. Drought Monitoring Framework

Operational drought monitoring integrates multiple satellite-derived indicators into a composite assessment. The key components are:

NDVI Anomaly

The departure of current NDVI from the long-term mean for the same period. Computed as$\Delta NDVI = NDVI_{current} - \overline{NDVI}_{climatology}$. Negative anomalies indicate vegetation stress relative to normal conditions. Typically uses MODIS 16-day composites (MOD13Q1) with a baseline period of 20+ years.

Standardized Precipitation Index (SPI)

A probability-based metric that normalizes precipitation data to a standard Gaussian distribution. SPI values below −1.0 indicate moderate drought, below −1.5 severe drought, and below −2.0 extreme drought. Can be computed at multiple timescales (1, 3, 6, 12 months) to capture different drought characteristics.

LST Anomaly

Elevated land surface temperature amplifies evapotranspiration and accelerates soil moisture depletion. MODIS-derived LST anomalies ($\Delta T_{LST}$) above +2–3 K signal thermal stress. Combined with low NDWI, LST anomalies identify flash drought conditions that develop within weeks.

Vegetation Condition Index (VCI)

VCI normalizes NDVI between historical minimum and maximum values for each pixel:

$$VCI = \frac{NDVI_{current} - NDVI_{min}}{NDVI_{max} - NDVI_{min}} \times 100\%$$

VCI below 35% indicates drought conditions. This approach accounts for differences in vegetation density across regions, making it suitable for continental-scale monitoring.

5. GEE Drought Pipeline (Illustrative)

The following Google Earth Engine script demonstrates a complete drought monitoring pipeline that combines NDVI anomaly detection, Land Surface Temperature, and the Vegetation Condition Index. This is reference code intended for the GEE Code Editor.

// ── GEE Drought Monitoring Pipeline ──
// Illustrative code for Google Earth Engine Code Editor

// 1. Define region and time period
var region = ee.Geometry.Rectangle([32.0, -4.0, 42.0, 5.0]); // East Africa
var year = 2023;
var month = 8;

// 2. MODIS NDVI - current and climatology
var modisNDVI = ee.ImageCollection('MODIS/061/MOD13Q1')
  .filter(ee.Filter.calendarRange(month, month, 'month'));

var ndviCurrent = modisNDVI
  .filter(ee.Filter.calendarRange(year, year, 'year'))
  .select('NDVI').mean().multiply(0.0001);

var ndviClimatology = modisNDVI
  .filter(ee.Filter.calendarRange(2003, 2022, 'year'))
  .select('NDVI').mean().multiply(0.0001);

var ndviAnomaly = ndviCurrent.subtract(ndviClimatology);

// 3. MODIS LST
var modisLST = ee.ImageCollection('MODIS/061/MOD11A2')
  .filter(ee.Filter.calendarRange(month, month, 'month'));

var lstCurrent = modisLST
  .filter(ee.Filter.calendarRange(year, year, 'year'))
  .select('LST_Day_1km').mean().multiply(0.02).subtract(273.15);

var lstClimatology = modisLST
  .filter(ee.Filter.calendarRange(2003, 2022, 'year'))
  .select('LST_Day_1km').mean().multiply(0.02).subtract(273.15);

var lstAnomaly = lstCurrent.subtract(lstClimatology);

// 4. Vegetation Condition Index
var ndviMin = modisNDVI
  .filter(ee.Filter.calendarRange(2003, 2022, 'year'))
  .select('NDVI').reduce(ee.Reducer.min()).multiply(0.0001);

var ndviMax = modisNDVI
  .filter(ee.Filter.calendarRange(2003, 2022, 'year'))
  .select('NDVI').reduce(ee.Reducer.max()).multiply(0.0001);

var vci = ndviCurrent.subtract(ndviMin)
  .divide(ndviMax.subtract(ndviMin))
  .multiply(100);

// 5. Composite drought severity
// Drought = low VCI + negative NDVI anomaly + positive LST anomaly
var droughtSeverity = vci.lt(35)
  .and(ndviAnomaly.lt(-0.05))
  .and(lstAnomaly.gt(2));

// 6. Visualization
Map.centerObject(region, 6);
Map.addLayer(ndviAnomaly.clip(region),
  {min: -0.3, max: 0.3, palette: ['d73027','fc8d59','fee08b','d9ef8b','91cf60','1a9850']},
  'NDVI Anomaly');
Map.addLayer(lstAnomaly.clip(region),
  {min: -5, max: 10, palette: ['313695','4575b4','abd9e9','fee090','fdae61','f46d43','a50026']},
  'LST Anomaly (C)');
Map.addLayer(vci.clip(region),
  {min: 0, max: 100, palette: ['8b0000','ff4500','ffa500','ffff00','7cfc00','228b22']},
  'VCI (%)');
Map.addLayer(droughtSeverity.selfMask().clip(region),
  {palette: ['ff0000']}, 'Severe Drought Areas');

// 7. Print area statistics
var droughtArea = droughtSeverity.selfMask()
  .multiply(ee.Image.pixelArea()).divide(1e6)
  .reduceRegion({reducer: ee.Reducer.sum(), geometry: region, scale: 1000});
print('Drought area (km2):', droughtArea);

Pipeline Explanation

●Step 1–2: Load MODIS NDVI (MOD13Q1, 250 m, 16-day) for the target month. Compute the 20-year climatological mean as baseline.
●Step 3: Load MODIS LST (MOD11A2, 1 km, 8-day). Convert from scaled DN to Celsius and compute anomaly.
●Step 4: Compute VCI by normalizing current NDVI between pixel-wise historical extremes.
●Step 5: Fuse all three indicators: pixels with VCI < 35%, NDVI anomaly < −0.05, and LST anomaly > +2 K are classified as severe drought.

6. NDVI Time Series & Phenology

Long-term NDVI time series reveal vegetation phenology — the seasonal cycle of green-up, peak productivity, senescence, and dormancy. Drought events appear as departures from this normal cycle. The following simulation generates a multi-year NDVI time series with a drought event and performs seasonal decomposition.

NDVI Time Series with Drought Event Detection

Python

script.py101 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

np.random.seed(42)

# Generate 6 years of 16-day NDVI composites (like MODIS)
n_years = 6
composites_per_year = 23  # 365/16
n_obs = n_years * composites_per_year
t = np.arange(n_obs)
doy = np.tile(np.linspace(1, 365, composites_per_year), n_years)
years = np.repeat(np.arange(2018, 2018 + n_years), composites_per_year)

# Base phenology: double logistic model (green-up + senescence)
def phenology(doy, ndvi_min=0.15, ndvi_max=0.72, sos=90, eos=280, rate=0.05):
    greenup = ndvi_min + (ndvi_max - ndvi_min) / (1 + np.exp(-rate * (doy - sos)))
    senescence = ndvi_min + (ndvi_max - ndvi_min) / (1 + np.exp(rate * (doy - eos)))
    return np.minimum(greenup, senescence)

ndvi_base = phenology(doy)

# Add drought in year 2021 (index 3): reduce peak NDVI by 35%
drought_mask = (years == 2021) & (doy > 120) & (doy < 300)
ndvi = ndvi_base.copy()
ndvi[drought_mask] *= 0.65

# Add realistic noise
ndvi += np.random.normal(0, 0.02, n_obs)
ndvi = np.clip(ndvi, 0, 1)

# Compute climatological mean and std (excluding drought year)
clim_mean = np.zeros(composites_per_year)
clim_std = np.zeros(composites_per_year)
for i in range(composites_per_year):
    mask = np.arange(i, n_obs, composites_per_year)
    non_drought = [m for m in mask if years[m] != 2021]
    clim_mean[i] = np.mean(ndvi[non_drought])
    clim_std[i] = np.std(ndvi[non_drought])

# Anomaly for 2021
idx_2021 = years == 2021
anomaly_2021 = ndvi[idx_2021] - clim_mean

fig, axes = plt.subplots(3, 1, figsize=(13, 10))

# Panel 1: Full time series
ax = axes[0]
dates = t
for yr in range(2018, 2024):
    mask = years == yr
    color = '#ef4444' if yr == 2021 else '#38bdf8'
    lw = 2.5 if yr == 2021 else 1.0
    alpha = 1.0 if yr == 2021 else 0.5
    ax.plot(t[mask], ndvi[mask], color=color, linewidth=lw, alpha=alpha, label=str(yr))
ax.set_ylabel('NDVI', fontsize=12)
ax.set_title('Multi-Year NDVI Time Series (Drought in 2021)', fontsize=13, fontweight='bold', color='white')
ax.legend(ncol=6, fontsize=9, loc='upper right')
ax.grid(True, alpha=0.3)
ax.set_facecolor('#0a0a0a')

# Panel 2: 2021 vs climatology
ax2 = axes[1]
x_doy = np.linspace(1, 365, composites_per_year)
ax2.fill_between(x_doy, clim_mean - 2*clim_std, clim_mean + 2*clim_std, alpha=0.2, color='#38bdf8', label='Mean +/- 2 sigma')
ax2.plot(x_doy, clim_mean, '--', color='#38bdf8', linewidth=2, label='Climatological mean')
ax2.plot(x_doy, ndvi[idx_2021], color='#ef4444', linewidth=2.5, label='2021 (drought)')
ax2.set_ylabel('NDVI', fontsize=12)
ax2.set_xlabel('Day of Year', fontsize=12)
ax2.set_title('2021 vs Climatology', fontsize=13, fontweight='bold', color='white')
ax2.legend(fontsize=9)
ax2.grid(True, alpha=0.3)
ax2.set_facecolor('#0a0a0a')

# Panel 3: VCI
vci = (ndvi[idx_2021] - np.tile(clim_mean - 2*clim_std, 1)[:composites_per_year]) /       (np.tile(clim_mean + 2*clim_std, 1)[:composites_per_year] - np.tile(clim_mean - 2*clim_std, 1)[:composites_per_year]) * 100
ax3 = axes[2]
colors_vci = ['#ef4444' if v < 35 else '#f59e0b' if v < 50 else '#22c55e' for v in vci]
ax3.bar(x_doy, vci, width=14, color=colors_vci, alpha=0.8)
ax3.axhline(y=35, color='#ef4444', linestyle='--', linewidth=1.5, label='Drought threshold (35%)')
ax3.set_ylabel('VCI (%)', fontsize=12)
ax3.set_xlabel('Day of Year', fontsize=12)
ax3.set_title('Vegetation Condition Index - 2021', fontsize=13, fontweight='bold', color='white')
ax3.legend(fontsize=9)
ax3.grid(True, alpha=0.3)
ax3.set_facecolor('#0a0a0a')
ax3.set_ylim(0, 110)

for ax in axes:
    ax.tick_params(colors='white')
    for spine in ax.spines.values():
        spine.set_color('#334155')

plt.tight_layout()
plt.savefig('output.png', dpi=130, bbox_inches='tight', facecolor='#0a0a0a')
plt.close()
print("Drought year 2021: peak NDVI reduced by ~35%")
print(f"Min VCI in 2021: {vci.min():.1f}%")
print(f"Composites below drought threshold (VCI<35%): {np.sum(vci < 35)} of {len(vci)}")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

7. Key Takeaways

✓NDVI, EVI, NDWI, and NBR each exploit different spectral contrasts to monitor vegetation health, water stress, and burn severity.

✓Land Surface Temperature from thermal IR requires Planck inversion and emissivity correction; the split-window technique improves accuracy to ~1 K.

✓Multi-indicator drought monitoring (NDVI anomaly + VCI + LST anomaly) outperforms any single index by capturing both meteorological and agricultural drought.

✓MODIS provides the longest consistent record (20+ years) at 250 m; Sentinel-2 offers 10 m resolution for field-scale drought assessment.

✓Google Earth Engine enables planetary-scale drought analysis without downloading any data locally.