Drought monitoring — SPI from monthly precipitation¶
The Standardized Precipitation Index (SPI) is the WMO-recommended meteorological-drought metric. It expresses the precipitation total over a window (e.g. 3 or 12 months) as a z-score against the long-term distribution of that same window. Negative SPI → dry, positive → wet, |SPI| > 1.5 → severe.
Domain context. SPI's appeal is that it normalizes across climates — a -2 SPI value means roughly the same drought severity in Spain as in the Sahel even though the absolute precipitation totals are wildly different.
This notebook computes SPI-12 (rolling 12-month sums) for a Mediterranean region using ERA5 monthly precipitation.
Setup¶
Consolidated imports for the whole notebook. pyramids provides
Dataset (raster reading); earthlens provides the unified EarthLens
entry point plus AggregationConfig for the monthly-sum aggregation.
from pathlib import Path
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from pyramids.dataset import Dataset
from earthlens import AggregationConfig, EarthLens
Step 1 — pull 25 years of monthly precipitation¶
Box: 1° around Madrid, Spain (40°–41°N, 4°–3°W). 25 years (1999–2023)
is barely enough for a robust SPI fit; operational drought monitors
use 50+. Precipitation is a flux → op="auto" → sum.
Output directory¶
ECMWF is a file-writing backend, so we first create the folder the downloaded NetCDF and its aggregated GeoTIFFs will land in.
OUT = Path("data/era5-madrid-spi")
OUT.mkdir(parents=True, exist_ok=True)
Build the request¶
Describe the ERA5 pull — source, monthly cadence, 25-year window, the
single-levels monthly-means dataset, total-precipitation, and the 1° box
around Madrid. Nothing is fetched yet.
earthlens = EarthLens(
data_source="ecmwf",
cadence="monthly",
start="1999-01-01",
end="2023-12-01",
dataset="reanalysis-era5-single-levels-monthly-means",
variables=[
"total-precipitation",
],
aoi=[-4.0, 40.0, -3.0, 41.0],
path=str(OUT),
)
Download and aggregate¶
Run the download as its own step. Precipitation is a flux, so op="auto"
resolves to a monthly sum; the 1MS frequency keeps one value per month.
earthlens.download(aggregate=AggregationConfig(freq="1MS", op="auto"))
Step 2 — compute SPI-12¶
- Spatial-mean each month → one number per month.
- Rolling-12-month sum → 12-month accumulated precipitation.
- Fit each calendar-month's distribution of accumulations to a normal distribution and compute the z-score (full SPI uses a gamma distribution and CDF transform; we'll use the simpler normal z-score here for clarity).
- Plot the monthly SPI-12 series.
Read the monthly precipitation series¶
Read each aggregated GeoTIFF, take its spatial mean, and convert from metres to millimetres — one precipitation total per month.
agg = OUT / "aggregated"
paths = sorted(agg.glob("total_precipitation_1MS_*.tif"))
monthly_mm = np.array(
[
np.nanmean(Dataset.read_file(str(p)).read_array()) * 1000 # m -> mm
for p in paths
]
)
months = pd.date_range("1999-01-01", periods=len(monthly_mm), freq="MS")
Rolling-12 accumulation and SPI z-score¶
Sum each trailing 12-month window, then standardize per calendar month: subtract that month's mean accumulation and divide by its standard deviation to get the SPI-12 z-score.
p12 = pd.Series(monthly_mm, index=months).rolling(12).sum()
# Per-calendar-month z-score on the rolling-12 series.
spi = p12.copy() * np.nan
for m in range(1, 13):
mask = p12.index.month == m
vals = p12[mask]
spi[mask] = (vals - vals.mean()) / vals.std()
Inspect the most recent year¶
Tabulate the last 12 months — monthly precipitation, the 12-month accumulation, and the resulting SPI-12 — as a quick sanity check.
tail = pd.DataFrame(
{
"P [mm/mo]": monthly_mm[-12:].round(1),
"P-12 [mm]": p12.tail(12).round(0).to_numpy(),
"SPI-12": spi.tail(12).round(2).to_numpy(),
},
index=months[-12:],
)
tail
Step 3 — plot SPI-12 with severity bands¶
Standard SPI categories: |SPI| < 1 normal; 1.0–1.5 mild; 1.5–2.0 moderate; > 2.0 severe / extreme.
Plot the SPI-12 series¶
Bar each month (red for dry, blue for wet) and shade the mild / moderate / severe dry bands so drought episodes stand out at a glance.
fig, ax = plt.subplots(figsize=(11, 5))
ax.bar(
spi.index,
spi.values,
color=["tab:red" if v < 0 else "tab:blue" for v in spi.values],
width=20,
alpha=0.7,
)
ax.axhspan(-1.5, -1.0, color="orange", alpha=0.10, label="Mild dry")
ax.axhspan(-2.0, -1.5, color="red", alpha=0.10, label="Moderate dry")
ax.axhspan(-3.0, -2.0, color="darkred", alpha=0.15, label="Severe dry")
ax.axhline(0, color="black", lw=0.5)
ax.set_ylabel("SPI-12")
ax.set_ylim(-3, 3)
ax.set_title("SPI-12 over Madrid bbox, 1999–2023 (ERA5 precipitation)")
ax.legend(loc="upper right")
ax.grid(alpha=0.3)
plt.tight_layout()
plt.show()
Notes¶
- Normal vs gamma SPI. Operational SPI fits a gamma distribution (better for skewed precipitation) then transforms to a standard normal via the CDF — see McKee et al. (1993) and the WMO SPI user guide. The normal z-score used here is a fine first-pass.
- Window length. SPI-1 / SPI-3 / SPI-12 measure short / medium / long droughts respectively. Agriculture often watches SPI-3; reservoirs care about SPI-12 / SPI-24.
- Combine with temperature for SPEI. The Standardized
Precipitation-Evapotranspiration Index uses $P - PET$ rather than
$P$ alone, capturing the warming-driven increase in atmospheric
demand. Compute PET from ERA5
2m-temperature/surface-net-solar-radiationand follow the same z-score recipe. - Spatial drought maps. Apply this same recipe pixel-by-pixel rather than after spatial-mean. The result is one map of SPI-12 per month — the natural input to e.g. crop-failure attribution studies.