Showcase — the Horn of Africa drought (CHIRPS rainfall)¶

Between 2020 and 2023 the Horn of Africa endured its worst drought in over 40 years: five consecutive rainy seasons failed, and the 2022 season was the driest on record, contributing to an estimated 43,000 deaths in Somalia. CHIRPS (Climate Hazards Center) is the standard quasi-global rainfall dataset for drought monitoring; the earthlens CHC backend fetches it over anonymous FTP.

What this notebook does¶

1. Download CHIRPS monthly rainfall for the March–May long rains of every year 2010–2023.
2. Build a seasonal-rainfall time series for the region and compare each year against the long-term average (a simple standardized anomaly, the basis of the SPI).
3. Map a normal season (2018) against the failed 2022 season and the deficit.
4. Cross-check on the ground with MODIS NDVI greenness over the same window, via Google Earth Engine.

Live CHIRPS FTP download — no credentials required (this step takes a few minutes). The NDVI section needs GEE_SERVICE_ACCOUNT / GEE_SERVICE_KEY.

Setup¶

First the imports and a bit of noise-suppression. pyramids provides the in-memory raster datasets that load() returns (we call read_array() on them); earthlens provides the unified EarthLens entry point.

In [1]:

import os
import tempfile
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from dotenv import find_dotenv, load_dotenv
from loguru import logger

from earthlens import EarthLens

logger.disable("earthlens")
logger.disable("pyramids")

import os import tempfile from pathlib import Path import matplotlib.pyplot as plt import numpy as np import pandas as pd from dotenv import find_dotenv, load_dotenv from loguru import logger from earthlens import EarthLens logger.disable("earthlens") logger.disable("pyramids")

Credentials¶

The CHIRPS section needs no credentials (anonymous FTP). The optional NDVI section reads its Google Earth Engine credentials from a .env file at the repository root (loaded with python-dotenv; your PATH is left untouched). Create that .env with:

GEE_SERVICE_ACCOUNT="your-sa@your-project.iam.gserviceaccount.com"
GEE_SERVICE_KEY="/path/to/your/service-account-key.json"

In [2]:

load_dotenv(find_dotenv(usecwd=True))

load_dotenv(find_dotenv(usecwd=True))

Out[2]:

True

1 · Download CHIRPS long-rains rainfall¶

CHIRPS is a raster dataset, so the CHC backend's load() returns in-memory pyramids datasets — one per monthly file. We pull the March–May long rains for each year and sum the three months into a single seasonal total over the Horn of Africa.

Region of interest¶

The bounding box covers Somalia, Ethiopia and Kenya — the core of the drought. We keep LAT / LON as module-level constants so the NDVI section can reuse the same extent.

In [3]:

LAT, LON = [-5.0, 12.0], [38.0, 52.0]  # Horn of Africa (Somalia / Ethiopia / Kenya)
# One bounded output root for the whole run; pass it to load() so each
# year writes to a known directory instead of a throwaway temp dir.
OUT_DIR = Path(tempfile.mkdtemp(prefix="horn-of-africa-"))

LAT, LON = [-5.0, 12.0], [38.0, 52.0] # Horn of Africa (Somalia / Ethiopia / Kenya) # One bounded output root for the whole run; pass it to load() so each # year writes to a known directory instead of a throwaway temp dir. OUT_DIR = Path(tempfile.mkdtemp(prefix="horn-of-africa-"))

A helper for one season¶

mam_total() builds a CHC request for one year's March–May window, loads the monthly rasters into memory, reads each into an array, masks the no-data sentinels, and sums the three months. The construct → load steps are kept on separate lines for readability.

In [4]:

def mam_total(year):
    """Sum the Mar-May CHIRPS monthly rasters for one year into a season total (mm)."""
    chirps = EarthLens(
        data_source="chc",
        cadence="monthly",
        variables=["precipitation"],
        start=f"{year}-03-01",
        end=f"{year}-05-31",
        aoi=[LON[0], LAT[0], LON[1], LAT[1]],
        path=str(OUT_DIR / "chirps" / str(year)),
    )
    rasters = chirps.load(progress_bar=False)
    grids = [np.asarray(r.read_array(), dtype="float32") for r in rasters]
    grids = [g[0] if g.ndim == 3 else g for g in grids]
    arr = np.where(np.array(grids) < 0, np.nan, grids)
    return np.nansum(arr, axis=0)  # mm summed over Mar+Apr+May

def mam_total(year): """Sum the Mar-May CHIRPS monthly rasters for one year into a season total (mm).""" chirps = EarthLens( data_source="chc", cadence="monthly", variables=["precipitation"], start=f"{year}-03-01", end=f"{year}-05-31", aoi=[LON[0], LAT[0], LON[1], LAT[1]], path=str(OUT_DIR / "chirps" / str(year)), ) rasters = chirps.load(progress_bar=False) grids = [np.asarray(r.read_array(), dtype="float32") for r in rasters] grids = [g[0] if g.ndim == 3 else g for g in grids] arr = np.where(np.array(grids) < 0, np.nan, grids) return np.nansum(arr, axis=0) # mm summed over Mar+Apr+May

Download every season (2010–2023)¶

Loop the helper over all 14 years. This is the live FTP download and takes a few minutes; the result is one season-total array per year.

In [5]:

YEARS = list(range(2010, 2024))
seasons = {}
for y in YEARS:
    seasons[y] = mam_total(y)
print(f"downloaded {len(seasons)} long-rains seasons (2010-2023)")

YEARS = list(range(2010, 2024)) seasons = {} for y in YEARS: seasons[y] = mam_total(y) print(f"downloaded {len(seasons)} long-rains seasons (2010-2023)")

2026-06-13 18:05:45 | INFO | pyramids.base.config | Logging is configured.

downloaded 14 long-rains seasons (2010-2023)

2 · Seasonal rainfall time series and standardized anomaly¶

Averaging each season over the region gives one rainfall total per year. Expressed as a standardized anomaly (year − mean) / std — the idea behind the SPI — the failed seasons stand out as strongly negative (drought) bars.

Regional means and the SPI-like z-score¶

Reduce each season to its regional mean, then standardize the 14-year series into a z-score.

In [6]:

yrs = sorted(seasons)
totals = np.array([np.nanmean(seasons[y]) for y in yrs])
z = (totals - totals.mean()) / totals.std()

yrs = sorted(seasons) totals = np.array([np.nanmean(seasons[y]) for y in yrs]) z = (totals - totals.mean()) / totals.std()

Plot rainfall and anomaly side by side¶

Left: the raw seasonal rainfall with the 2010–23 mean. Right: the standardized anomaly, red for below-average (drought) years.

In [7]:

fig, ax = plt.subplots(1, 2, figsize=(14, 4))

ax[0].bar(yrs, totals, color="tab:blue")
ax[0].axhline(totals.mean(), color="k", ls="--", label="2010-23 mean")
ax[0].set(ylabel="mean MAM rainfall (mm)", title="Long-rains rainfall by year")
ax[0].legend()

colors = ["firebrick" if v < 0 else "tab:green" for v in z]
ax[1].bar(yrs, z, color=colors)
ax[1].axhline(0, color="k", lw=0.8)
ax[1].set(
    ylabel="standardized anomaly (σ)",
    title="Standardized rainfall anomaly (SPI-like)",
)

plt.tight_layout()
plt.show()

fig, ax = plt.subplots(1, 2, figsize=(14, 4)) ax[0].bar(yrs, totals, color="tab:blue") ax[0].axhline(totals.mean(), color="k", ls="--", label="2010-23 mean") ax[0].set(ylabel="mean MAM rainfall (mm)", title="Long-rains rainfall by year") ax[0].legend() colors = ["firebrick" if v < 0 else "tab:green" for v in z] ax[1].bar(yrs, z, color=colors) ax[1].axhline(0, color="k", lw=0.8) ax[1].set( ylabel="standardized anomaly (σ)", title="Standardized rainfall anomaly (SPI-like)", ) plt.tight_layout() plt.show()

The 2021 and 2022 seasons were among the most severe; print their totals and anomalies.

In [8]:

for y in (2021, 2022):
    if y in seasons:
        print(f"{y}: {totals[yrs.index(y)]:.0f} mm  ({z[yrs.index(y)]:+.1f}σ)")

for y in (2021, 2022): if y in seasons: print(f"{y}: {totals[yrs.index(y)]:.0f} mm ({z[yrs.index(y)]:+.1f}σ)")

2021: 92 mm  (-0.6σ)
2022: 63 mm  (-1.6σ)

3 · Maps: a normal season vs the failed 2022 season¶

The time series collapses each season to a single number; the maps show where the deficit fell. We compare 2018 (near-normal) against 2022 (failed) and the per-pixel difference.

Align the two seasons¶

Crop both arrays to their common shape and pick a shared colour scale from the 98th percentile so the two maps are directly comparable.

In [9]:

wet, dry = seasons[2018], seasons[2022]
h, w = min(wet.shape[0], dry.shape[0]), min(wet.shape[1], dry.shape[1])
wet, dry = wet[:h, :w], dry[:h, :w]
vm = float(np.nanpercentile(np.concatenate([wet.ravel(), dry.ravel()]), 98))

wet, dry = seasons[2018], seasons[2022] h, w = min(wet.shape[0], dry.shape[0]), min(wet.shape[1], dry.shape[1]) wet, dry = wet[:h, :w], dry[:h, :w] vm = float(np.nanpercentile(np.concatenate([wet.ravel(), dry.ravel()]), 98))

Render the three maps¶

2018 rainfall, 2022 rainfall, and the 2022 − 2018 deficit.

In [10]:

fig, ax = plt.subplots(1, 3, figsize=(16, 5))

ax[0].imshow(wet, cmap="YlGnBu", vmin=0, vmax=vm)
ax[0].set_title("MAM rainfall — 2018 (mm)")
ax[1].imshow(dry, cmap="YlGnBu", vmin=0, vmax=vm)
ax[1].set_title("MAM rainfall — 2022 (mm)")
im = ax[2].imshow(dry - wet, cmap="BrBG", vmin=-vm, vmax=vm)
ax[2].set_title("deficit (2022 − 2018)")

for a in ax:
    a.axis("off")
fig.colorbar(im, ax=ax[2], fraction=0.046, label="mm")
plt.tight_layout()
plt.show()

fig, ax = plt.subplots(1, 3, figsize=(16, 5)) ax[0].imshow(wet, cmap="YlGnBu", vmin=0, vmax=vm) ax[0].set_title("MAM rainfall — 2018 (mm)") ax[1].imshow(dry, cmap="YlGnBu", vmin=0, vmax=vm) ax[1].set_title("MAM rainfall — 2022 (mm)") im = ax[2].imshow(dry - wet, cmap="BrBG", vmin=-vm, vmax=vm) ax[2].set_title("deficit (2022 − 2018)") for a in ax: a.axis("off") fig.colorbar(im, ax=ax[2], fraction=0.046, label="mm") plt.tight_layout() plt.show()

Quantify the regional drop between the two seasons.

In [11]:

print(
    f"mean MAM rainfall  2018: {np.nanmean(wet):.0f} mm   2022: {np.nanmean(dry):.0f} mm   "
    f"({(np.nanmean(dry)/np.nanmean(wet)-1)*100:+.0f}%)"
)

print( f"mean MAM rainfall 2018: {np.nanmean(wet):.0f} mm 2022: {np.nanmean(dry):.0f} mm " f"({(np.nanmean(dry)/np.nanmean(wet)-1)*100:+.0f}%)" )

mean MAM rainfall  2018: 159 mm   2022: 63 mm   (-60%)

4 · Vegetation response (MODIS NDVI)¶

Failed rains show up on the ground as browning vegetation. Pulling MODIS NDVI (the standard greenness index) for the same March–May window over the Horn — this time from a different provider, Google Earth Engine — shows how the 2022 rainfall deficit translated into a collapse in greenness, the basis of the Vegetation Condition Index used in famine early-warning.

GEE credentials¶

Read the Earth Engine service-account credentials from the environment (loaded from .env above). If they are absent the NDVI section is skipped.

In [12]:

SERVICE_ACCOUNT = os.environ.get("GEE_SERVICE_ACCOUNT", "")
SERVICE_KEY = os.environ.get("GEE_SERVICE_KEY", "")

SERVICE_ACCOUNT = os.environ.get("GEE_SERVICE_ACCOUNT", "") SERVICE_KEY = os.environ.get("GEE_SERVICE_KEY", "")

A helper for one NDVI window¶

mam_ndvi() builds a GEE request for the same LAT / LON extent and March–May window, authenticates, loads the single composite raster, and returns the scaled NDVI array. The construct → authenticate → load steps are kept on separate lines.

In [13]:

def mam_ndvi(year):
    """Mean MODIS NDVI over the Horn for one year's Mar-May window (scaled to -1..1)."""
    scene = EarthLens(
        data_source="gee",
        start=f"{year}-03-01",
        end=f"{year}-05-31",
        cadence="raw",
        dataset="MODIS/061/MOD13A2",
        variables=["NDVI"],
        aoi=[LON[0], LAT[0], LON[1], LAT[1]],
        scale=5000,
        reducer="mean",
        path=str(OUT_DIR / "modis" / str(year)),
    )
    scene.authenticate(service_account=SERVICE_ACCOUNT, service_key=SERVICE_KEY)
    rasters = scene.load(progress_bar=False)
    if not rasters:
        return None
    arr = np.asarray(rasters[0].read_array(), dtype="float32")
    return (arr[0] if arr.ndim == 3 else arr) * 0.0001  # MODIS NDVI scale factor

def mam_ndvi(year): """Mean MODIS NDVI over the Horn for one year's Mar-May window (scaled to -1..1).""" scene = EarthLens( data_source="gee", start=f"{year}-03-01", end=f"{year}-05-31", cadence="raw", dataset="MODIS/061/MOD13A2", variables=["NDVI"], aoi=[LON[0], LAT[0], LON[1], LAT[1]], scale=5000, reducer="mean", path=str(OUT_DIR / "modis" / str(year)), ) scene.authenticate(service_account=SERVICE_ACCOUNT, service_key=SERVICE_KEY) rasters = scene.load(progress_bar=False) if not rasters: return None arr = np.asarray(rasters[0].read_array(), dtype="float32") return (arr[0] if arr.ndim == 3 else arr) * 0.0001 # MODIS NDVI scale factor

Load NDVI for both seasons¶

Run the helper for 2018 and 2022. Missing credentials are caught so the rest of the notebook still runs.

In [14]:

try:
    nd18, nd22 = mam_ndvi(2018), mam_ndvi(2022)
except Exception as exc:
    nd18 = nd22 = None
    print("NDVI skipped (needs GEE credentials):", exc)

try: nd18, nd22 = mam_ndvi(2018), mam_ndvi(2022) except Exception as exc: nd18 = nd22 = None print("NDVI skipped (needs GEE credentials):", exc)

Map the NDVI change¶

If both seasons loaded, align them and render 2018 greenness, 2022 greenness, and the change — the browning under the drought.

In [15]:

if nd18 is not None and nd22 is not None:
    h, w = min(nd18.shape[0], nd22.shape[0]), min(nd18.shape[1], nd22.shape[1])
    nd18, nd22 = nd18[:h, :w], nd22[:h, :w]

    fig, ax = plt.subplots(1, 3, figsize=(16, 5))
    ax[0].imshow(nd18, cmap="YlGn", vmin=0, vmax=0.6)
    ax[0].set_title("NDVI — MAM 2018")
    ax[1].imshow(nd22, cmap="YlGn", vmin=0, vmax=0.6)
    ax[1].set_title("NDVI — MAM 2022")
    im = ax[2].imshow(nd22 - nd18, cmap="BrBG", vmin=-0.2, vmax=0.2)
    ax[2].set_title("NDVI change (2022 − 2018)")
    for a in ax:
        a.axis("off")
    fig.colorbar(im, ax=ax[2], fraction=0.046)
    plt.tight_layout()
    plt.show()
    print(
        f"mean NDVI  2018 {np.nanmean(nd18):.3f} -> 2022 {np.nanmean(nd22):.3f} "
        f"({(np.nanmean(nd22)-np.nanmean(nd18)):+.3f}) - vegetation browning under the drought"
    )
else:
    print("no NDVI - set GEE credentials and rerun")

if nd18 is not None and nd22 is not None: h, w = min(nd18.shape[0], nd22.shape[0]), min(nd18.shape[1], nd22.shape[1]) nd18, nd22 = nd18[:h, :w], nd22[:h, :w] fig, ax = plt.subplots(1, 3, figsize=(16, 5)) ax[0].imshow(nd18, cmap="YlGn", vmin=0, vmax=0.6) ax[0].set_title("NDVI — MAM 2018") ax[1].imshow(nd22, cmap="YlGn", vmin=0, vmax=0.6) ax[1].set_title("NDVI — MAM 2022") im = ax[2].imshow(nd22 - nd18, cmap="BrBG", vmin=-0.2, vmax=0.2) ax[2].set_title("NDVI change (2022 − 2018)") for a in ax: a.axis("off") fig.colorbar(im, ax=ax[2], fraction=0.046) plt.tight_layout() plt.show() print( f"mean NDVI 2018 {np.nanmean(nd18):.3f} -> 2022 {np.nanmean(nd22):.3f} " f"({(np.nanmean(nd22)-np.nanmean(nd18)):+.3f}) - vegetation browning under the drought" ) else: print("no NDVI - set GEE credentials and rerun")

mean NDVI  2018 0.187 -> 2022 0.144 (-0.043) - vegetation browning under the drought

Recap¶

A 14-year CHIRPS time series turns the drought into numbers: the standardized anomaly flags the consecutive failed long rains, and the 2018-vs-2022 maps show where the deficit fell hardest. The MODIS NDVI cross-check confirms the rainfall deficit as on-the-ground browning. This standardized-anomaly approach is the core of the SPI, the most widely used meteorological-drought index.

Try it yourself¶

• Add the October–December short rains for the full picture of consecutive failures.
• Extend the baseline to 1981 (CHIRPS' start) for a true climatological SPI.
• See the CHC backend reference.