Showcase — the 2022 Pakistan floods (Sentinel-1 SAR)¶
The 2022 monsoon submerged roughly a third of Pakistan. Optical sensors are blocked by monsoon cloud, so Sentinel-1 SAR is the standard tool: smooth open water reflects radar away and reads very dark. This notebook renders the SAR composites with pyramids' Dataset.plot(), picks a water threshold automatically with Otsu's method, and maps the newly-flooded land.
What this notebook does¶
- Fetch a pre-monsoon and a peak-flood Sentinel-1 VH composite over Sindh.
- Render both with
Dataset.plot(band=0, cmap='gray'). - Otsu threshold the bimodal SAR histogram, then map and measure the new water.
Needs
GEE_SERVICE_ACCOUNT/GEE_SERVICE_KEYandpyramids-gis[viz].
Setup¶
First the imports and a bit of noise-suppression. pyramids provides Dataset (reading + plotting); earthlens provides the unified EarthLens entry point.
import os
import tempfile
import warnings
import matplotlib.pyplot as plt
import numpy as np
from dotenv import find_dotenv, load_dotenv
from loguru import logger
from pyramids.dataset import Dataset # pyramids' own reading + plotting
from earthlens import EarthLens
warnings.filterwarnings("ignore")
logger.remove()
2026-06-13 18:15:10 | INFO | pyramids.base.config | Logging is configured.
Credentials¶
The notebook reads its credentials from a .env file at the repository root (loaded with python-dotenv; your PATH is left untouched). Create that .env with at least:
GEE_SERVICE_ACCOUNT="your-sa@your-project.iam.gserviceaccount.com"
GEE_SERVICE_KEY="/path/to/your/service-account-key.json"
load_dotenv(find_dotenv(usecwd=True))
SERVICE_ACCOUNT = os.environ["GEE_SERVICE_ACCOUNT"] # Earth Engine service-account email (.env)
SERVICE_KEY = os.environ["GEE_SERVICE_KEY"] # path to the service-account JSON key (.env)
1 · Fetch the Sentinel-1 VH composites¶
Sentinel-1 VH backscatter (dB) is a single radar channel shown in grey: open water is specular and so appears dark — that is what makes flooding visible. We define the AOI over Sindh and the shared request parameters once, then fetch a pre-monsoon (before) and a peak-flood (after) median composite. Each fetch is split into construct → authenticate → load so the steps stay easy to read.
AOI = [67.8, 26.6, 68.8, 27.6]
C = "COPERNICUS/S1_GRD"
The pre-monsoon baseline is a mean VH composite over July 2022, before the flood peak.
before_scene = EarthLens(
data_source="gee",
cadence="raw",
dataset=C,
variables=["VH"],
aoi=AOI,
start="2022-07-01",
end="2022-07-31",
scale=60,
reducer="mean",
)
before_scene.authenticate(service_account=SERVICE_ACCOUNT, service_key=SERVICE_KEY)
before_rasters = before_scene.load()
before = before_rasters[0]
COPERNICUS/S1_GRD [VH]: 0%| | 0/1 [00:00<?, ?img/s]
COPERNICUS/S1_GRD [VH]: 100%|██████████| 1/1 [00:08<00:00, 8.73s/img]
COPERNICUS/S1_GRD [VH]: 100%|██████████| 1/1 [00:08<00:00, 8.73s/img]
The peak-flood composite covers late August to mid-September 2022, when inundation was at its worst.
after_scene = EarthLens(
data_source="gee",
cadence="raw",
dataset=C,
variables=["VH"],
aoi=AOI,
start="2022-08-25",
end="2022-09-20",
scale=60,
reducer="mean",
)
after_scene.authenticate(service_account=SERVICE_ACCOUNT, service_key=SERVICE_KEY)
after_rasters = after_scene.load()
after = after_rasters[0]
print("composites:", bool(before), bool(after))
COPERNICUS/S1_GRD [VH]: 0%| | 0/1 [00:00<?, ?img/s]
COPERNICUS/S1_GRD [VH]: 100%|██████████| 1/1 [00:09<00:00, 9.66s/img]
COPERNICUS/S1_GRD [VH]: 100%|██████████| 1/1 [00:09<00:00, 9.66s/img]
composites: True True
2 · Before and peak flood — Dataset.plot()¶
pyramids renders each VH composite (dB) side by side; the dark fabric is open water. The before/after pair makes the new inundation jump out.
if before and after:
fig, ax = plt.subplots(1, 2, figsize=(12, 6))
before.plot(band=0, ax=ax[0], cmap="gray", title="VH — Jul 2022 (before)")
after.plot(band=0, ax=ax[1], cmap="gray", title="VH — Aug/Sep 2022 (flood)")
plt.show()
3 · Otsu threshold on the SAR histogram¶
The SAR histogram is bimodal — a dark water mode and a brighter land mode. otsu() is a dependency-free Otsu implementation that returns the dB value maximising between-class variance, i.e. the split between the two modes.
def otsu(values, bins=256):
v = values[np.isfinite(values)]
hist, edges = np.histogram(v, bins=bins)
centres = (edges[:-1] + edges[1:]) / 2
w = np.cumsum(hist)
wb = w / w[-1]
mb = np.cumsum(hist * centres) / np.maximum(w, 1)
mt = mb[-1]
between = (
wb
* (1 - wb)
* (mb - np.where(wb < 1, (mt - mb * wb) / np.maximum(1 - wb, 1e-9), mb)) ** 2
)
return centres[np.nanargmax(between)]
Apply the threshold to the peak-flood composite and plot the histogram: water falls left of the line, land to the right.
if before and after:
a = np.asarray(after.read_array(), dtype="float32")
a = a if a.ndim == 2 else a[0]
thr = otsu(a.ravel())
fig, ax = plt.subplots(figsize=(8, 3.4))
ax.hist(a[np.isfinite(a)].ravel(), bins=120, color="0.6")
ax.axvline(thr, color="tab:blue", lw=2, label=f"Otsu threshold = {thr:.1f} dB")
ax.set(
xlabel="VH backscatter (dB)",
ylabel="pixels",
title="SAR histogram: water (left) vs land (right)",
)
ax.legend()
plt.show()
4 · Flood map and inundated area¶
New water is everything below the Otsu threshold now that was above it before — i.e. land that turned to water. The flood overlay needs pixel-aligned masking, so this map is drawn with matplotlib (the standalone SAR stills above use Dataset.plot).
if before and after:
b = np.asarray(before.read_array(), dtype="float32")
b = b if b.ndim == 2 else b[0]
new_water = (a < thr) & ~(b < thr)
px_km2 = (60 * 60) / 1e6
def db(x):
lo, hi = np.nanpercentile(x, 2), np.nanpercentile(x, 98)
return np.clip((x - lo) / (hi - lo), 0, 1)
fig, ax = plt.subplots(figsize=(7, 7))
ax.imshow(db(a), cmap="gray")
ov = np.zeros((*new_water.shape, 4))
ov[new_water] = [1, 0, 0, 0.55]
ax.imshow(ov)
ax.set_title("newly-flooded land (red)")
ax.axis("off")
plt.show()
print(f"newly-flooded area in this tile: {new_water.sum() * px_km2:,.0f} km^2")
newly-flooded area in this tile: 2,388 km^2
5 · Flood extent over time — the recession¶
A single before/after pair shows the peak, but the flood is a process: it surged in late August and drained slowly over the following months. Re-applying the Otsu water test to a series of Sentinel-1 composites traces the flooded area from the pre-monsoon baseline through the peak and the long recession — the monitoring view emergency agencies track week to week.
First fix the permanent-water baseline (rivers and lakes that are wet even in the dry pre-monsoon scene), and list the monitoring windows to sample.
if before:
b = np.asarray(before.read_array(), dtype="float32")
b = b if b.ndim == 2 else b[0]
base_water = b < otsu(b.ravel()) # permanent rivers/lakes
windows = [
("Jul (pre)", "2022-07-01", "2022-07-31"),
("late Aug", "2022-08-25", "2022-09-05"),
("mid Sep", "2022-09-10", "2022-09-25"),
("early Oct", "2022-10-01", "2022-10-18"),
("Nov", "2022-11-05", "2022-11-25"),
]
px_km2 = (60 * 60) / 1e6
For each window, fetch the VH composite, run the Otsu water test, and record the newly-flooded area relative to the permanent-water baseline. Each fetch follows the same construct → authenticate → load form.
if before:
labels, areas = [], []
for lab, s, e in windows:
window_scene = EarthLens(
data_source="gee",
cadence="raw",
dataset=C,
variables=["VH"],
aoi=AOI,
start=s,
end=e,
scale=60,
reducer="mean",
)
window_scene.authenticate(service_account=SERVICE_ACCOUNT, service_key=SERVICE_KEY)
window_rasters = window_scene.load()
gp = window_rasters[0]
if not gp:
continue
a = np.asarray(gp.read_array(), dtype="float32")
a = a if a.ndim == 2 else a[0]
new_water = (a < otsu(a.ravel())) & ~base_water
labels.append(lab)
areas.append(float(new_water.sum() * px_km2))
COPERNICUS/S1_GRD [VH]: 0%| | 0/1 [00:00<?, ?img/s]
COPERNICUS/S1_GRD [VH]: 100%|██████████| 1/1 [00:12<00:00, 12.33s/img]
COPERNICUS/S1_GRD [VH]: 100%|██████████| 1/1 [00:12<00:00, 12.33s/img]
COPERNICUS/S1_GRD [VH]: 0%| | 0/1 [00:00<?, ?img/s]
COPERNICUS/S1_GRD [VH]: 100%|██████████| 1/1 [00:06<00:00, 6.15s/img]
COPERNICUS/S1_GRD [VH]: 100%|██████████| 1/1 [00:06<00:00, 6.15s/img]
COPERNICUS/S1_GRD [VH]: 0%| | 0/1 [00:00<?, ?img/s]
COPERNICUS/S1_GRD [VH]: 100%|██████████| 1/1 [00:08<00:00, 8.25s/img]
COPERNICUS/S1_GRD [VH]: 100%|██████████| 1/1 [00:08<00:00, 8.25s/img]
COPERNICUS/S1_GRD [VH]: 0%| | 0/1 [00:00<?, ?img/s]
COPERNICUS/S1_GRD [VH]: 100%|██████████| 1/1 [00:07<00:00, 7.87s/img]
COPERNICUS/S1_GRD [VH]: 100%|██████████| 1/1 [00:07<00:00, 7.87s/img]
COPERNICUS/S1_GRD [VH]: 0%| | 0/1 [00:00<?, ?img/s]
COPERNICUS/S1_GRD [VH]: 100%|██████████| 1/1 [00:07<00:00, 7.74s/img]
COPERNICUS/S1_GRD [VH]: 100%|██████████| 1/1 [00:07<00:00, 7.74s/img]
Plot the flooded-area time series: the surge to the peak and the slow recession that follows.
if before and areas:
fig, ax = plt.subplots(figsize=(8, 4))
ax.plot(range(len(areas)), areas, "-o", color="tab:blue")
ax.fill_between(range(len(areas)), 0, areas, alpha=0.2)
ax.set_xticks(range(len(labels)))
ax.set_xticklabels(labels)
ax.set(
ylabel="newly-flooded area (km²)",
title="Flood extent over time — peak and recession (Sindh tile)",
)
plt.show()
print(
"flooded area (km²) by window:",
dict(zip(labels, [round(v) for v in areas])),
)
elif not before:
print("no baseline — set GEE credentials and rerun")
flooded area (km²) by window: {'Jul (pre)': 0, 'late Aug': 2108, 'mid Sep': 1502, 'early Oct': 1125, 'Nov': 1109}
Recap¶
The SAR stills are rendered with Dataset.plot(); the flood-overlay map stays in matplotlib because it composites a pixel-aligned mask.
Try it yourself¶
- Slide the
AOIalong the Indus, or compare VH vs VV. - See the GEE backend reference.