Your own reservoir monitor¶

Take-home notebook · Workshop 1 · IC-TWEP 2026 · QuestFeed¶

Give the notebook the coordinates of any reservoir on Earth. It will:

  1. find the reservoir in the 1984–2021 JRC Global Surface Water record;
  2. process every Sentinel-2 image of the last months with exactly the method used for the workshop's ten-year records (shown in full below, so you can read and adapt it);
  3. produce a monthly open-water series, a map and files you can download for a report.

Everything is read from open archives; no account other than Google (for Colab) is needed. A 12-month run on a medium-sized reservoir takes about 5–10 minutes; longer periods take proportionally longer.

Ready.
In [2]:
# @title 1 · Your reservoir
NAME = "Beni Haroun Reservoir"  # @param {type:"string"}
LAT = 36.555  # @param {type:"number"}
LON = 6.270  # @param {type:"number"}
MONTHS = 12  # @param {type:"slider", min:6, max:60, step:6}
RES = 20  # @param [10, 20, 40, 80] {type:"raw"}

pcat = Client.open("https://planetarycomputer.microsoft.com/api/stac/v1", modifier=pc.sign_inplace)
half = 0.3
box_ = [LON - half, LAT - half, LON + half, LAT + half]
jrc = load(list(pcat.search(collections=["jrc-gsw"], bbox=box_).items()), bands=["occurrence"], bbox=box_,
           crs="EPSG:4326", resolution=0.00025, resampling="nearest")["occurrence"].isel(time=0)
occ = jrc.values
lab, n = ndimage.label((occ >= 10) & (occ <= 100), structure=np.ones((3, 3)))
iy, ix = np.abs(jrc.latitude.values - LAT).argmin(), np.abs(jrc.longitude.values - LON).argmin()
yy, xx = np.nonzero(lab)
k = lab[yy, xx][np.argmin((yy - iy) ** 2 + (xx - ix) ** 2)]
body = lab == k
rows, cols = np.nonzero(body)
lat_v, lon_v = jrc.latitude.values, jrc.longitude.values
pad = 0.01
BBOX = [lon_v[cols].min() - pad, lat_v[rows].min() - pad, lon_v[cols].max() + pad, lat_v[rows].max() + pad]
max_km2 = body.sum() * (0.00025 * 111.32) ** 2 * np.cos(np.radians(LAT))
print(f"{NAME}: water body found {np.hypot((yy - iy), (xx - ix)).min() * 0.00025 * 111:.1f} km from your point; "
      f"largest surface 1984–2021 ≈ {max_km2:,.1f} km²")
CRS = f"EPSG:{(32600 if LAT >= 0 else 32700) + int((LON + 180) // 6) + 1}"

Beni Haroun Reservoir: water body found 0.0 km from your point; largest surface 1984–2021 ≈ 37.3 km²

The method (identical to the workshop's ten-year records)¶

  • clear pixels from Sentinel-2's scene classification (no cloud, shadow, cirrus or snow);
  • water = MNDWI > 0; plus very dark water (SWIR1 < 0.03, NIR < 0.05, green < 0.10) inside the reservoir's historical footprint; plus, on images degraded by glint or haze, pixels in the footprint that stay dark in the near- and short-wave infrared;
  • only water connected to the reservoir is kept;
  • images that see much less water than the month's best clear image are set aside (an image that sees far more water than the reservoir has ever held is not used as that reference); a pixel is water in a month if it was water on most remaining clear images; pixels never seen in a month are taken from the month before or after; months need ≥ 95 % of the reservoir observed;
  • a one-month dip that recovers at once is flagged, because it is almost always an image problem, not the reservoir.
In [3]:
# @title 2 · Run the monitor (about 5–10 minutes for 12 months)
footprint_occ = load(list(pcat.search(collections=["jrc-gsw"], bbox=BBOX).items()), bands=["occurrence"], bbox=BBOX,
                     crs=CRS, resolution=RES, resampling="nearest")["occurrence"].isel(time=0).values
lab2, _ = ndimage.label((footprint_occ >= 10) & (footprint_occ <= 100), structure=np.ones((3, 3)))
big = np.bincount(lab2.ravel())[1:].argmax() + 1 if lab2.max() else 0
comp = lab2 == big
footprint = comp | (ndimage.binary_dilation(comp, iterations=max(1, int(1000 / RES))) & (footprint_occ > 0) & (footprint_occ <= 100))
core = footprint & (footprint_occ >= 75) & (footprint_occ <= 100)
zone = ndimage.binary_dilation(footprint, iterations=max(1, int(300 / RES)))
EXCL = [0, 1, 3, 8, 9, 10, 11]          # no data, saturated, cloud shadow, cloud, cirrus, snow
PX_KM2 = (RES / 1000) ** 2

es = Client.open("https://earth-search.aws.element84.com/v1")
end = pd.Timestamp.now(tz="UTC").normalize()
start = (end - pd.DateOffset(months=MONTHS)).replace(day=1)
found = list(es.search(collections=["sentinel-2-l2a"], bbox=BBOX, datetime=f"{start:%Y-%m-%d}/{end:%Y-%m-%d}",
                       query={"eo:cloud_cover": {"lt": 70}}).items())
best_item = {}                          # one item per tile and acquisition: the most recently processed
for it in found:
    key = (it.properties.get("grid:code"), it.datetime.strftime("%Y%m%dT%H%M"))
    rank = it.properties.get("updated") or it.properties.get("created") or ""
    if key not in best_item or rank > best_item[key][0]:
        best_item[key] = (rank, it)
items = [v[1] for v in best_item.values()]
# where tiles overlap, every date is delivered twice: keep one tile if it covers the whole area on its own
from shapely.geometry import box as _box, shape as _shape
aoi = _box(*BBOX)
by_tile = {}
for it in items:
    by_tile.setdefault(it.properties.get("grid:code"), []).append(it)
def _epsg(it):
    return it.properties.get("proj:epsg") or int(str(it.properties.get("proj:code", "0")).split(":")[-1] or 0)
def _rank(tile):
    its_ = by_tile[tile]
    cover = np.median([_shape(i.geometry).intersection(aoi).area / aoi.area for i in its_])
    return (cover >= 0.99, f"EPSG:{_epsg(its_[0])}" == CRS, len(its_))
best_tile = max(by_tile, key=_rank)
if len(by_tile) > 1 and _rank(best_tile)[0]:
    items = by_tile[best_tile]
print(f"{len(items)} Sentinel-2 images to process")

months = [str(m) for m in pd.period_range(start, end, freq="M")]
freqs, used, last_water, last_day, t0 = {}, {}, None, None, time.time()
for mon in months:
    its = [i for i in items if i.datetime.strftime("%Y-%m") == mon]
    if not its:
        continue
    ds = load(its, bands=["green", "nir08", "swir16", "scl"], bbox=BBOX, crs=CRS, resolution=RES, groupby="solar_day",
              resampling={"*": "average", "scl": "nearest"}, pool=16, fail_on_error=False)
    dates = []
    for t in range(ds.sizes["time"]):
        g, nir, sw = (np.where(ds[b].values[t] > 0, ds[b].values[t] * 1e-4, np.nan) for b in ("green", "nir08", "swir16"))
        valid = ~np.isin(ds["scl"].values[t], EXCL) & np.isfinite(g) & np.isfinite(sw)
        with np.errstate(invalid="ignore", divide="ignore"):
            m = (g - sw) / (g + sw)
        cpx = m[valid & core]
        if cpx.size < 100:
            cpx = m[valid & footprint]
        water = (valid & (m > 0)) | (valid & footprint & (sw < 0.03) & (nir < 0.05) & (g < 0.10))
        if cpx.size and np.nanmedian(cpx) < 0.2:                       # glint or haze on this image
            water |= valid & footprint & (sw < 0.12) & (nir < 0.15) & (m > -0.3)
        lab3, n3 = ndimage.label(water & zone, structure=np.ones((3, 3), bool))
        if n3:
            keep = np.unique(lab3[footprint & (lab3 > 0)]); water = np.isin(lab3, keep[keep > 0])
        vz = valid & zone
        if vz.sum():
            dates.append((str(ds.time.values[t])[:10], vz, water & zone, (water & zone).sum() / vz.sum(), vz.sum() / zone.sum(),
                          (water & zone).sum() * PX_KM2))
    cand = [(f, w) for (_, _, _, f, c, w) in dates if c >= 0.5]
    sane = [f for f, w in cand if w <= 1.2 * max_km2]                      # far above the largest surface = artefact
    ref = max(sane or [f for f, _ in cand], default=None)
    nv = np.zeros(zone.shape, np.uint16); nw = np.zeros(zone.shape, np.uint16); n_used = 0
    for day, vz, wz, f, c, _w in dates:
        if ref and f < 0.6 * ref:                                      # sees far less water than the best image
            continue
        nv += vz; nw += wz; n_used += 1
        if c > 0.98:
            last_water, last_day = wz, day
    freqs[mon] = np.where(nv > 0, nw / np.maximum(nv, 1), np.nan)
    used[mon] = (len(dates), n_used)
    print(f"{mon}: {len(dates)} images, {n_used} used  ({time.time() - t0:.0f} s)")

records = []
for k, mon in enumerate(months):
    if mon not in freqs:
        records.append(dict(month=mon, images=0, images_used=0, observed_share=0.0, water_km2=np.nan)); continue
    filled = freqs[mon].copy()
    for j in (k - 1, k + 1):                                           # gaps filled from the month before / after
        if 0 <= j < len(months) and months[j] in freqs:
            gap = ~np.isfinite(filled) & np.isfinite(freqs[months[j]])
            filled[gap] = freqs[months[j]][gap]
    obs = np.isfinite(filled) & zone
    cover = obs.sum() / zone.sum()
    area = ((filled >= 0.5) & obs).sum() * PX_KM2
    records.append(dict(month=mon, images=used[mon][0], images_used=used[mon][1], observed_share=round(cover, 3),
                        water_km2=round(area, 3) if cover >= 0.95 else np.nan))
series = pd.DataFrame(records)
v = series["water_km2"].to_numpy(float)
flag = np.zeros(len(v), bool)
for i in np.flatnonzero(np.isfinite(v)):                               # one-month dips that recover at once
    nb = np.r_[v[max(0, i - 2):i], v[i + 1:i + 3]]; nb = nb[np.isfinite(nb)]
    nxt = v[i + 1:i + 3]; nxt = nxt[np.isfinite(nxt)]
    flag[i] = len(nb) >= 2 and len(nxt) >= 1 and v[i] < 0.75 * np.median(nb) and v[i] < 0.8 * nxt[0]
series["flagged_dip"] = flag
print(f"\nDone in {time.time() - t0:.0f} s: {int(np.isfinite(v).sum())} of {len(months)} months reported, "
      f"{int(flag.sum())} flagged.")

81 Sentinel-2 images to process

2025-10: 7 images, 6 used  (57 s)

2025-11: 7 images, 7 used  (79 s)

2025-12: 8 images, 8 used  (105 s)

2026-01: 4 images, 3 used  (137 s)

2026-02: 2 images, 2 used  (159 s)

2026-03: 3 images, 3 used  (193 s)

2026-04: 5 images, 3 used  (220 s)

2026-05: 7 images, 7 used  (262 s)

2026-06: 9 images, 8 used  (298 s)

2026-07: 10 images, 10 used  (332 s)

2026-08: 9 images, 9 used  (375 s)

2026-09: 7 images, 7 used  (400 s)

2026-10: 3 images, 3 used  (422 s)

Done in 422 s: 13 of 13 months reported, 1 flagged.
In [4]:
# @title 3 · Results and downloadable report files
series["date"] = pd.to_datetime(series["month"])
ok = series[series["water_km2"].notna() & ~series["flagged_dip"]]
bad = series[series["flagged_dip"]]
fig = go.Figure(go.Scatter(x=ok["date"], y=ok["water_km2"], mode="lines+markers", name="open water",
                           line=dict(color=COLORS["water"], width=3)))
if len(bad):
    fig.add_trace(go.Scatter(x=bad["date"], y=bad["water_km2"], mode="markers", name="flagged dip",
                             marker=dict(symbol="x", size=9, color=COLORS["muted"])))
fig.add_hline(y=max_km2, line_dash="dot", line_color=COLORS["muted"], annotation_text="largest surface 1984–2021")
fig.update_layout(title=f"{NAME} — monthly open-water area (Sentinel-2)", yaxis_title="km²", height=420, **LAYOUT)
fig.show()
series.drop(columns="date").to_csv("reservoir_monitor.csv", index=False)
if last_water is not None:
    fig2, ax = plt.subplots(figsize=(8, 6))
    ax.imshow(np.where(footprint, 0.15, np.nan), cmap="Greys", vmin=0, vmax=1)
    ax.imshow(np.ma.masked_where(~last_water, last_water), cmap=ListedColormap([COLORS["water_light"]]))
    ax.set_title(f"{NAME}: open water on {last_day} (grey = 1984–2021 footprint)"); ax.axis("off")
    fig2.savefig("reservoir_monitor_map.png", dpi=150, bbox_inches="tight"); plt.show()
print("Files written: reservoir_monitor.csv" + (", reservoir_monitor_map.png" if last_water is not None else ""))
try:
    from google.colab import files  # in Colab, offer the files for download
    files.download("reservoir_monitor.csv")
except Exception:
    pass
No description has been provided for this image 
Files written: reservoir_monitor.csv, reservoir_monitor_map.png

Make it operational¶

  • Run this notebook every month, or run the same code on a schedule on any server.
  • Add rainfall and the drought index for the catchment (see the live notebook, Step 4) to the same report.
  • Alert when the area falls below a threshold you choose — for example the lowest month of the past ten years.

Data: Sentinel-2 L2A (contains modified Copernicus Sentinel data, ESA; Earth Search on AWS) · JRC Global Surface Water (Pekel et al. 2016; Microsoft Planetary Computer).

QuestFeed · IC-TWEP 2026 take-home materials · questfeed.com