Showcase — the Rhine & Danube low-flow crisis of 2026, from four angles¶
In the summer of 2026 Europe's two great rivers fell to historic lows:
- The Rhine at Kaub — the classic bottleneck gauge — dropped to a navigable depth of ~26–29 cm, brushing the 25 cm record of 2018, and it happened in mid-July, weeks before the usual low-water season (Freight Perspectives).
- The Danube at Budapest hit a record-low 23 cm; cargo shipping nearly stopped and long-sunk WWII wrecks re-emerged near Novi Sad (NBC).
- The Rhine at Lobith recorded 771 m³/s on 15 July 2026 — the lowest July flow on record (Deltares).
We reconstruct the event through earthlens from four complementary angles.
What this notebook does¶
- The flow collapse — GloFAS modeled river discharge at the Kaub (Rhine) and Budapest (Danube)
gauges, 2026 against the 2023–25 baseline (
ecmwf). - Drought spreading — the Copernicus EDO Combined Drought Indicator animated over both basins
(
drought). - From orbit — Sentinel-2 imagery of the Danube at Budapest: before/after, an NDWI water-extent
animation, and an exposed-riverbed map (
gee). - The driver — ERA5-Land cumulative precipitation, 2026 vs recent years, over the basins (
gee).
One deliberate choice: the flow story (1, 2) sits on the Rhine, whose Kaub/Lobith reaches are leveed, dredged shipping channels — their discharge collapses but their width barely changes. The visual story (3) sits on the Danube at Budapest, an un-engineered reach where falling water actually exposes sand and shoals. (CHIRPS is not used — it stops at 50°N, north of which the Rhine sits.)
Runs against the
feat/ecmwf-glofas-historicalbranch. Needs the CDS token (GloFAS) and a GEE service account (Sentinel-2); the EDO layers need no credentials.
Setup¶
First the imports. pyramids reads the GeoTIFFs and NetCDFs (Dataset / DatasetCollection /
NetCDF); earthlens provides the unified EarthLens entry point; the rest are plotting and array
helpers.
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.netcdf import NetCDF
from earthlens.core import EarthLens
warnings.filterwarnings("ignore")
logger.remove()
plt.rcParams["figure.dpi"] = 80
Credentials come from a repo-root .env (loaded with python-dotenv); all downloads land under one
output directory.
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"]
OUT = Path("out") / "rhine_danube_lowflow_2026"
OUT.mkdir(parents=True, exist_ok=True)
The two river gauges we study, and the box covering both basins that the EDO maps are clipped to.
KAUB = (50.09, 7.77) # Rhine bottleneck gauge
BUDAPEST = (47.50, 19.05) # Danube, record-low 23 cm in 2026
GAUGES = {"Rhine-Kaub": KAUB, "Danube-Budapest": BUDAPEST}
BASINS = {"lat_lim": [43.0, 52.0], "lon_lim": [4.0, 22.0]}
1 · The flow collapse — GloFAS river discharge¶
GloFAS is the Copernicus modeled river-discharge reanalysis (m³/s). We compare 2026 with the 2023–25
baseline at each gauge. Two streams feed this: the consolidated reanalysis (final, ~2–3 months
behind) for the 2023–25 baseline, and the intermediate stream (ERA5T, ~2–5 days behind) for 2026,
which reaches the recent weeks the consolidated product has not caught up to. Both carry the NetCDF
variable avg_dis.
The request settings: the discharge variable, the baseline years, and a small ±0.15° box around each gauge (GloFAS is on a 0.05° grid, so that is a handful of pixels).
GLOFAS_VAR = "average-river-discharge-in-the-last-24-hours"
HALF = 0.15 # half-width of the gauge box, in degrees (GloFAS is on a 0.05 deg grid)
GLOFAS_YEARS = [2023, 2024, 2025, 2026] # one Jun-Jul window per year
Download one June–July window per year at each gauge. Earlier years come from the final
consolidated reanalysis; 2026 (which the consolidated stream has not caught up to) comes from the
intermediate stream and is capped at 15 July for its ~2-week latency. Each retrieve is cached under
out/, so re-running reuses the NetCDFs instead of re-hitting the EWDS queue. (One request per year —
not a single multi-year span — so the CDS request stays a clean Jun–Jul selection.)
glofas_paths = {}
for gauge, (lat, lon) in GAUGES.items():
box = {"lat_lim": [lat - HALF, lat + HALF], "lon_lim": [lon - HALF, lon + HALF]}
for year in GLOFAS_YEARS:
dataset = (
"cems-glofas-historical-intermediate"
if year == 2026
else "cems-glofas-historical"
)
end = "2026-07-15" if year == 2026 else f"{year}-07-31"
ddir = OUT / "glofas" / f"{gauge}_{year}"
cached = sorted(ddir.glob("*.nc"))
if cached:
glofas_paths[(gauge, year)] = cached
continue
ddir.mkdir(parents=True, exist_ok=True)
glofas_paths[(gauge, year)] = EarthLens(
data_source="ecmwf",
variables={dataset: [GLOFAS_VAR]},
start=f"{year}-06-01",
end=end,
temporal_resolution="daily",
path=str(ddir),
skip_constraints=True, # validator false-negatives on multi-month consolidated
**box,
).download()
{k: len(v) for k, v in glofas_paths.items()}
Each NetCDF holds avg_dis on a (time, lat, lon) grid. Averaging the gauge box over space and
over the June–July window gives one mean discharge per (gauge, year).
records = []
for (gauge, year), paths in glofas_paths.items():
for path in paths:
nc = NetCDF.read_file(str(path))
arr = np.asarray(nc.read_array("avg_dis"), dtype="float32") # (time, y, x)
nc.close()
valid = arr[arr > 0]
mean = float(np.nanmean(valid)) if valid.size else np.nan
records.append({"gauge": gauge, "year": year, "discharge": mean})
flow = pd.DataFrame(records)
flow.pivot_table(index="year", columns="gauge", values="discharge").round(1)
Plot each gauge's June–July mean discharge by year, with 2026 highlighted. If 2026 is the driest, its bar sits clearly below the baseline years.
fig, axes = plt.subplots(1, 2, figsize=(13, 4.5))
for ax, gauge in zip(axes, GAUGES):
g = flow[flow.gauge == gauge].sort_values("year")
colors = ["#c1440e" if y == 2026 else "#4a90d9" for y in g.year]
ax.bar(g.year.astype(str), g.discharge, color=colors)
ax.set(title=f"{gauge} - mean Jun-Jul discharge", ylabel="discharge (m^3/s)")
ax.tick_params(axis="x", rotation=45)
plt.tight_layout()
plt.show()
2 · Drought spreading — Copernicus EDO Combined Drought Indicator¶
The European Drought Observatory's dedicated Low-Flow Index (edo-lfinx-lgs) is the ideal layer
here, but it is not published for 2026 (its WCS returns a server error for 2026 dates — verified
live). So we use the Combined Drought Indicator (edo-cdiad), which does reach 2026: a categorical
severity index blending precipitation, soil-moisture and vegetation anomalies, on a 10-day cadence and
free to use. We animate it over both basins across the drought's progression.
Download the Combined Drought Indicator over the two-basin box for April–June 2026 (cached under
out/). CDI is not real-time — verified live, it is published through the 21 June 2026 dekad, so the end
date is capped there.
cdi_dir = OUT / "edo_cdi"
cdi_cached = sorted(cdi_dir.glob("*.tif"))
if cdi_cached:
cdi_paths = cdi_cached
else:
cdi_paths = EarthLens(
data_source="drought",
dataset="edo-cdiad",
variables=[],
start="2026-04-01",
end="2026-06-21",
path=str(cdi_dir),
**BASINS,
).download(progress_bar=False)
len(cdi_paths)
Read the dekads into one time-stacked collection and derive a single colour scale from the whole stack, so every animation frame is directly comparable.
dc_cdi = DatasetCollection.read_multiple_files(str(cdi_dir), date=False)
cdi_frames = [
np.asarray(dc_cdi.iloc(i).read_array(), "float32")
for i in range(dc_cdi.time_length)
]
cdi_stack = np.stack(cdi_frames)
cdi_vmin, cdi_vmax = float(np.nanmin(cdi_stack)), float(np.nanmax(cdi_stack))
cdi_vmin, cdi_vmax
Animate the drought severity spreading across the Rhine and Danube basins.
cdi_gif = OUT / "edo_cdi.gif"
dc_cdi.plot(
band=0,
cmap="YlOrRd",
vmin=cdi_vmin,
vmax=cdi_vmax,
figsize=(6, 5),
title="Rhine + Danube basins - EDO Combined Drought Indicator, Apr-Jun 2026",
).save_animation(str(cdi_gif), fps=1.5)
plt.close("all")
display(Image(filename=str(cdi_gif)))
3 · From orbit — the Danube at Budapest¶
The Danube through Budapest is un-leveed enough that falling flow exposes real sand and shoals. We pull a cloud-median Sentinel-2 composite per month through the drought and derive NDWI (open water).
The reach and the monthly windows. Windows are ~4 weeks wide so the per-pixel median has enough looks to reject cloud (Central-European summers are not reliably clear). The pixel size is 20 m — a compromise that keeps each true-colour composite under Earth Engine's 50 MB single-URL export cap while still resolving the channel and its sandbars.
DANUBE_AOI = [18.98, 47.40, 19.14, 47.58] # lon_min, lat_min, lon_max, lat_max
S2 = "COPERNICUS/S2_SR_HARMONIZED"
S2_SCALE = 20 # metres per pixel (keeps the GEE URL export under 50 MB)
WINDOWS = [
("2026-04-01", "2026-04-28"),
("2026-05-01", "2026-05-28"),
("2026-06-01", "2026-06-28"),
("2026-07-01", "2026-07-24"),
]
rgb_dir = OUT / "danube_rgb"
ndwi_dir = OUT / "danube_ndwi"
rgb_dir.mkdir(parents=True, exist_ok=True)
ndwi_dir.mkdir(parents=True, exist_ok=True)
For each window, fetch a true-colour (+NIR) median composite, save it, and derive an NDWI raster
(green − NIR)/(green + NIR). We also record each composite's mean visible brightness — a cloud proxy
used in the next cell. Everything is cached, so re-runs skip the Earth Engine pulls.
labels = []
visible_means = []
for start, end in WINDOWS:
tag = start.replace("-", "")
rgb_path = rgb_dir / f"budapest_{tag}.tif"
if not rgb_path.exists():
scene = EarthLens(
data_source="gee",
cadence="raw",
path=tempfile.mkdtemp(),
dataset=S2,
variables=["B4", "B3", "B2", "B8"],
aoi=DANUBE_AOI,
start=start,
end=end,
scale=S2_SCALE,
reducer="median",
)
scene.authenticate(service_account=SERVICE_ACCOUNT, service_key=SERVICE_KEY)
outputs = scene.download(progress_bar=False)
if not outputs:
continue
shutil.copy(str(outputs[0]), rgb_path)
a = np.asarray(Dataset.read_file(str(rgb_path)).read_array(), dtype="float32")
ndwi = (a[1] - a[3]) / (a[1] + a[3] + 1e-6)
src = Dataset.read_file(str(rgb_path))
Dataset.create_from_array(
arr=ndwi[None, :, :], geo=src.geotransform, epsg=src.epsg
).to_file(str(ndwi_dir / f"ndwi_{tag}.tif"))
labels.append(start)
visible_means.append(float(a[:3].mean()))
labels
Flag cloud-contaminated composites: a cloudy median stays bright and flat-spectrum, so any window more than 1.5x brighter than the clearest one is dropped from the quantitative panels (it still appears in the animation). This catches a hazy spring composite that would otherwise corrupt the NDWI.
visible_means = np.array(visible_means)
clear = visible_means <= 1.5 * visible_means.min()
clear_idx = [i for i, ok in enumerate(clear) if ok]
dc_rgb = DatasetCollection.read_multiple_files(str(rgb_dir), date=False)
dc_ndwi = DatasetCollection.read_multiple_files(str(ndwi_dir), date=False)
pd.DataFrame({"window": labels, "mean_visible": visible_means.round(0), "clear": clear})
Before / after true colour — the first vs last cloud-clear window.
if len(clear_idx) >= 2:
first_i, last_i = clear_idx[0], clear_idx[-1]
fig, axes = plt.subplots(1, 2, figsize=(11, 5.5))
rgb_opts = {"rgb": [0, 1, 2], "percentile": 2}
for ax, i in zip(axes, (first_i, last_i)):
dc_rgb.iloc(i).plot(rgb_options=rgb_opts, ax=ax, fig=fig)
ax.set_title(labels[i])
ax.set_xticks([])
ax.set_yticks([])
fig.suptitle("Danube at Budapest - Sentinel-2 true colour, before vs after", y=1.02)
plt.tight_layout()
plt.show()
else:
print("not enough cloud-clear windows for a before/after panel")
Animate the NDWI water extent across all windows.
ndwi_gif = OUT / "danube_ndwi.gif"
dc_ndwi.plot(
band=0,
cmap="Blues",
figsize=(5, 5),
title="Danube at Budapest - NDWI (open water), Apr-Jul 2026",
).save_animation(str(ndwi_gif), fps=1)
plt.close("all")
display(Image(filename=str(ndwi_gif)))
Turn the NDWI into a number: count open-water pixels (each S2_SCALE m across) over the
cloud-clear windows to get the channel's wetted area per window.
px_km2 = (S2_SCALE * S2_SCALE) / 1e6
clear_labels = [labels[i] for i in clear_idx]
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)
areas.append(float(np.sum(ndwi > 0) * px_km2))
fig, ax = plt.subplots(figsize=(7, 4))
ax.bar(clear_labels, areas, color="tab:blue")
for j, v in enumerate(areas):
ax.text(j, v, f"{v:.2f}", ha="center", va="bottom")
ax.set(ylabel="wetted area (km^2)", title="Danube at Budapest - wetted area")
ax.tick_params(axis="x", rotation=30)
plt.tight_layout()
plt.show()
dict(zip(clear_labels, [round(v, 2) for v in areas]))
Reading this honestly. The Danube through Budapest is a wide, embanked urban reach, so its open-water width is far less sensitive to discharge than an un-engineered braided river — the wetted area moves only modestly, and the numbers carry Sentinel-2 mixed-pixel and residual-haze noise. Read it as a trend (the drop into the record-low July), not a precise gauge. The definitive flow evidence is the GloFAS discharge in section 1.
4 · The driver — precipitation deficit (ERA5-Land)¶
The missing rain behind the low flow. CHIRPS stops at 50°N (north of much of the Rhine), so we use ERA5-Land daily precipitation (hosted on Earth Engine) and compare 2026 against recent years over the two basins. If 2026 is genuinely drier, its cumulative-precipitation curve sits below the others.
Request settings: four years of daily ERA5-Land precipitation over the basins. ERA5-Land's daily aggregate lags real time by ~2 weeks, so 2026 is capped at 15 July; the baseline years run the full 1 April – 15 July window.
precip_dir = OUT / "era5_precip"
PRECIP_YEARS = [2023, 2024, 2025, 2026]
LAST_2026 = "2026-07-15" # latest published ERA5-Land daily day
Download one daily precipitation tile per day per year over the basins (cached under out/).
era5 = {}
for year in PRECIP_YEARS:
ydir = precip_dir / str(year)
cached = sorted(ydir.glob("*.tif"))
if cached:
era5[year] = cached
continue
ydir.mkdir(parents=True, exist_ok=True)
end = LAST_2026 if year == 2026 else f"{year}-07-15"
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=11000.0,
path=ydir,
export_via="url",
**BASINS,
)
job.authenticate(service_account=SERVICE_ACCOUNT, service_key=SERVICE_KEY)
era5[year] = job.download(progress_bar=False)
{y: len(p) for y, p in era5.items()}
Reduce each daily tile to a basin-mean precipitation (mm), indexed by day-of-season (days since 1 April) so the four years line up regardless of calendar date.
rows = []
for year, paths in era5.items():
for p in sorted(paths):
day = pd.to_datetime(Path(p).stem.split("_")[-1])
arr = np.asarray(Dataset.read_file(str(p)).read_array(), dtype="float32")
rows.append(
{
"year": year,
"season_day": (day - pd.Timestamp(f"{year}-04-01")).days,
"precip_mm": float(np.nanmean(arr)) * 1000.0,
}
)
precip_df = pd.DataFrame(rows).sort_values(["year", "season_day"])
precip_df.groupby("year")["precip_mm"].sum().round(1)
Plot cumulative precipitation per year, 2026 highlighted — a curve below the others is a drier season.
fig, ax = plt.subplots(figsize=(9, 5))
totals = {}
for year, g in precip_df.groupby("year"):
g = g.sort_values("season_day")
cum = g["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(g["season_day"], cum, label=str(year), **style)
ax.set(
xlabel="days since 1 April",
ylabel="cumulative precipitation (mm)",
title="Rhine + Danube basins - cumulative ERA5-Land precipitation",
)
ax.legend(title="year")
ax.grid(alpha=0.3)
plt.tight_layout()
plt.show()
{y: round(t, 1) for y, t in totals.items()}
Recap¶
Four independent lines of evidence on one event:
ecmwf→ GloFAS discharge — the modeled flow collapse at Kaub and Budapest, 2026 vs the 2023–25 baseline.drought→ EDO Combined Drought Indicator — drought severity spreading across both basins.gee→ Sentinel-2 — the Danube at Budapest before/after, with a wetted-area trend into July.gee→ ERA5-Land precipitation — 2026's cumulative rainfall against recent years, the deficit driver.
Sources¶
- Freight Perspectives — Rhine at Kaub 2026
- NBC — record-low Danube
- Deltares — Rhine/Meuse low water 2026
- ECMWF — GloFAS river discharge
Try it yourself¶
- Add the Danube at Novi Sad (
45.25, 19.85) or the Waal near Nijmegen as a second optical reach. - Swap
edo-lfinx-lgsforgdo-lfinx-sms/-mds, or add GRACE terrestrial-water-storage viagee.