Showcase — the Netherlands' 2026 water shortage, from orbit and from the ground¶
On 24 July 2026 the Netherlands officially declared a water shortage, raising the national drought response to Level 2. The signals behind that declaration are real and dramatic:
- The Rhine at Lobith — the gauge where the river enters the Netherlands — recorded a daily-average flow of 771 m³/s on 15 July 2026, the lowest July flow ever measured, rivalling the historic 1976 drought year (Deltares).
- National river levels and the precipitation deficit are both at "1-in-20-year" severity, according to KNMI and the national water-distribution committee (DutchNews).
- The Twentekanaal (Gelderland) closed to shipping from 25 July 2026 — the Eefde locks were shut to stop irreparable damage to flood defenses, soil and nature, and farmers were banned from irrigating with its water (NL Times).
- The IJsselmeer is deliberately being kept artificially high as the country's freshwater buffer — areas that cannot draw from it are drying out fastest (NL Times).
- Houseboats sit stranded and riverbanks lie exposed along the Rhine and Meuse (CNN).
Two data-source decisions follow directly from those facts:
- CHIRPS/CHIRTS are the wrong tool here — that catalog is quasi-global, 50°S–50°N, and the Netherlands sits at 51–53.5°N, outside the band. Precipitation and soil moisture come from ERA5-Land instead.
- The IJsselmeer and the Twentekanaal are the wrong satellite targets. Both are level-controlled (the IJsselmeer is kept artificially high; the Twentekanaal's locks were shut specifically to stop its level dropping), so a water-extent animation there would show a flat, unconvincing signal. The Rhine at Lobith — the exact reach the 771 m³/s record refers to — is at least un-leveed and worth testing, even though (spoiler, see Section 3) being a maintained shipping channel it turns out to be a subtler signal than expected too.
What this notebook does¶
- The official drought signal — Copernicus European Drought Observatory (EDO) Combined
Drought Indicator, animated over the Netherlands from April to June 2026 (
droughtbackend, no credentials). - Precipitation deficit, quantified — daily ERA5-Land precipitation over the same area, for
four consecutive years (2023–2026), so "lowest in years" is a number pulled from data rather
than a quote from the news (
geebackend). - Seeing it from orbit — Sentinel-2 true colour and NDWI over the Rhine at Lobith across the
same season: a before/after comparison, an animated water-extent GIF, a wetted-area chart, and an
honest read of what that chart does (and doesn't) show (
geebackend).
Needs
GEE_SERVICE_ACCOUNT/GEE_SERVICE_KEYandpyramids-gis[viz]. Thedroughtbackend needs no credentials.
Setup¶
pyramids supplies Dataset / DatasetCollection for reading, plotting and animating rasters;
earthlens supplies the unified EarthLens entry point.
import os
import shutil
import tempfile
import warnings
from pathlib import Path
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from IPython.display import Image, display
from loguru import logger
from pyramids.dataset import Dataset, DatasetCollection
from pyramids.plot import ColorBar
from earthlens.core import EarthLens
warnings.filterwarnings("ignore")
logger.remove()
plt.rcParams["figure.dpi"] = 80
OUT = Path("out") / "netherlands_drought_2026"
OUT.mkdir(parents=True, exist_ok=True)
Credentials¶
Read from a .env file at the repository root (loaded with python-dotenv; your PATH is left
untouched):
GEE_SERVICE_ACCOUNT="your-sa@your-project.iam.gserviceaccount.com"
GEE_SERVICE_KEY="/path/to/your/service-account-key.json"
from dotenv import find_dotenv, load_dotenv
load_dotenv(find_dotenv(usecwd=True))
SERVICE_ACCOUNT = os.environ["GEE_SERVICE_ACCOUNT"]
SERVICE_KEY = os.environ["GEE_SERVICE_KEY"]
1 · The official drought signal — Copernicus EDO Combined Drought Indicator¶
The drought backend reaches the Copernicus European Drought Observatory over its GetCoverage WCS
— no credentials, no discovery handshake. edo-cdiad is the Combined Drought Indicator, a
categorical severity index blending precipitation, soil moisture and vegetation anomalies. A single
start/end request over the whole drought's 10-day-release cadence returns one GeoTIFF per dekad
covering the country.
CDI is not real-time: verified live against the Copernicus endpoint, as of this writing it is
published through the 21 June 2026 dekad (a DATE_OUT_OF_RANGE error for anything past that), so
the animation below stops there — a few weeks before the 24 July water-shortage declaration, but
already showing the drought that led to it. Section 3 picks the story back up with Sentinel-2 imagery
all the way to 24 July, where the satellite record has no such lag.
NL_LAT, NL_LON = [50.7, 53.6], [3.3, 7.3]
edo_dir = OUT / "edo_cdi"
TIF_GLOB = "*.tif" # reused by the ERA5-Land and CDI cache checks below
edo_cached = sorted(edo_dir.glob(TIF_GLOB))
if edo_cached: # reuse a previous run's rasters instead of re-hitting Copernicus
edo_paths = edo_cached
else:
edo_paths = EarthLens(
data_source="drought",
dataset="edo-cdiad",
variables=[],
start="2026-04-01",
end="2026-06-21", # CDI's most recent published dekad as of this writing
lat_lim=NL_LAT,
lon_lim=NL_LON,
path=str(edo_dir),
).download(progress_bar=False)
len(edo_paths)
Animate the dekad-by-dekad progression¶
The colour stretch is taken from the data's own min/max across the whole stack (not a guessed range), so every frame shares one consistent scale.
dc_cdi = DatasetCollection.from_files(str(edo_dir))
cdi_labels = [p.stem.split("_")[-1] for p in sorted(edo_dir.glob(TIF_GLOB))]
cdi_stack = np.stack(
[
np.asarray(dc_cdi.iloc(i).read_array(), dtype="float32")
for i in range(dc_cdi.time_length)
]
)
cdi_vmin, cdi_vmax = float(np.nanmin(cdi_stack)), float(np.nanmax(cdi_stack))
cdi_vmin, cdi_vmax
gif = OUT / "edo_cdi.gif"
dc_cdi.plot(
band=0,
cmap="YlOrRd",
vmin=cdi_vmin,
vmax=cdi_vmax,
figsize=(5, 5),
title="Netherlands - Copernicus EDO Combined Drought Indicator, Apr-Jun 2026",
).save_animation(str(gif), fps=1.5)
plt.close("all")
display(Image(filename=str(gif)))
The worst dekad, as a single hero map¶
The most recent (and most severe) dekad on its own, with a labelled colorbar.
fig, ax = plt.subplots(figsize=(7, 6))
dc_cdi.iloc(dc_cdi.time_length - 1).plot(
band=0,
ax=ax,
fig=fig,
cmap="YlOrRd",
vmin=cdi_vmin,
vmax=cdi_vmax,
colorbar=ColorBar(label="Combined Drought Indicator"),
)
ax.set_title(f"Netherlands drought severity - dekad of {cdi_labels[-1]}")
ax.set_xlabel("lon")
ax.set_ylabel("lat")
plt.tight_layout()
plt.show()
2 · Precipitation deficit, quantified — ERA5-Land via Earth Engine¶
ERA5-Land's daily aggregates are hosted on Earth Engine (ECMWF/ERA5_LAND/DAILY_AGGR), so we reach
them through the gee key and sidestep the CDS queue entirely — the same trick the heat-wave
notebook uses. To turn "lowest in years" into a number
instead of a quote, we pull the same 1 April window across four consecutive years (2023–2026) and
compare cumulative precipitation directly. (This daily-aggregate collection carries
total_precipitation_sum but not a soil-moisture band — ECMWF/ERA5_LAND/HOURLY has
volumetric_soil_water_layer_1 for anyone who wants to add that angle.)
ERA5-Land's daily aggregate also isn't instantaneous: verified live, the collection's most recent image as of this writing is 15 July 2026, about nine days behind. So 2026 stops at 15 July while 2023–2025 (fully in the past) run the complete window to 24 July — still a fair comparison, since we only ever compare each year's cumulative total at the same day-of-season.
precip_dir = OUT / "era5_land"
YEARS = [2023, 2024, 2025, 2026]
LAST_AVAILABLE_2026 = "2026-07-15" # verified live against ECMWF/ERA5_LAND/DAILY_AGGR
era5_daily = {}
for year in YEARS:
year_dir = precip_dir / str(year)
cached = sorted(year_dir.glob(TIF_GLOB))
if cached: # reuse a previous run's tiles instead of re-pulling from Earth Engine
era5_daily[year] = cached
print(year, len(cached), "daily tiles (cached)")
continue
end = LAST_AVAILABLE_2026 if year == 2026 else f"{year}-07-24"
job = EarthLens(
data_source="gee",
dataset="ECMWF/ERA5_LAND/DAILY_AGGR",
variables=["total_precipitation_sum"],
start=f"{year}-04-01",
end=end,
temporal_resolution="daily",
scale=9000.0,
path=year_dir,
export_via="url",
lat_lim=NL_LAT,
lon_lim=NL_LON,
)
job.authenticate(service_account=SERVICE_ACCOUNT, service_key=SERVICE_KEY)
era5_daily[year] = job.download(progress_bar=False)
print(year, len(era5_daily[year]), "daily tiles")
Reduce each daily tile to a country-mean¶
Average each day's precipitation tile spatially over the NL box, and index by day-of-season (days since 1 April) so all four years line up regardless of the calendar date.
records = []
for year, paths in era5_daily.items():
for p in sorted(paths):
day_str = p.stem.split("_")[-1]
day = pd.to_datetime(day_str)
season_day = (day - pd.Timestamp(f"{year}-04-01")).days
arr = np.asarray(Dataset.read_file(p).read_array(), dtype="float32")
records.append(
{
"year": year,
"season_day": season_day,
"precip_mm": float(np.nanmean(arr)) * 1000.0,
}
)
era5_df = pd.DataFrame(records).sort_values(["year", "season_day"])
era5_df.groupby("year")["precip_mm"].mean().round(3)
Cumulative precipitation, four years side by side¶
Compared at the same day-of-season, 2026 and 2025 land almost on top of each other — 138.0 mm and 136.3 mm — while 2023 reaches 153.5 mm and 2024 reaches 262.6 mm. So the honest reading is not "2026 is the driest year" but "2026 is the second of two consecutive springs at roughly half of 2024", with 2025 marginally lower still. That is what makes the 2026 declaration notable: the shortage arrived on top of an already-depleted 2025, not after a wet year.
fig, ax = plt.subplots(figsize=(9, 5))
# 2026 stops at the last available ERA5-Land day, so its season is shorter than
# the others'. Compare every year at the same day-of-season, otherwise 2026 looks
# drier simply for having fewer days in the sum.
common_last = int(era5_df.groupby("year")["season_day"].max().min())
totals = {}
for year, grp in era5_df.groupby("year"):
grp = grp.sort_values("season_day")
# Truncate before the cumulative sum so the curve and the total describe the
# same window as the title.
grp = grp[grp["season_day"] <= common_last]
cum = grp["precip_mm"].cumsum()
totals[year] = float(cum.iloc[-1])
style = {"lw": 3, "color": "#c1440e"} if year == 2026 else {"lw": 1.5, "alpha": 0.6}
ax.plot(grp["season_day"], cum, label=str(year), **style)
ax.set_xlabel("days since 1 April")
ax.set_ylabel("cumulative precipitation (mm)")
ax.set_title(
f"Netherlands - cumulative ERA5-Land precipitation, "
f"1 Apr to day {common_last} of season"
)
ax.legend(title="year")
ax.grid(alpha=0.3)
plt.tight_layout()
plt.show()
{year: round(total, 1) for year, total in totals.items()}
3 · Seeing it from orbit — the Rhine at Lobith¶
Lobith is the exact gauge behind the 771 m³/s record, and — unlike the IJsselmeer or the
Twentekanaal — a free-flowing reach with no locks or level control, so it's a reasonable place to
test whether low discharge shows up as a visibly narrower channel. We pull a near-monthly,
cloud-median Sentinel-2 composite (true colour B4,B3,B2 + NIR B8) through the drought's
progression and derive NDWI (open water) from each. The windows are widened to ~4 weeks (rather than
2) specifically to give the per-pixel median enough independent looks to reject cloud — the Dutch
spring/summer sky is not reliably clear even during a precipitation drought, and a narrower window
left two of five composites visibly cloud-contaminated in an earlier pass of this notebook.
LOBITH_AOI = [6.03, 51.82, 6.18, 51.89] # lon_min, lat_min, lon_max, lat_max
S2 = "COPERNICUS/S2_SR_HARMONIZED"
WINDOWS = [
("2026-04-01", "2026-04-28"),
("2026-05-01", "2026-05-28"),
("2026-06-01", "2026-06-28"),
("2026-07-01", "2026-07-15"),
("2026-07-16", "2026-07-24"),
]
rgb_dir, ndwi_dir = OUT / "lobith_rgb", OUT / "lobith_ndwi"
rgb_dir.mkdir(parents=True, exist_ok=True)
ndwi_dir.mkdir(parents=True, exist_ok=True)
window_labels = []
visible_means = [] # mean B4/B3/B2 reflectance per window - a cloud-contamination proxy
for start, end in WINDOWS:
tag = start.replace("-", "")
rgb_path = rgb_dir / f"lobith_{tag}.tif"
if (
rgb_path.exists()
): # reuse a previous run's composite; skip the Earth Engine pull
p = str(rgb_path)
else:
scene = EarthLens(
data_source="gee",
cadence="raw",
path=tempfile.mkdtemp(), # scratch only; the OS reclaims it, nothing here is kept
dataset=S2,
variables=["B4", "B3", "B2", "B8"],
aoi=LOBITH_AOI,
start=start,
end=end,
scale=10,
reducer="median",
)
scene.authenticate(service_account=SERVICE_ACCOUNT, service_key=SERVICE_KEY)
outputs = scene.download(progress_bar=False)
if not outputs:
continue
p = str(outputs[0])
shutil.copy(p, rgb_path)
a = np.asarray(Dataset.read_file(p).read_array(), dtype="float32")
ndwi = (a[1] - a[3]) / (a[1] + a[3] + 1e-6) # (Green - NIR)/(Green + NIR)
src = Dataset.read_file(p)
ndwi_ds = Dataset.create_from_array(
arr=ndwi[None, :, :], geo=src.geotransform, epsg=src.epsg
)
ndwi_ds.to_file(str(ndwi_dir / f"ndwi_{tag}.tif"))
window_labels.append(start)
visible_means.append(float(a[:3].mean()))
len(window_labels), window_labels
Flag cloud-contaminated composites¶
A median composite over a mostly-overcast window stays bright and flat-spectrum (clouds reflect strongly across red/green/blue alike), unlike this reach's actual land/water mix. Rather than trust every window equally, flag any whose mean visible reflectance is far above the group's clearest window — those get excluded from the quantitative panels below, though they still appear (visibly cloudy) in the true-colour animation.
visible_means = np.array(visible_means)
CLEAR = visible_means <= 2.0 * visible_means.min()
clear_idx = [i for i, ok in enumerate(CLEAR) if ok]
pd.DataFrame(
{
"window": window_labels,
"mean_visible_reflectance": visible_means.round(0),
"clear": CLEAR,
}
)
dc_rgb = DatasetCollection.from_files(str(rgb_dir)) if window_labels else None
dc_ndwi = DatasetCollection.from_files(str(ndwi_dir)) if window_labels else None
print(f"{len(window_labels)} usable windows:", window_labels)
print(f"{len(clear_idx)} cloud-clear windows:", [window_labels[i] for i in clear_idx])
Before / after — true colour¶
The first clear window (spring, before the worst of the drought) against the last clear one — whichever of the flagged-clean windows is latest — rather than a fixed calendar date that might land on a cloudy composite.
if dc_rgb and len(clear_idx) >= 2:
first_i, last_i = clear_idx[0], clear_idx[-1]
fig, axes = plt.subplots(1, 2, figsize=(11, 5.5))
for ax, i, tag in zip(
axes, (first_i, last_i), (window_labels[first_i], window_labels[last_i])
):
dc_rgb.iloc(i).plot(
rgb_options={"rgb": [0, 1, 2], "percentile": 2}, ax=ax, fig=fig
)
ax.set_title(tag)
ax.set_xticks([])
ax.set_yticks([])
fig.suptitle("Rhine at Lobith - Sentinel-2 true colour, before vs after", y=1.02)
plt.tight_layout()
plt.show()
else:
print("not enough cloud-free windows for a before/after panel")
Animated true colour and animated NDWI¶
The same pattern as the GERD reservoir-filling showcase: step through the composites and watch for any visible change in the channel's width or the exposure of its sandbanks.
if dc_rgb:
gif = OUT / "lobith_true_colour.gif"
dc_rgb.plot(
rgb_options={"rgb": [0, 1, 2], "percentile": 2},
figsize=(5, 5),
title="Rhine at Lobith - Sentinel-2 true colour, Apr-Jul 2026",
).save_animation(str(gif), fps=1)
plt.close("all")
display(Image(filename=str(gif)))
else:
print("no scenes to animate")
if dc_ndwi:
gif = OUT / "lobith_ndwi.gif"
dc_ndwi.plot(
band=0,
cmap="Blues",
figsize=(5, 5),
title="Rhine at Lobith - NDWI (open water), Apr-Jul 2026",
).save_animation(str(gif), fps=1)
plt.close("all")
display(Image(filename=str(gif)))
else:
print("no scenes to animate")
Wetted area per window¶
Counting NDWI > 0 pixels turns the animation into a number — the channel's open-water area at each window, scaled by the 10 m x 10 m Sentinel-2 pixel. Only the cloud-clear windows flagged above are charted — a cloud-contaminated composite's NDWI count is not a real measurement of the channel.
if dc_rgb and clear_idx:
px_km2 = (10 * 10) / 1e6
wet_areas = []
for i in clear_idx:
a = np.asarray(dc_rgb.iloc(i).read_array(), dtype="float32")
ndwi = (a[1] - a[3]) / (a[1] + a[3] + 1e-6)
wet_areas.append(float(np.sum(ndwi > 0) * px_km2))
clear_labels = [window_labels[i] for i in clear_idx]
fig, ax = plt.subplots(figsize=(8, 4))
ax.bar(clear_labels, wet_areas, color="tab:blue")
for i, v in enumerate(wet_areas):
ax.text(i, v, f"{v:.2f}", ha="center", va="bottom")
ax.set(
ylabel="wetted area (km^2)",
title="Rhine at Lobith - wetted area, cloud-clear windows only",
)
plt.xticks(rotation=30, ha="right")
plt.tight_layout()
plt.show()
print(
"wetted area (km^2):", {k: round(v, 2) for k, v in zip(clear_labels, wet_areas)}
)
else:
print("no cloud-clear windows to measure")
Reading the result honestly¶
Only two windows survived the cloud filter, and between them the wetted area grew rather than shrank — the direction is the point, and the plot above carries the magnitude. That is a real result, not a bug, and it rejects the "the channel visibly narrows" hypothesis for this reach — worth saying plainly rather than forcing the opposite conclusion.
Two things explain it. First, timing: the 771 m³/s record is specifically a July anomaly — Alpine snowmelt still buffers Rhine discharge in spring, so an April-vs-July comparison isn't guaranteed to run in one direction. Second, and more fundamentally: Lobith is the main-stem international shipping channel — deep, leveed and engineered to a maintained navigation depth — so it does not respond to a discharge drop the way an unregulated, braided reach would. A gauge reading of "lowest July flow on record" does not automatically translate into a satellite-visible geomorphology change at this particular location. The EDO drought-index animation and the four-year precipitation comparison above remain the more sensitive instruments for this event; the Sentinel-2 imagery here is best read qualitatively (the pale sand/gravel margins along the banks are real, permanent training works, not a drought artefact) rather than as a precise area trend. A smaller, unregulated reach — the Waal's braided floodplain near Nijmegen, for instance — would be a better bet for a visible signal, and is exactly the kind of thing to try in "Try it yourself" below.
Recap¶
Three independent lines of evidence, three different visualization languages, one event — and one honestly negative result:
drought→ Copernicus EDO Combined Drought Indicator — an animated map of the drought severity spreading over the Netherlands from April to June 2026 (the most recent published dekad).gee→ ERA5-Land — a four-year (2023–2026) comparison turning "lowest in years" into an actual cumulative-precipitation number (2026's 138 mm and 2025's 136 mm against 154 mm and 263 mm in 2023 and 2024, all compared at the same day-of-season), not just a quote.gee→ Sentinel-2 true colour + NDWI — a before/after panel, an animated water-extent GIF, and a wetted-area chart over the Rhine at Lobith, which rejects the assumption that this particular reach visibly narrows — a useful, honest boundary on what medium-resolution optical imagery of a leveed shipping channel can and can't show.
Sources¶
- NL Times — water shortage declared, drought response Level 2
- Deltares — 17 questions: drought and low water levels in Rhine and Meuse
- DutchNews — risk of serious water shortages as river levels fall
- NL Times — Twentekanaal closed to shipping traffic
- NL Times — water shortage looms as drought continues
- CNN — Europe's rivers are shrinking amid heat and drought
Try it yourself¶
- Swap
edo-cdiadforedo-smand(soil-moisture anomaly) oredo-spgTS(monthly SPI) — see the drought backend reference. - Point
LOBITH_AOIat another un-leveed reach (the Waal near Nijmegen, the Meuse near Maastricht) and compare wetted-area curves. - Extend
YEARSfurther back to see how 2026 stacks up against a longer climatology.