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.

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.

In [1]:

Copied!





from pathlib import Path
from earthlens import EarthLens, AggregationConfig

OUT = Path("data/era5-madrid-spi")
OUT.mkdir(parents=True, exist_ok=True)

earthlens = EarthLens(
    data_source="ecmwf",
    temporal_resolution="monthly",
    start="1999-01-01",
    end="2023-12-01",
    variables={
        "reanalysis-era5-single-levels-monthly-means": [
            "total-precipitation",
        ],
    },
    lat_lim=[40.0, 41.0],
    lon_lim=[-4.0, -3.0],
    path=str(OUT),
)
earthlens.download(aggregate=AggregationConfig(freq="1MS", op="auto"))
from pathlib import Path
from earthlens import EarthLens, AggregationConfig

OUT = Path("data/era5-madrid-spi")
OUT.mkdir(parents=True, exist_ok=True)

earthlens = EarthLens(
    data_source="ecmwf",
    temporal_resolution="monthly",
    start="1999-01-01",
    end="2023-12-01",
    variables={
        "reanalysis-era5-single-levels-monthly-means": [
            "total-precipitation",
        ],
    },
    lat_lim=[40.0, 41.0],
    lon_lim=[-4.0, -3.0],
    path=str(OUT),
)
earthlens.download(aggregate=AggregationConfig(freq="1MS", op="auto"))

2026-05-10 01:48:39.347 | INFO     | earthlens.ecmwf.backend:download:536 - Download ECMWF reanalysis-era5-single-levels-monthly-means/total-precipitation data for period 1999-01-01 00:00:00 till 2023-12-01 00:00:00

2026-05-10 01:48:39.980 | INFO     | earthlens.ecmwf.backend:_api:724 - Requesting reanalysis-era5-single-levels-monthly-means from CDS; this may take several minutes

2026-05-10 01:48:41,142 INFO Request ID is 03cd6935-e248-4d41-b9e6-626fffd05e7a

2026-05-10 01:48:41,300 INFO status has been updated to accepted

2026-05-10 01:49:03,014 INFO status has been updated to running

2026-05-10 01:49:14,477 INFO status has been updated to successful

2026-05-10 01:49:15 | INFO | pyramids.base.config | Logging is configured.

2026-05-10 01:49:18.970 | INFO     | earthlens.ecmwf.backend:download:575 - ECMWF download summary: all 1 variables succeeded ([('reanalysis-era5-single-levels-monthly-means', 'total-precipitation')])

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.

In [2]:

Copied!





import numpy as np
import pandas as pd
from pyramids.dataset import Dataset

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")

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()

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
import numpy as np
import pandas as pd
from pyramids.dataset import Dataset

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")

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()

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

Out[2]:

	P [mm/mo]	P-12 [mm]	SPI-12
2023-01-01	1.2	19.0	0.67
2023-02-01	0.3	19.0	0.57
2023-03-01	0.4	15.0	-0.72
2023-04-01	0.3	13.0	-1.50
2023-05-01	1.9	15.0	-0.87
2023-06-01	2.5	17.0	-0.22
2023-07-01	0.1	17.0	-0.21
2023-08-01	0.0	16.0	-0.36
2023-09-01	3.5	18.0	0.25
2023-10-01	4.1	21.0	0.92
2023-11-01	1.5	20.0	0.73
2023-12-01	1.2	17.0	-0.12

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.

In [3]:

Copied!





import matplotlib.pyplot as plt

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()
import matplotlib.pyplot as plt

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()

No description has been provided for this image

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-radiation and 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.