Building a Drought Index Pipeline with NDMI

Drought is an anomaly, so the pipeline compares each new NDMI value against what is normal for that pixel at that time of year:

ndmi = (nir - swir16) / (nir + swir16)
period = ndmi.time.dt.dayofyear // 16                        # 16-day periods
base = hist_ndmi.groupby(hist_ndmi.time.dt.dayofyear // 16)
z = (ndmi.groupby(period) - base.mean()).groupby(period) / base.std()

A single NDMI value of 0.1 is healthy for a dry grassland in August and alarming for an irrigated field in June; only the anomaly carries the drought signal. This page belongs to spectral index calculation pipelines in Satellite Processing Workflows & Index Pipelines.

From Raw Index to Drought Class

Four stages of a moisture drought pipeline Clear observations produce NDMI values. Several years of history produce a per-pixel, per-period baseline of mean and standard deviation. Each new value is converted to a standardised anomaly against its period. Thresholds on the anomaly produce drought severity classes. Index → baseline → anomaly → class NDMI NIR vs SWIR-1 baseline mean, std per period anomaly z-score severity none → extreme The baseline is built once and refreshed yearly; the anomaly runs on every new scene.

NDMI — the normalised difference of near-infrared and shortwave-infrared around 1.6 µm — tracks the water content of leaves, because liquid water absorbs strongly in the shortwave infrared. It responds to moisture stress earlier than NDVI, which reacts only once leaves discolour or drop. That lead time is the reason to build a moisture index into a drought pipeline rather than relying on greenness alone.

Environment & Setup

Package Version pin Used for
xarray >=2023.1 Grouping by period, anomalies
odc-stac >=0.3.9 Loading multi-year stacks
dask >=2023.1 Lazy baseline computation
rioxarray >=0.15 Writing baseline and class rasters 
pip install "xarray>=2023.1" "odc-stac>=0.3.9" "dask>=2023.1" "rioxarray>=0.15"

Complete Working Example

import numpy as np
import xarray as xr

PERIOD_DAYS = 16
CLASSES = [(-np.inf, -2.0, 4), (-2.0, -1.5, 3), (-1.5, -1.0, 2), (-1.0, -0.5, 1), (-0.5, np.inf, 0)]


def ndmi(ds: xr.Dataset) -> xr.DataArray:
    nir, swir = ds.nir.astype("float32") * 1e-4, ds.swir16.astype("float32") * 1e-4
    clear = ~ds.scl.isin([0, 1, 3, 8, 9, 10, 11])
    return ((nir - swir) / (nir + swir)).where(clear).rename("ndmi")


def period_of(t: xr.DataArray) -> xr.DataArray:
    return (t.dt.dayofyear - 1) // PERIOD_DAYS


def build_baseline(hist: xr.DataArray, min_obs: int = 5) -> xr.Dataset:
    # composite each period of each year first, so cloudy years do not dominate
    per = hist.resample(time=f"{PERIOD_DAYS}D").median()
    g = per.groupby(period_of(per.time))
    base = xr.Dataset({"mean": g.mean(), "std": g.std(), "n": g.count()})
    return base.where(base.n >= min_obs)


def anomaly(obs: xr.DataArray, base: xr.Dataset) -> xr.DataArray:
    p = period_of(obs.time)
    b = base.sel(time=p)                     # 'time' is the groupby coordinate name
    z = (obs - b["mean"]) / b["std"].clip(min=0.02)
    return z.rename("ndmi_z")


def classify(z: xr.DataArray) -> xr.DataArray:
    out = xr.full_like(z, 255, dtype="uint8")
    for lo, hi, code in CLASSES:
        out = out.where(~((z > lo) & (z <= hi)), code)
    return out.where(z.notnull(), 255).rename("drought_class")

Two details matter for a trustworthy baseline. First, composite each period within each year before computing statistics; otherwise a year with twelve clear scenes in a period outweighs a year with one, and the baseline reflects cloud climatology as much as moisture. Second, clip the standard deviation at a small floor — in irrigated or evergreen areas it can be tiny, and dividing by it turns noise into apparent extreme drought.

Choosing the Baseline Length

More years steady the baseline — up to a point The uncertainty of the per-pixel baseline standard deviation falls quickly between two and five years of history and slowly after that. Beyond about ten years, land-use change and sensor differences begin to matter more than sampling error. Five to eight years of Sentinel-2 is a practical baseline. Baseline uncertainty vs years of history 5–8 years 2 yr 15 yr Longer records need harmonised Landsat, which adds cross-sensor error.

Sentinel-2 offers since 2017 a record long enough for a five- to eight-year baseline, which is enough to estimate a mean reliably and a standard deviation reasonably. Extending further back needs Landsat, which in turn needs harmonisation so that a sensor change is not mistaken for a drought; see harmonising Landsat and Sentinel-2 reflectance. A baseline that includes a major drought year is also biased towards dryness, so some programmes exclude known extreme years or use the median and a robust spread in place of mean and standard deviation.

Interpreting Anomalies Carefully

A negative NDMI anomaly means the canopy holds less water than usual for the date — but not always because of drought. Harvest, mowing, fire and clearance all produce the same signal, abruptly, and a phenology shift of two weeks produces a transient anomaly at green-up or senescence with no stress at all. Two practical filters help. Require an anomaly to persist over two or more consecutive periods before reporting it, which rejects most management events and phenology shifts. And mask land-cover classes where the index is not meaningful — built-up areas, water, bare rock — using a land-cover map, so the drought map only speaks for vegetation.

Publishing the Output

A drought product is usually consumed as a map and a table. Write the class raster as a COG with a colour table so it renders correctly without styling, and the z-score as float32 for analysts. Aggregate the classes to administrative or catchment units using the methods in zonal statistics and vector–raster integration — the fraction of each unit’s vegetated area in each class is the number decision-makers read. Keep every period’s outputs rather than overwriting, so trends through the season can be shown.

Verification

A normal year should look normal Applying the pipeline to a year inside the baseline period should give z-scores centred on zero with a standard deviation near one. A shift indicates a bias in scaling or masking, and a spread much larger than one indicates the standard deviation floor is too low. z-score histogram for a baseline year 0 Mean near 0 and spread near 1 means the baseline and scaling agree. 
import numpy as np

z = anomaly(ndmi(ds_2021), base)
print("mean", float(z.mean()), "std", float(z.std()))
assert abs(float(z.mean())) < 0.2 and 0.7 < float(z.std()) < 1.4
cls = classify(z)
print(np.unique(cls.values, return_counts=True))

Common Errors

Irrigated fields always show extreme drought

The baseline standard deviation is tiny there. Raise the floor on the standard deviation.

The whole scene shifts after early 2022

Sentinel-2’s processing baseline added an offset. Apply the offset before computing NDMI.

Harvested fields appear as drought

A single-period anomaly after harvest. Require persistence over two periods and mask cropland harvest windows.

Baseline has holes in winter

Too few clear observations per period. Lengthen the period or lower min_obs and flag low-confidence pixels.

Frequently Asked Questions

Q: Why use NDMI rather than NDVI for drought? NDMI responds to leaf water content, which drops before leaves lose chlorophyll, so it signals moisture stress earlier. NDVI remains useful for confirming impact once canopy greenness declines.

Q: How long a baseline do I need? Five to eight years is practical with Sentinel-2. Shorter baselines give unstable standard deviations; longer ones need harmonised Landsat.

Q: Which SWIR band should NDMI use? The band near 1.6 µm — Sentinel-2 B11 or Landsat 8/9 band 6. Using the 2.2 µm band gives a related but different index.

Q: Why standardise instead of using raw NDMI? Because normal NDMI varies enormously between land covers and through the season. Standardising against each pixel’s own history makes a single threshold meaningful everywhere.

Q: Can this run operationally? Yes. The baseline is computed once a year; each new scene needs only the index, a lookup of the period’s statistics and a classification, which is cheap enough to run on every acquisition.