Notebook 04 — Change Detection¶

Computes pixel-wise before/after change maps for both NDWI and SAM water masks.

Change classes (encoding: before × 2 + after):

class	name	meaning
0	Stable Land	non-water in both scenes
1	Water Gain	new water (flooding / inundation)
2	Water Loss	drained / dried
3	Stable Water	permanent water bodies

Inputs : masks/
Outputs : change/change_ndwi.nc, change/change_sam.nc, updated scene_metadata.json

1. Imports & Config¶

In [1]:

Copied!





import json, tomllib, warnings, math
import numpy as np
import pandas as pd
import xarray as xr
import rioxarray
import hvplot.xarray
import holoviews as hv
from pathlib import Path

warnings.filterwarnings("ignore")
hv.extension("bokeh")
import json, tomllib, warnings, math
import numpy as np
import pandas as pd
import xarray as xr
import rioxarray
import hvplot.xarray
import holoviews as hv
from pathlib import Path

warnings.filterwarnings("ignore")
hv.extension("bokeh")

In [2]:

Copied!





with open(Path("config.toml"), "rb") as f:
    cfg = tomllib.load(f)

DATA_ROOT  = Path(cfg["data"]["root"])
MASKS_DIR  = DATA_ROOT / "masks"
CHANGE_DIR = DATA_ROOT / "change"
META_FILE  = DATA_ROOT / "scene_metadata.json"
CHANGE_DIR.mkdir(parents=True, exist_ok=True)

PIXEL_RES  = cfg["processing"]["pixel_res_m"]
CROP_SIZE  = cfg["processing"]["crop_size"]
PIXEL_KM2  = (PIXEL_RES ** 2) / 1e6   # 10m → 0.0001 km²

CHANGE_CLASSES = {0:"Stable Land", 1:"Water Gain", 2:"Water Loss", 3:"Stable Water"}
CHANGE_COLORS  = {0:"#d4d4d4", 1:"#d62728", 2:"#aec7e8", 3:"#08519c"}

meta = json.loads(META_FILE.read_text())
print(f"Before : {meta['before']['date']}  After : {meta['after']['date']}")
print(f"Model  : {meta['classification']['model']}")

with open(Path("config.toml"), "rb") as f:
    cfg = tomllib.load(f)

DATA_ROOT  = Path(cfg["data"]["root"])
MASKS_DIR  = DATA_ROOT / "masks"
CHANGE_DIR = DATA_ROOT / "change"
META_FILE  = DATA_ROOT / "scene_metadata.json"
CHANGE_DIR.mkdir(parents=True, exist_ok=True)

PIXEL_RES  = cfg["processing"]["pixel_res_m"]
CROP_SIZE  = cfg["processing"]["crop_size"]
PIXEL_KM2  = (PIXEL_RES ** 2) / 1e6   # 10m → 0.0001 km²

CHANGE_CLASSES = {0:"Stable Land", 1:"Water Gain", 2:"Water Loss", 3:"Stable Water"}
CHANGE_COLORS  = {0:"#d4d4d4", 1:"#d62728", 2:"#aec7e8", 3:"#08519c"}

meta = json.loads(META_FILE.read_text())
print(f"Before : {meta['before']['date']}  After : {meta['after']['date']}")
print(f"Model  : {meta['classification']['model']}")

Before : 2024-08-28  After : 2024-10-22
Model  : facebook/sam2.1-hiera-base-plus

In [3]:

Copied!

meta
meta

Out[3]:

{'before': {'id': 'S2B_17SLV_20240828_0_L2A',
  'date': '2024-08-28',
  'tile': '17SLV',
  'cloud_cover': 16.833709,
  'bbox': [-82.8, 35.4, -82.3, 35.8]},
 'after': {'id': 'S2A_17SLV_20241022_0_L2A',
  'date': '2024-10-22',
  'tile': '17SLV',
  'cloud_cover': 0.001261,
  'bbox': [-82.8, 35.4, -82.3, 35.8]},
 'catalog': 'aws-element84',
 'collection': 'sentinel-2-l2a',
 'preprocessing': {'bands': ['blue', 'green', 'red', 'nir08', 'swir16'],
  'cloud_pct_before': 4.25,
  'cloud_pct_after': 0.0},
 'classification': {'ndwi_thresh_before': -0.5824,
  'ndwi_thresh_after': -0.4963,
  'model': 'facebook/sam2.1-hiera-base-plus',
  'crop_size': 2048,
  'chip_size': 512,
  'water_seg_threshold': 0.0,
  'water_pct': {'ndwi_before': 14.8772,
   'ndwi_after': 10.6904,
   'sam_before': 0.1239,
   'sam_after': 0.5572}}}

2. Load Water Masks from Notebook 03¶

In [4]:

Copied!





def load(fname, chunks=None):
    kw = dict(chunks=chunks) if chunks else {}
    return xr.open_dataarray(MASKS_DIR / fname, **kw).squeeze(drop=True)

# NDWI — full scene
water_b_ndwi = load("before_water_ndwi.nc", {"y":1024,"x":1024})
water_a_ndwi = load("after_water_ndwi.nc",  {"y":1024,"x":1024})

# SAM — center crop
water_b_sam  = load("before_water_sam.nc")
water_a_sam  = load("after_water_sam.nc")

print("NDWI masks — shape:", dict(water_b_ndwi.sizes))
print("SAM  masks — shape:", dict(water_b_sam.sizes))
def load(fname, chunks=None):
    kw = dict(chunks=chunks) if chunks else {}
    return xr.open_dataarray(MASKS_DIR / fname, **kw).squeeze(drop=True)

# NDWI — full scene
water_b_ndwi = load("before_water_ndwi.nc", {"y":1024,"x":1024})
water_a_ndwi = load("after_water_ndwi.nc",  {"y":1024,"x":1024})

# SAM — center crop
water_b_sam  = load("before_water_sam.nc")
water_a_sam  = load("after_water_sam.nc")

print("NDWI masks — shape:", dict(water_b_ndwi.sizes))
print("SAM  masks — shape:", dict(water_b_sam.sizes))

NDWI masks — shape: {'y': 4510, 'x': 4600}
SAM  masks — shape: {'y': 2048, 'x': 2048}

3. Spatial Alignment Check¶

In [5]:

Copied!





def aligned(a, b, tol=1.0):
    return (a.sizes == b.sizes and
            a.rio.crs == b.rio.crs and
            np.allclose(a.x.values, b.x.values, atol=tol) and
            np.allclose(a.y.values, b.y.values, atol=tol))

# NDWI
if aligned(water_b_ndwi, water_a_ndwi):
    print("NDWI: already aligned ✓")
    water_a_ndwi_al = water_a_ndwi
else:
    print("NDWI: reprojecting after → before …")
    water_a_ndwi_al = water_a_ndwi.rio.reproject_match(
        water_b_ndwi, resampling=rioxarray.enums.Resampling.nearest
    )

# SAM
if aligned(water_b_sam, water_a_sam):
    print("SAM : already aligned ✓")
    water_a_sam_al = water_a_sam
else:
    print("SAM : reprojecting after → before …")
    water_a_sam_al = water_a_sam.rio.reproject_match(
        water_b_sam, resampling=rioxarray.enums.Resampling.nearest
    )
def aligned(a, b, tol=1.0):
    return (a.sizes == b.sizes and
            a.rio.crs == b.rio.crs and
            np.allclose(a.x.values, b.x.values, atol=tol) and
            np.allclose(a.y.values, b.y.values, atol=tol))

# NDWI
if aligned(water_b_ndwi, water_a_ndwi):
    print("NDWI: already aligned ✓")
    water_a_ndwi_al = water_a_ndwi
else:
    print("NDWI: reprojecting after → before …")
    water_a_ndwi_al = water_a_ndwi.rio.reproject_match(
        water_b_ndwi, resampling=rioxarray.enums.Resampling.nearest
    )

# SAM
if aligned(water_b_sam, water_a_sam):
    print("SAM : already aligned ✓")
    water_a_sam_al = water_a_sam
else:
    print("SAM : reprojecting after → before …")
    water_a_sam_al = water_a_sam.rio.reproject_match(
        water_b_sam, resampling=rioxarray.enums.Resampling.nearest
    )

NDWI: already aligned ✓
SAM : already aligned ✓

4. Change Maps¶

In [6]:

Copied!





def change_map(before, after, name):
    """Encode 4-class change: before×2 + after (0–3)."""
    cm = (before * 2 + after).rename(name)
    return cm.compute()

change_ndwi = change_map(water_b_ndwi, water_a_ndwi_al, "change_ndwi")
change_sam  = change_map(water_b_sam,  water_a_sam_al,  "change_sam")
print("Change maps computed.")

def class_stats(change_da, pixel_km2=PIXEL_KM2):
    valid = ~np.isnan(change_da.values)
    n = valid.sum()
    rows = []
    for cls_id, cls_name in CHANGE_CLASSES.items():
        count = int((change_da.values[valid] == cls_id).sum())
        rows.append({
            "Class": cls_name,
            "Pixels": count,
            "Area km²": round(count * pixel_km2, 3),
            "%": round(100.0*count/n, 3) if n else 0.0,
        })
    return pd.DataFrame(rows)

print("\n── NDWI Change (full scene) ──────────────────────────────────────────")
df_ndwi = class_stats(change_ndwi)
print(df_ndwi.to_string(index=False))

print("\n── SAM Change (center crop) ──────────────────────────────────────────")
df_sam = class_stats(change_sam)
print(df_sam.to_string(index=False))
def change_map(before, after, name):
    """Encode 4-class change: before×2 + after (0–3)."""
    cm = (before * 2 + after).rename(name)
    return cm.compute()

change_ndwi = change_map(water_b_ndwi, water_a_ndwi_al, "change_ndwi")
change_sam  = change_map(water_b_sam,  water_a_sam_al,  "change_sam")
print("Change maps computed.")

def class_stats(change_da, pixel_km2=PIXEL_KM2):
    valid = ~np.isnan(change_da.values)
    n = valid.sum()
    rows = []
    for cls_id, cls_name in CHANGE_CLASSES.items():
        count = int((change_da.values[valid] == cls_id).sum())
        rows.append({
            "Class": cls_name,
            "Pixels": count,
            "Area km²": round(count * pixel_km2, 3),
            "%": round(100.0*count/n, 3) if n else 0.0,
        })
    return pd.DataFrame(rows)

print("\n── NDWI Change (full scene) ──────────────────────────────────────────")
df_ndwi = class_stats(change_ndwi)
print(df_ndwi.to_string(index=False))

print("\n── SAM Change (center crop) ──────────────────────────────────────────")
df_sam = class_stats(change_sam)
print(df_sam.to_string(index=False))

Change maps computed.

── NDWI Change (full scene) ──────────────────────────────────────────
       Class   Pixels  Area km²      %
 Stable Land 17913801  1791.380 90.184
  Water Gain   282046    28.205  1.420
  Water Loss   760259    76.026  3.827
Stable Water   907433    90.743  4.568

── SAM Change (center crop) ──────────────────────────────────────────
       Class  Pixels  Area km²      %
 Stable Land 4168722   416.872 99.390
  Water Gain   20386     2.039  0.486
  Water Loss    2213     0.221  0.053
Stable Water    2983     0.298  0.071

5. Key Findings¶

In [7]:

Copied!





def _safe(v):
    return None if (v is None or (isinstance(v, float) and (math.isnan(v) or math.isinf(v)))) else round(v, 4)

def summary(df, label):
    s = df.set_index("Class")["Area km²"]
    wg, wl, sw = s["Water Gain"], s["Water Loss"], s["Stable Water"]
    total_b = wg + sw
    net     = wg - wl
    pct     = 100.0 * net / total_b if total_b > 0 else float("nan")
    print(f"\n{label}")
    print(f"  Water Gain : {wg:.3f} km²   Water Loss: {wl:.3f} km²   Stable Water: {sw:.3f} km²")
    print(f"  Net change : {net:+.3f} km²  ({pct:+.1f}% relative to before)")
    return {"water_gain_km2": _safe(wg), "water_loss_km2": _safe(wl),
            "stable_water_km2": _safe(sw), "net_change_km2": _safe(net), "net_pct": _safe(pct)}

print("=" * 60)
print("  CHANGE DETECTION SUMMARY")
print(f"  {meta['before']['date']} → {meta['after']['date']}")
print("=" * 60)
ndwi_summary = summary(df_ndwi, "NDWI (full scene):")
sam_summary  = summary(df_sam,  f"SAM ({cfg['model']['name']}) (crop):")
def _safe(v):
    return None if (v is None or (isinstance(v, float) and (math.isnan(v) or math.isinf(v)))) else round(v, 4)

def summary(df, label):
    s = df.set_index("Class")["Area km²"]
    wg, wl, sw = s["Water Gain"], s["Water Loss"], s["Stable Water"]
    total_b = wg + sw
    net     = wg - wl
    pct     = 100.0 * net / total_b if total_b > 0 else float("nan")
    print(f"\n{label}")
    print(f"  Water Gain : {wg:.3f} km²   Water Loss: {wl:.3f} km²   Stable Water: {sw:.3f} km²")
    print(f"  Net change : {net:+.3f} km²  ({pct:+.1f}% relative to before)")
    return {"water_gain_km2": _safe(wg), "water_loss_km2": _safe(wl),
            "stable_water_km2": _safe(sw), "net_change_km2": _safe(net), "net_pct": _safe(pct)}

print("=" * 60)
print("  CHANGE DETECTION SUMMARY")
print(f"  {meta['before']['date']} → {meta['after']['date']}")
print("=" * 60)
ndwi_summary = summary(df_ndwi, "NDWI (full scene):")
sam_summary  = summary(df_sam,  f"SAM ({cfg['model']['name']}) (crop):")

============================================================
  CHANGE DETECTION SUMMARY
  2024-08-28 → 2024-10-22
============================================================

NDWI (full scene):
  Water Gain : 28.205 km²   Water Loss: 76.026 km²   Stable Water: 90.743 km²
  Net change : -47.821 km²  (-40.2% relative to before)

SAM (facebook/sam2.1-hiera-base-plus) (crop):
  Water Gain : 2.039 km²   Water Loss: 0.221 km²   Stable Water: 0.298 km²
  Net change : +1.818 km²  (+77.8% relative to before)

6. Visualisation¶

In [8]:

Copied!





opts = dict(x="x", y="y", rasterize=True, geo=True,
            cmap=[CHANGE_COLORS[i] for i in range(4)],
            clim=(-0.5, 3.5), colorbar=True,
            clabel="0=Stable Land  1=Water Gain  2=Water Loss  3=Stable Water")

p_ndwi = change_ndwi.hvplot.image(
    title=f"NDWI Change  {meta['before']['date']} → {meta['after']['date']}",
    width=700, height=550, **opts)
p_sam  = change_sam.hvplot.image(
    title=f"SAM Change (crop)  {meta['before']['date']} → {meta['after']['date']}",
    width=600, height=550, **opts)

(p_ndwi + p_sam).cols(2)
opts = dict(x="x", y="y", rasterize=True, geo=True,
            cmap=[CHANGE_COLORS[i] for i in range(4)],
            clim=(-0.5, 3.5), colorbar=True,
            clabel="0=Stable Land  1=Water Gain  2=Water Loss  3=Stable Water")

p_ndwi = change_ndwi.hvplot.image(
    title=f"NDWI Change  {meta['before']['date']} → {meta['after']['date']}",
    width=700, height=550, **opts)
p_sam  = change_sam.hvplot.image(
    title=f"SAM Change (crop)  {meta['before']['date']} → {meta['after']['date']}",
    width=600, height=550, **opts)

(p_ndwi + p_sam).cols(2)

WARNING:param.Image00347: Image dimension(s) x and y are not evenly sampled to relative tolerance of 0.001. Please use the QuadMesh element for irregularly sampled data or set a higher tolerance on hv.config.image_rtol or the rtol parameter in the Image constructor.
WARNING:param.Image00397: Image dimension(s) x and y are not evenly sampled to relative tolerance of 0.001. Please use the QuadMesh element for irregularly sampled data or set a higher tolerance on hv.config.image_rtol or the rtol parameter in the Image constructor.
WARNING:param.Image00347: Image dimension(s) x and y are not evenly sampled to relative tolerance of 0.001. Please use the QuadMesh element for irregularly sampled data or set a higher tolerance on hv.config.image_rtol or the rtol parameter in the Image constructor.
WARNING:param.Image00397: Image dimension(s) x and y are not evenly sampled to relative tolerance of 0.001. Please use the QuadMesh element for irregularly sampled data or set a higher tolerance on hv.config.image_rtol or the rtol parameter in the Image constructor.

Out[8]:

BokehModel(combine_events=True, render_bundle={'docs_json': {'11fd14e8-5bee-4186-be59-9f01d16ba07b': {'version…

7. Save Results¶

In [9]:

Copied!





change_ndwi.astype("float32").to_netcdf(CHANGE_DIR / "change_ndwi.nc")
print(f"Saved change_ndwi.nc  {(CHANGE_DIR/'change_ndwi.nc').stat().st_size/1e6:.0f} MB")
change_sam.astype("float32").to_netcdf(CHANGE_DIR / "change_sam.nc")
print(f"Saved change_sam.nc   {(CHANGE_DIR/'change_sam.nc').stat().st_size/1e6:.0f} MB")

meta["change_detection"] = {
    "ndwi":  ndwi_summary,
    "sam":   sam_summary,
    "model": cfg["model"]["name"],
}
META_FILE.write_text(json.dumps(meta, indent=2))
print(f"Updated {META_FILE}")
print("\n✅ Notebook 04 complete.")
change_ndwi.astype("float32").to_netcdf(CHANGE_DIR / "change_ndwi.nc")
print(f"Saved change_ndwi.nc  {(CHANGE_DIR/'change_ndwi.nc').stat().st_size/1e6:.0f} MB")
change_sam.astype("float32").to_netcdf(CHANGE_DIR / "change_sam.nc")
print(f"Saved change_sam.nc   {(CHANGE_DIR/'change_sam.nc').stat().st_size/1e6:.0f} MB")

meta["change_detection"] = {
    "ndwi":  ndwi_summary,
    "sam":   sam_summary,
    "model": cfg["model"]["name"],
}
META_FILE.write_text(json.dumps(meta, indent=2))
print(f"Updated {META_FILE}")
print("\n✅ Notebook 04 complete.")

Saved change_ndwi.nc  83 MB
Saved change_sam.nc   17 MB
Updated /home/bny/data/sentinel2_helene/scene_metadata.json

✅ Notebook 04 complete.