Ready — 9 reservoirs available: Beni Haroun — Algeria · Hammam Boughrara (Tlemcen) — Algeria · Beni Bahdel (Tlemcen) — Algeria · Al Massira — Morocco · Lake Nasser — Egypt · Lake Volta — Ghana · Bui — Ghana · King Fahd Dam — Saudi Arabia · Wadi Dayqah — Oman

# @title Step 1 · Choose a reservoir
RESERVOIR = "Beni Haroun — Algeria"  # @param ["Beni Haroun — Algeria", "Hammam Boughrara (Tlemcen) — Algeria", "Beni Bahdel (Tlemcen) — Algeria", "Al Massira — Morocco", "Lake Nasser — Egypt", "Lake Volta — Ghana", "Bui — Ghana", "King Fahd Dam — Saudi Arabia", "Wadi Dayqah — Oman"]
site = SITE_ID[RESERVOIR]
info = fetch_json(f"{site}/summary.json")
catchment = fetch_vector(f"{site}/catchment.geojson")
extent = fetch_vector(f"{site}/reservoir_max_extent.geojson")

facts = [
    ("Country", info["country"]),
    ("Largest water surface seen 1984–2021", f'{info["reservoir_max_extent_km2"]:,.1f} km²'),
    ("Upstream catchment", f'{info["catchment_km2"]:,.0f} km²'),
    ("Mean rainfall over the catchment, 1981–2026", f'{info["rainfall_mean_annual_mm"]:,.0f} mm per year'),
    ("Sentinel-2 monthly record", f'{info["months_reported"]} months, {info["first_year"]} to {info["latest_month"][:4]}'),
]
rows = "".join(f"<tr><td style='padding:3px 16px 3px 0;color:#6B7280'>{k}</td>"
               f"<td style='padding:3px 0;font-weight:600'>{v}</td></tr>" for k, v in facts)
display(HTML(f"<div style='font-family:Inter,Arial,sans-serif'><h3 style='margin:4px 0 8px'>{info['name']}</h3>"
             f"<table>{rows}</table></div>"))

m = folium.Map(tiles=None, control_scale=True, height=520)
folium.TileLayer(tiles="https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{z}/{y}/{x}",
                 attr="Imagery © Esri, Maxar, Earthstar Geographics", name="Satellite imagery").add_to(m)
folium.TileLayer("cartodbpositron", name="Light basemap").add_to(m)
folium.GeoJson(catchment, name="Upstream catchment",
               style_function=lambda f: {"color": "#FFC107", "weight": 2, "fillOpacity": 0.04}).add_to(m)
folium.GeoJson(extent, name="Largest water surface, 1984–2021",
               style_function=lambda f: {"color": "#4FC3F7", "weight": 1.5, "fillColor": "#4FC3F7",
                                         "fillOpacity": 0.35}).add_to(m)
folium.LayerControl(collapsed=False).add_to(m)
b = catchment.total_bounds if info["catchment_km2"] < 60000 else extent.total_bounds
m.fit_bounds([[b[1], b[0]], [b[3], b[2]]])
display(m)
if info["catchment_km2"] >= 60000:
    print("Zoom out to see the whole catchment — it is very large.")

# @title Step 2 · Map today's water (about 30 seconds)
from pystac_client import Client
from odc.stac import load, configure_rio

configure_rio(cloud_defaults=True, aws={"aws_unsigned": True})
grid = info["grid"]
catalog = Client.open("https://earth-search.aws.element84.com/v1")   # public catalogue, no account needed

today = pd.Timestamp.now(tz="UTC").normalize()
items = []
for days, max_cloud in [(45, 10), (90, 20), (180, 30)]:
    found = catalog.search(collections=["sentinel-2-l2a"], bbox=grid["bbox"],
                           datetime=f"{today - pd.Timedelta(days=days):%Y-%m-%d}/{today:%Y-%m-%d}",
                           query={"eo:cloud_cover": {"lt": max_cloud}}).item_collection()
    if len(found):
        items = list(found)
        break
if not items:
    raise RuntimeError("No sufficiently clear Sentinel-2 image in the last six months for this reservoir.")
latest = max(it.datetime for it in items).date()
day_items = [it for it in items if it.datetime.date() == latest]

img = load(day_items, bands=["red", "green", "blue", "nir08", "swir16", "scl"], bbox=grid["bbox"],
           crs=grid["crs"], resolution=grid["res"], groupby="solar_day", resampling="average").isel(time=0)
refl = {b: np.where(img[b].values > 0, img[b].values * 1e-4, np.nan)
        for b in ["red", "green", "blue", "nir08", "swir16"]}
clear = ~np.isin(img["scl"].values, [0, 1, 3, 8, 9, 10, 11])      # drop cloud, shadow, cirrus, snow
with np.errstate(invalid="ignore", divide="ignore"):
    mndwi = (refl["green"] - refl["swir16"]) / (refl["green"] + refl["swir16"])

# Historical footprint and analysis zone (same grid as the image)
occ = fetch_raster(f"{site}/jrc_occurrence.tif")[0][0]
zone = fetch_raster(f"{site}/water_frequency.tif")[0][0] != 255
footprint, core = (occ > 0) & (occ != 255), (occ >= 75) & (occ != 255)

water = clear & (mndwi > 0)
water |= clear & footprint & (refl["swir16"] < 0.03) & (refl["nir08"] < 0.05) & (refl["green"] < 0.10)  # very dark water
core_px = mndwi[clear & core]
if core_px.size > 100 and np.nanmedian(core_px) < 0.2:          # glint or haze on this image
    water |= clear & footprint & (refl["swir16"] < 0.12) & (refl["nir08"] < 0.15) & (mndwi > -0.3)
water &= zone
from scipy import ndimage                                      # keep only water connected to the reservoir
lab, n = ndimage.label(water, structure=np.ones((3, 3), bool))
if n:
    keep = np.unique(lab[footprint & (lab > 0)])
    water = np.isin(lab, keep[keep > 0])

area = water.sum() * (grid["res"] / 1000) ** 2
clear_share = (clear & zone).sum() / zone.sum()

rgb = np.dstack([refl["red"], refl["green"], refl["blue"]])
lo, hi = np.nanpercentile(rgb, 1), np.nanpercentile(rgb, 99.5)
rgb = np.nan_to_num(np.clip((rgb - lo) / (hi - lo), 0, 1))
fig, ax = plt.subplots(1, 2, figsize=(13, 5.2))
ax[0].imshow(rgb)
ax[0].set_title(f"Sentinel-2 true colour · {latest:%d %b %Y}")
ax[1].imshow(rgb.mean(axis=2), cmap="gray", vmin=0, vmax=1)
ax[1].imshow(np.ma.masked_where(~water, water), cmap=ListedColormap([COLORS["water_light"]]), alpha=0.85,
             interpolation="nearest")
ax[1].set_title(f"Open water detected · {area:,.1f} km²")
for a in ax:
    a.axis("off")
plt.tight_layout()
plt.show()

share = 100 * area / info["reservoir_max_extent_km2"]
print(f"{info['name']} on {latest:%d %B %Y}: {area:,.1f} km² of open water, {share:.0f} % of the largest "
      f"surface recorded in 1984–2021 ({info['reservoir_max_extent_km2']:,.1f} km²).")
if clear_share < 0.95:
    print(f"Clouds hid {100 * (1 - clear_share):.0f} % of the reservoir on this image, so this is a lower bound.")
print("Surface area is not volume: a 10 % drop in area can mean a much larger drop in storage in a steep valley.")

Beni Haroun Reservoir on 02 October 2026: 32.3 km² of open water, 87 % of the largest surface recorded in 1984–2021 (37.3 km²).
Surface area is not volume: a 10 % drop in area can mean a much larger drop in storage in a steep valley.

# @title Step 3 · Ten years of monthly change
monthly = fetch_csv(f"{site}/reservoir_monthly.csv")
monthly["date"] = pd.to_datetime(monthly["month"])
good = monthly[monthly["reported"] & ~monthly["flag_quality"]].copy()
setaside = monthly[monthly["reported"] & monthly["flag_quality"]]
trend = good.set_index("date")["water_area_km2"].rolling("365D", min_periods=6).mean()

fig = go.Figure()
fig.add_trace(go.Scatter(x=good["date"], y=good["water_area_km2"], mode="lines+markers", name="Monthly open-water area",
                         line=dict(color=COLORS["water_light"], width=1.5), marker=dict(size=5, color=COLORS["water"])))
fig.add_trace(go.Scatter(x=trend.index, y=trend.values, mode="lines", name="12-month average",
                         line=dict(color=COLORS["water"], width=3)))
if len(setaside):
    fig.add_trace(go.Scatter(x=setaside["date"], y=setaside["water_area_km2"], mode="markers",
                             name="Set aside (image quality)", marker=dict(symbol="x", size=8, color=COLORS["muted"])))
fig.add_hline(y=info["reservoir_max_extent_km2"], line_dash="dot", line_color=COLORS["muted"],
              annotation_text="largest surface 1984–2021", annotation_position="top left")
fig.update_layout(title=f"{info['name']} — open-water area from Sentinel-2", yaxis_title="km²", height=440, **LAYOUT)
fig.show()

hi_row = good.loc[good["water_area_km2"].idxmax()]
lo_row = good.loc[good["water_area_km2"].idxmin()]
first12, last12 = good.head(12)["water_area_km2"].mean(), good.tail(12)["water_area_km2"].mean()
print(f"Largest month: {hi_row['water_area_km2']:,.1f} km² ({hi_row['date']:%B %Y}) · "
      f"smallest: {lo_row['water_area_km2']:,.1f} km² ({lo_row['date']:%B %Y})")
print(f"First 12 months average {first12:,.1f} km² → latest 12 months average {last12:,.1f} km² "
      f"({100 * (last12 - first12) / first12:+.0f} %)")

freq, f_tr, _ = fetch_raster(f"{site}/water_frequency.tif")
occ = fetch_raster(f"{site}/jrc_occurrence.tif")[0][0]
freq = np.where(freq[0] == 255, np.nan, freq[0].astype(float))
occ = np.where(occ == 255, np.nan, occ.astype(float))
cmap = plt.get_cmap("Blues").copy()
cmap.set_bad("#F3F4F6")
fig, ax = plt.subplots(1, 2, figsize=(13, 5.4))
for a, arr, title in [(ax[0], freq, "Sentinel-2, 2017–2026\n(% of months)"),
                      (ax[1], occ, "Landsat, 1984–2021\n(% of observations)")]:
    im = a.imshow(np.where(arr > 0, arr, np.nan), cmap=cmap, vmin=0, vmax=100, interpolation="nearest")
    a.set_title(title)
    a.axis("off")
fig.suptitle("How often each pixel was covered by water", fontsize=13, fontweight="bold")
fig.colorbar(im, ax=ax, shrink=0.75, label="% of time water")
plt.show()
px = abs(f_tr.a) / 1000
print(f"Permanent water (≥ 90 % of months): {np.nansum(freq >= 90) * px**2:,.1f} km² · "
      f"zone that floods and dries (10–90 %): {np.nansum((freq >= 10) & (freq < 90)) * px**2:,.1f} km²")

Largest month: 38.8 km² (February 2019) · smallest: 29.1 km² (December 2022)
First 12 months average 34.9 km² → latest 12 months average 34.5 km² (-1 %)

Permanent water (≥ 90 % of months): 30.9 km² · zone that floods and dries (10–90 %): 5.9 km²

# @title Step 4 · Rainfall and drought upstream
rain = fetch_csv(f"{site}/rainfall_monthly.csv")
rain["date"] = pd.to_datetime(rain["month"])

fig = make_subplots(rows=2, cols=1, shared_xaxes=True, vertical_spacing=0.09, row_heights=[0.42, 0.58],
                    subplot_titles=("Monthly rainfall over the catchment (mm)",
                                    "12-month Standardized Precipitation Index (SPI-12)"))
fig.add_trace(go.Bar(x=rain["date"], y=rain["precip_mm"], name="Monthly rainfall", marker_color=COLORS["rain"],
                     opacity=0.8), row=1, col=1)
spi = rain["spi12"]
fig.add_trace(go.Bar(x=rain["date"], y=spi.where(spi >= 0), name="Wetter than normal", marker_color=COLORS["wet"]),
              row=2, col=1)
fig.add_trace(go.Bar(x=rain["date"], y=spi.where(spi < 0), name="Drier than normal", marker_color=COLORS["dry"]),
              row=2, col=1)
for level in (-1, -1.5, -2):
    fig.add_hline(y=level, line_dash="dot", line_color="#B45309", line_width=1, row=2, col=1)
fig.update_layout(title=f"{info['name']} catchment — rainfall and drought, 1981 to today", height=560, bargap=0,
                  **LAYOUT)
fig.show()

latest_spi = rain.dropna(subset=["spi12"]).iloc[-1]
state = ("extreme drought" if latest_spi.spi12 <= -2 else "severe drought" if latest_spi.spi12 <= -1.5 else
         "moderate drought" if latest_spi.spi12 <= -1 else "near normal" if latest_spi.spi12 < 1 else "wet")
print(f"Latest SPI-12 ({pd.Timestamp(latest_spi.month):%B %Y}): {latest_spi.spi12:+.2f} — {state}.")

# Does the reservoir follow the rain? Correlate monthly area with SPI-12 shifted by 0–12 months.
area_s = good.set_index("date")["water_area_km2"]
spi_s = rain.set_index("date")["spi12"]
scores = []
for lag in range(13):
    pair = pd.concat([area_s, spi_s.shift(lag)], axis=1, join="inner").dropna()
    if len(pair) >= 24:
        scores.append((lag, pair.corr().iloc[0, 1], len(pair)))
lag, r, n = max(scores, key=lambda s: s[1])

fig = make_subplots(specs=[[{"secondary_y": True}]])
fig.add_trace(go.Scatter(x=area_s.index, y=area_s.values, name="Open-water area (km²)",
                         line=dict(color=COLORS["water"], width=2.5)), secondary_y=False)
shifted = spi_s.shift(lag).loc[area_s.index.min():]
fig.add_trace(go.Scatter(x=shifted.index, y=shifted.values, name=f"SPI-12, {lag} month{'' if lag == 1 else 's'} earlier",
                         line=dict(color=COLORS["dry"], width=1.8, dash="dot")), secondary_y=True)
fig.update_yaxes(title_text="km²", secondary_y=False)
fig.update_yaxes(title_text="SPI-12", secondary_y=True)
fig.update_layout(title=f"Reservoir area and upstream drought index (best match: {lag} month lag, r = {r:.2f})",
                  height=420, **LAYOUT)
fig.show()
print(f"Best match: the reservoir follows the 12-month rainfall index {lag} month{'' if lag == 1 else 's'} later "
      f"(Pearson r = {r:.2f}, {n} months).")
print("A strong link means storage is driven by rainfall in the catchment; a weak one points to dam operations, "
      "irrigation withdrawals, inflow from far upstream, or groundwater.")

Latest SPI-12 (September 2026): +0.02 — near normal.

Best match: the reservoir follows the 12-month rainfall index 0 months later (Pearson r = 0.53, 115 months).
A strong link means storage is driven by rainfall in the catchment; a weak one points to dam operations, irrigation withdrawals, inflow from far upstream, or groundwater.

# @title Step 5 · Water-quality signals
q = monthly[monthly["reported"] & monthly["ndci_median"].notna()]
fig = make_subplots(rows=2, cols=1, shared_xaxes=True, vertical_spacing=0.1,
                    subplot_titles=("Chlorophyll indicator — NDCI (higher = greener water)",
                                    "Turbidity indicator — NDTI (higher = more sediment)"))
fig.add_trace(go.Scatter(x=q["date"], y=q["ndci_median"], mode="lines+markers", name="NDCI",
                         line=dict(color=COLORS["accent"])), row=1, col=1)
fig.add_trace(go.Scatter(x=q["date"], y=q["ndti_median"], mode="lines+markers", name="NDTI",
                         line=dict(color="#8D6E63")), row=2, col=1)
fig.update_layout(title=f"{info['name']} — monthly median over open water", height=520, **LAYOUT)
fig.show()

season = q.assign(m=q["date"].dt.month).groupby("m")["ndci_median"].mean()
print(f"On average the water is greenest in {pd.Timestamp(2000, int(season.idxmax()), 1):%B} and clearest "
      f"in {pd.Timestamp(2000, int(season.idxmin()), 1):%B} (NDCI, {len(q)} months with glint-free images).")

qa = fetch_raster(f"{site}/latest_quality.tif")[0]
rgb_img = np.moveaxis(fetch_raster(f"{site}/latest_rgb.tif")[0], 0, -1)
ndci = np.where(qa[0] == 255, np.nan, qa[0] / 127 - 1)
ndti = np.where(qa[1] == 255, np.nan, qa[1] / 127 - 1)
freq_q = fetch_raster(f"{site}/water_frequency.tif")[0][0]
outline = np.where((freq_q != 255) & (freq_q >= 10), 1.0, np.nan)      # reservoir footprint, for context
fig, ax = plt.subplots(1, 3, figsize=(16, 5))
ax[0].imshow(rgb_img)
ax[0].set_title("True colour · " + " and ".join(info.get("latest_image_dates") or [info["latest_image_date"]]))
for a in ax[1:]:
    a.imshow(outline, cmap=ListedColormap(["#E5E7EB"]), interpolation="nearest")
im1 = ax[1].imshow(ndci, cmap="YlGn", vmin=-0.2, vmax=0.4, interpolation="nearest")
ax[1].set_title("NDCI — chlorophyll indicator")
im2 = ax[2].imshow(ndti, cmap="YlOrBr", vmin=-0.4, vmax=0.2, interpolation="nearest")
ax[2].set_title("NDTI — turbidity indicator")
for a in ax:
    a.axis("off")
fig.colorbar(im1, ax=ax[1], shrink=0.7)
fig.colorbar(im2, ax=ax[2], shrink=0.7)
plt.tight_layout()
plt.show()

On average the water is greenest in May and clearest in January (NDCI, 112 months with glint-free images).

# @title Step 6 · AI change radar, 2018 → 2025
from scipy import ndimage
from pyproj import Transformer

try:
    arr, c_tr, c_crs = fetch_raster(f"{site}/ai_change.tif")
except Exception:
    arr = None
if arr is None:
    print("The AI change layer is not available for this reservoir.")
else:
    change = np.where(arr[0] == 255, np.nan, arr[0] / 250.0)      # 0 = unchanged … 1 = completely different
    thr = max(0.3, float(np.nanpercentile(change, 98)))
    patches, n = ndimage.label(np.nan_to_num(change) >= thr)
    sizes = ndimage.sum(np.ones_like(change), patches, index=np.arange(1, n + 1)) if n else np.array([])
    top = (np.argsort(sizes)[::-1][:5] + 1) if n else []
    to_ll = Transformer.from_crs(c_crs, "EPSG:4326", always_xy=True)
    res_m = abs(c_tr.a)

    fig, ax = plt.subplots(figsize=(10, 6))
    cm = plt.get_cmap("magma").copy()
    cm.set_bad("#F3F4F6")
    im = ax.imshow(change, cmap=cm, vmin=0, vmax=0.6, interpolation="nearest")
    for k, lab in enumerate(top, start=1):
        r0, c0 = ndimage.center_of_mass(patches == lab)
        ax.annotate(str(k), (c0, r0), color="white", fontsize=12, fontweight="bold", ha="center", va="center",
                    bbox=dict(boxstyle="circle", fc=COLORS["accent"], ec="white"))
    ax.set_title(f"How much each place changed, 2018 → 2025 ({res_m:.0f} m pixels)")
    ax.axis("off")
    fig.colorbar(im, ax=ax, shrink=0.75, label="embedding change (1 − cosine similarity)")
    plt.show()

    # Interactive map: the change layer over satellite imagery, with the five largest hot spots
    dst_crs = "EPSG:4326"
    west, south, east, north = transform_bounds(c_crs, dst_crs, *rasterio.transform.array_bounds(*change.shape, c_tr))
    from rasterio.transform import from_bounds
    h, w = change.shape
    ll = np.full((h, w), np.nan, np.float32)
    reproject(change.astype(np.float32), ll, src_transform=c_tr, src_crs=c_crs, src_nodata=np.nan,
              dst_transform=from_bounds(west, south, east, north, w, h), dst_crs=dst_crs, dst_nodata=np.nan,
              resampling=Resampling.nearest)
    rgba = cm(np.clip(np.nan_to_num(ll) / 0.6, 0, 1))
    rgba[..., 3] = np.where(np.isfinite(ll) & (ll > 0.15), 0.75, 0)
    m = folium.Map(tiles=None, control_scale=True, height=520)
    folium.TileLayer(tiles="https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{z}/{y}/{x}",
                     attr="Imagery © Esri, Maxar, Earthstar Geographics", name="Satellite imagery").add_to(m)
    folium.raster_layers.ImageOverlay(rgba, bounds=[[south, west], [north, east]], name="AI change 2018 → 2025",
                                      mercator_project=True).add_to(m)
    for k, lab in enumerate(top, start=1):
        r0, c0 = ndimage.center_of_mass(patches == lab)
        x, y = c_tr * (c0 + 0.5, r0 + 0.5)
        lon, lat = to_ll.transform(x, y)
        area_ha = sizes[lab - 1] * res_m ** 2 / 1e4
        folium.Marker([lat, lon], tooltip=f"Hot spot {k}: about {area_ha:,.0f} ha changed",
                      icon=folium.DivIcon(html=f"<div style='background:{COLORS['accent']};color:white;border-radius:50%;"
                                               f"width:24px;height:24px;text-align:center;font-weight:bold;"
                                               f"line-height:24px;border:2px solid white'>{k}</div>")).add_to(m)
    folium.LayerControl(collapsed=False).add_to(m)
    m.fit_bounds([[south, west], [north, east]])
    display(m)
    print("Zoom into each numbered hot spot and switch the change layer on and off: what happened there?")

Zoom into each numbered hot spot and switch the change layer on and off: what happened there?

Step	What you will do	Data
1	Choose a reservoir and see its catchment	JRC Global Surface Water · HydroBASINS
2	Map today's open water from the live satellite archive	Sentinel-2 (ESA Copernicus)
3	Follow ten years of monthly change	Sentinel-2, 2017 → today
4	Link the reservoir to rainfall and drought upstream	CHIRPS rainfall, 1981 → today
5	Read water-quality signals from space	Sentinel-2 red-edge bands
6	Let an AI model flag what changed around the reservoir	AlphaEarth Foundations embeddings

Watching reservoirs from space¶

Open satellite data and AI for water security · IC-TWEP 2026 workshop¶

Step 1 — Choose a reservoir¶

Beni Haroun Reservoir

Step 2 — Map today's water from the open Sentinel-2 archive¶

Step 3 — Ten years of change, month by month¶

Step 4 — Rainfall and drought upstream¶

Step 5 — Water-quality signals from space¶

Step 6 — AI change radar: what changed around the reservoir?¶

What you just did¶

Data sources¶

Country	Algeria
Largest water surface seen 1984–2021	37.3 km²
Upstream catchment	7,574 km²
Mean rainfall over the catchment, 1981–2026	556 mm per year
Sentinel-2 monthly record	116 months, 2017 to 2026