Source code: Earthquake damage mapping with Sentinel-1 in Python: log ratio and z-scores in Antakya

"""Earthquake damage mapping with Sentinel-1 in Python: log ratio, pre-event z-scores and an agreement table.

Streams Sentinel-1 RTC (gamma0 VV and VH, 10 m) for Antakya, Hatay, Turkey, from the Microsoft Planetary
Computer (no account or key): a pre-event stack from October 2022 to 5 February 2023 and the first
post-event acquisition after the 6 February 2023 Kahramanmaras earthquakes, on three relative orbits.
Computes the log ratio and a z-score against the pre-event stack, masks to ESA WorldCover built-up land,
aggregates to a 200 m grid and scores the result against the HOT/OpenStreetMap "destroyed buildings"
reference (ODbL) with ROC AUC and a confusion table.
Run with:  python damage_mapping.py
"""
import io
import json
import platform
import time
import zipfile
from collections import defaultdict
from concurrent.futures import ThreadPoolExecutor
from pathlib import Path

import matplotlib

matplotlib.use("Agg")
import geopandas as gpd
import matplotlib.dates as mdates
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import planetary_computer
import pystac_client
import rasterio
import requests
import scipy
from matplotlib.colors import ListedColormap
from matplotlib.patches import Patch
from rasterio.enums import Resampling
from rasterio.transform import from_origin
from rasterio.vrt import WarpedVRT
from rasterio.warp import transform_bounds
from scipy import stats

HERE = Path(__file__).parent
OUT, DATA = HERE / "output", HERE / "data"
for d in (OUT, DATA):
    d.mkdir(exist_ok=True)
plt.rcParams.update({"figure.dpi": 150, "savefig.bbox": "tight", "font.size": 9})
pd.set_option("display.width", 140)

print("=== versions ===")
print(f"python {platform.python_version()} | numpy {np.__version__} | scipy {scipy.__version__} | "
      f"rasterio {rasterio.__version__} | geopandas {gpd.__version__} | pystac-client {pystac_client.__version__}")


def section(title: str) -> None:
    print(f"\n=== {title} ===", flush=True)


# --- Step 1: pre- and post-event Sentinel-1 RTC on a common 10 m grid -----------------
BBOX = (36.10, 36.15, 36.22, 36.25)        # lon/lat: Antakya (Hatay), Turkey
EVENT = pd.Timestamp("2023-02-06 01:17", tz="UTC")
ORBITS = {14: "ascending", 21: "descending", 94: "descending"}
CRS = "EPSG:32637"                          # UTM 37N
b = [round(v, -1) for v in transform_bounds(4326, CRS, *BBOX)]
TR = from_origin(b[0], b[3], 10, 10)
W, H = int((b[2] - b[0]) / 10), int((b[3] - b[1]) / 10)

section("Step 1: search and read Sentinel-1 RTC")
catalog = pystac_client.Client.open("https://planetarycomputer.microsoft.com/api/stac/v1",
                                    modifier=planetary_computer.sign_inplace)
print(f"grid: {W} x {H} pixels of 10 m in {CRS}")


def read_band(item, band: str) -> np.ndarray:
    """Warp one RTC asset onto the common grid; 0 is nodata in the COGs."""
    with rasterio.open(item.assets[band].href) as src, \
            WarpedVRT(src, crs=CRS, transform=TR, width=W, height=H, nodata=0,
                      resampling=Resampling.bilinear) as vrt:
        a = vrt.read(1).astype("float32")
    a[a <= 0] = np.nan
    return a


def load_orbit(orbit: int) -> dict:
    """Return {"dates", "vv" (t, y, x), "vh" (t, y, x)} for one relative orbit, cached as .npz."""
    cache = DATA / f"s1_orbit{orbit:03d}.npz"
    if not cache.exists():
        items = list(catalog.search(collections=["sentinel-1-rtc"], bbox=BBOX, datetime="2022-10-01/2023-02-19",
                                    query={"sat:relative_orbit": {"eq": orbit}}).items())
        jobs = [(i, band) for i in items for band in ("vv", "vh")]
        with ThreadPoolExecutor(8) as ex:
            arrays = list(ex.map(lambda j: read_band(*j), jobs))
        by_date = defaultdict(lambda: defaultdict(list))  # several slices can cover the town on one pass
        for (item, band), a in zip(jobs, arrays):
            by_date[item.datetime.strftime("%Y-%m-%dT%H:%M")][band].append(a)
        dates = sorted(by_date)
        merged = {band: np.stack([np.nanmax(np.stack(by_date[d][band]), 0) for d in dates]) for band in ("vv", "vh")}
        np.savez_compressed(cache, dates=np.array(dates), **merged)
    z = np.load(cache)
    return {"dates": pd.to_datetime(z["dates"]).tz_localize("UTC"), "vv": z["vv"], "vh": z["vh"]}


t0 = time.perf_counter()
s1 = {orbit: load_orbit(orbit) for orbit in ORBITS}
rows = []
for orbit, d in s1.items():
    pre = d["dates"][d["dates"] < EVENT]
    post = d["dates"][d["dates"] > EVENT]
    rows.append({"orbit": orbit, "pass": ORBITS[orbit], "pre_scenes": len(pre), "first_pre": f"{pre[0]:%d %b %Y}",
                 "last_pre": f"{pre[-1]:%d %b %Y}", "post_scenes": len(post), "first_post": f"{post[0]:%d %b %Y}",
                 "valid_%": round(100 * np.isfinite(d['vv'][-1]).mean(), 1)})
print(pd.DataFrame(rows).to_string(index=False))
print(f"read in {time.perf_counter() - t0:.1f} s")


# --- Step 2: the reference: destroyed buildings mapped in OpenStreetMap (HOT export) ---
section("Step 2: download the HOT/OSM destroyed-buildings reference")
REF = DATA / "hotosm_tur_destroyed_buildings.gpkg"
if not REF.exists():
    pkg = requests.get("https://data.humdata.org/api/3/action/package_show",
                       params={"id": "hotosm_tur_destroyed_buildings"}, timeout=60).json()["result"]
    url = next(r["url"] for r in pkg["resources"] if r["url"].endswith(".zip") and "gpkg" in r["url"])
    with zipfile.ZipFile(io.BytesIO(requests.get(url, timeout=300).content)) as zf:
        name = next(n for n in zf.namelist() if n.endswith(".gpkg"))
        REF.write_bytes(zf.read(name))
ref = gpd.read_file(REF)
print(f"{len(ref)} destroyed-building polygons in Turkey | sources: {ref['source'].value_counts().to_dict()}")
ref = ref.to_crs(CRS)
ref["geometry"] = ref.geometry.centroid
x0, y0 = TR.c, TR.f
ref = ref[(ref.geometry.x > x0) & (ref.geometry.x < x0 + W * 10) & (ref.geometry.y < y0) & (ref.geometry.y > y0 - H * 10)]
print(f"{len(ref)} of them inside the Antakya study area")

# --- Step 3: a built-up mask from ESA WorldCover 2021 ---------------------------------
section("Step 3: ESA WorldCover 2021 built-up mask")
wc_cache = DATA / "worldcover_2021.npy"
if not wc_cache.exists():
    wc_item = next(i for i in catalog.search(collections=["esa-worldcover"], bbox=BBOX).items() if "2021" in i.id)
    with rasterio.open(wc_item.assets["map"].href) as src, \
            WarpedVRT(src, crs=CRS, transform=TR, width=W, height=H, resampling=Resampling.nearest) as vrt:
        np.save(wc_cache, vrt.read(1))
wc = np.load(wc_cache)
built = wc == 50
print(f"built-up: {built.mean():.1%} of the study area ({built.sum() * 100 / 1e6:.1f} km2)")

# --- Step 4: log ratio and z-score against the pre-event stack ------------------------
section("Step 4: log ratio and pre-event z-score per orbit and polarisation")
change = {}
rows = []
for orbit, d in s1.items():
    pre = np.asarray(d["dates"] < EVENT)
    i_post = int(np.flatnonzero(~pre)[0])
    for pol in ("vv", "vh"):
        db = 10 * np.log10(d[pol])
        stack_db, post_db = db[pre], db[i_post]
        log_ratio = post_db - stack_db[-1]                 # dB difference = 10 log10(post / pre)
        mu, sd = stack_db.mean(0), stack_db.std(0, ddof=1)
        z = (post_db - mu) / sd
        z_pre = (stack_db[-1] - stack_db[:-1].mean(0)) / stack_db[:-1].std(0, ddof=1)  # a no-event control
        change[(orbit, pol)] = {"log_ratio": log_ratio, "z": z, "z_control": z_pre}
        rows.append({"orbit": orbit, "pol": pol.upper(), "post": f"{d['dates'][i_post]:%d %b}",
                     "built_sd_dB": round(float(np.nanmedian(sd[built])), 2),
                     "|LR|>3dB_%": round(100 * float(np.nanmean(np.abs(log_ratio[built]) > 3)), 1),
                     "|z|>2_%": round(100 * float(np.nanmean(np.abs(z[built]) > 2)), 1),
                     "|z|>2_%_control": round(100 * float(np.nanmean(np.abs(z_pre[built]) > 2)), 1)})
print(pd.DataFrame(rows).to_string(index=False))
print("control = last pre-event date scored against the earlier pre-event dates (no earthquake in between)")
LC = {10: "tree cover", 30: "grassland", 40: "cropland", 50: "built-up"}
rows = []
for code, name in LC.items():
    mask = wc == code
    ev = [np.nanmean(np.abs(c["z"][mask]) > 2) for c in change.values()]
    ct = [np.nanmean(np.abs(c["z_control"][mask]) > 2) for c in change.values()]
    rows.append({"land cover": name, "pixels": int(mask.sum()), "|z|>2 % event": round(100 * float(np.mean(ev)), 1),
                 "|z|>2 % control": round(100 * float(np.mean(ct)), 1)})
print("\nexceedance by WorldCover class (mean of the 6 orbit/polarisation pairs):")
print(pd.DataFrame(rows).to_string(index=False))

# --- Step 5: aggregate to a 200 m grid and attach the reference -----------------------
section("Step 5: aggregate to a 200 m grid")


def grid_stats(cell_px: int):
    h, w = H // cell_px, W // cell_px

    def agg(a):
        a = np.where(built, a, np.nan)[:h * cell_px, :w * cell_px].reshape(h, cell_px, w, cell_px)
        return np.nanmean(a, axis=(1, 3))

    built_frac = built[:h * cell_px, :w * cell_px].reshape(h, cell_px, w, cell_px).mean(axis=(1, 3))
    cols = ((ref.geometry.x - TR.c) / (10 * cell_px)).astype(int).to_numpy()
    rws = ((TR.f - ref.geometry.y) / (10 * cell_px)).astype(int).to_numpy()
    ok = (rws < h) & (cols < w)
    count = np.zeros((h, w))
    np.add.at(count, (rws[ok], cols[ok]), 1)
    z_index = np.mean([agg((np.abs(c["z"]) > 2).astype("float32")) for c in change.values()], axis=0)
    lr_index = np.mean([agg(np.abs(c["log_ratio"])) for c in change.values()], axis=0)
    ctrl_index = np.mean([agg((np.abs(c["z_control"]) > 2).astype("float32")) for c in change.values()], axis=0)
    per_combo = {k: agg((np.abs(c["z"]) > 2).astype("float32")) for k, c in change.items()}
    return built_frac, count, z_index, lr_index, ctrl_index, per_combo


def auc(score, label):
    """ROC AUC via the Mann-Whitney U statistic (probability a damaged cell outranks an undamaged one)."""
    r = stats.rankdata(score)
    n1, n0 = label.sum(), (~label).sum()
    return float((r[label].sum() - n1 * (n1 + 1) / 2) / (n1 * n0))


CELL = 20                                       # 20 x 10 m = 200 m
built_frac, count, z_index, lr_index, ctrl_index, per_combo = grid_stats(CELL)
sel = built_frac >= 0.5                         # score only cells that are at least half built-up
label = count[sel] >= 1
print(f"{sel.sum()} cells of 200 m are >= 50 % built-up; {label.sum()} contain >= 1 destroyed building "
      f"({int(count[sel].sum())} of {len(ref)} reference buildings fall in scored cells)")
print(f"damage index (mean share of built pixels with |z| > 2 over 6 orbit/polarisation pairs): "
      f"median {np.median(z_index[sel]):.3f}, 95th pct {np.percentile(z_index[sel], 95):.3f}")

section("Step 6: how well does each score rank the reference cells? (ROC AUC, 0.5 = random)")
rows = [{"score": "|log ratio|, mean of 6", "AUC": auc(lr_index[sel], label)},
        {"score": "share |z|>2, mean of 6", "AUC": auc(z_index[sel], label)}]
rows += [{"score": f"share |z|>2, orbit {o} {p.upper()}", "AUC": auc(v[sel], label)} for (o, p), v in per_combo.items()]
auc_table = pd.DataFrame(rows)
print(auc_table.round(3).to_string(index=False))
rho, p_rho = stats.spearmanr(z_index[sel], count[sel])
print(f"Spearman rho (damage index vs number of destroyed buildings per cell) = {rho:.3f}, p = {p_rho:.1e}")
sens = []
for cell_px in (10, 20, 40):
    bf, cnt, zi, li, _, _ = grid_stats(cell_px)
    s = bf >= 0.5
    for k in (1, 3):
        lab = cnt[s] >= k
        sens.append({"cell_m": cell_px * 10, "damaged_if_>=": k, "cells": int(s.sum()), "damaged": int(lab.sum()),
                     "AUC_logratio": round(auc(li[s], lab), 3), "AUC_zindex": round(auc(zi[s], lab), 3)})
print("\nsensitivity to cell size and reference threshold:")
print(pd.DataFrame(sens).to_string(index=False))

# --- Step 7: threshold without the reference, then build the agreement table ---------
section("Step 7: unsupervised thresholds and the agreement table")


def otsu(values: np.ndarray, bins: int = 64) -> float:
    """Otsu's threshold: maximise the between-class variance of a histogram."""
    hist, edges = np.histogram(values, bins)
    centres = (edges[:-1] + edges[1:]) / 2
    best, best_var = centres[0], -1.0
    for k in range(1, bins):
        w0, w1 = hist[:k].sum(), hist[k:].sum()
        if w0 == 0 or w1 == 0:
            continue
        m0, m1 = (hist[:k] * centres[:k]).sum() / w0, (hist[k:] * centres[k:]).sum() / w1
        if w0 * w1 * (m0 - m1) ** 2 > best_var:
            best, best_var = centres[k], w0 * w1 * (m0 - m1) ** 2
    return float(best)


def agreement(pred: np.ndarray, lab: np.ndarray) -> dict:
    tp, fp = int((pred & lab).sum()), int((pred & ~lab).sum())
    fn, tn = int((~pred & lab).sum()), int((~pred & ~lab).sum())
    n = tp + fp + fn + tn
    po = (tp + tn) / n
    pe = ((tp + fp) * (tp + fn) + (fn + tn) * (fp + tn)) / n**2
    return {"tp": tp, "fp": fp, "fn": fn, "tn": tn, "flagged": tp + fp, "precision": round(tp / max(tp + fp, 1), 3),
            "recall": round(tp / (tp + fn), 3), "f1": round(2 * tp / (2 * tp + fp + fn), 3),
            "overall_accuracy": round(po, 3), "kappa": round((po - pe) / (1 - pe), 3)}


v = z_index[sel]
thresholds = {"control 99th pct": float(np.percentile(ctrl_index[sel], 99)),
              "Otsu on event index": otsu(v),
              "event mean + 1 sd": float(v.mean() + v.std()),
              "event mean + 2 sd": float(v.mean() + 2 * v.std())}
rows = [{"rule": k, "threshold": round(t, 3), **agreement(v > t, label)} for k, t in thresholds.items()]
print(pd.DataFrame(rows).drop(columns=["tp", "fp", "fn", "tn"]).to_string(index=False))
threshold = thresholds["Otsu on event index"]
pred = v > threshold
m = agreement(pred, label)
tp, fp, fn, tn = m["tp"], m["fp"], m["fn"], m["tn"]
print(f"\nconfusion table at the Otsu threshold ({threshold:.3f}), 200 m built-up cells:")
conf = pd.DataFrame([[tp, fn, tp + fn], [fp, tn, fp + tn], [tp + fp, fn + tn, tp + fp + fn + tn]],
                    index=["reference: >= 1 destroyed", "reference: none mapped", "total"],
                    columns=["SAR: damaged", "SAR: not flagged", "total"])
print(conf.to_string())
print(f"odds ratio {(tp * tn) / (fp * fn):.1f}; Fisher exact p = {stats.fisher_exact([[tp, fn], [fp, tn]])[1]:.1e}")
metrics = {"threshold_otsu": round(threshold, 3), "cells": len(v), **{k: m[k] for k in
           ("tp", "fp", "fn", "tn", "precision", "recall", "f1", "overall_accuracy", "kappa")}}
# --- Step 8: figures and metrics ------------------------------------------------------
section("Step 8: figures")
EXT = (0, W / 100, 0, H / 100)                         # km from the south-west corner
km = (ref.geometry.x - TR.c) / 1000, (ref.geometry.y - (TR.f - H * 10)) / 1000
d21 = s1[21]
pre21 = np.asarray(d21["dates"] < EVENT)
i21 = int(np.flatnonzero(~pre21)[0])

fig, ax = plt.subplots(1, 3, figsize=(13, 4.6), sharey=True, layout="constrained")
panels = [(10 * np.log10(np.nanmean(d21["vv"][pre21], 0)), "gray", (-20, 5), "pre-event mean VV, 11 dates"),
          (10 * np.log10(d21["vv"][i21]), "gray", (-20, 5), f"post-event VV, {d21['dates'][i21]:%d %b %Y}"),
          (change[(21, "vv")]["log_ratio"], "RdBu_r", (-6, 6), "log ratio, 10 Feb vs 29 Jan")]
for a, (img, cmap, (lo, hi), title) in zip(ax, panels):
    im = a.imshow(img, cmap=cmap, vmin=lo, vmax=hi, extent=EXT)
    a.set_title(title)
    a.set_xlabel("km east")
    fig.colorbar(im, ax=a, shrink=0.8, label="dB")
ax[0].set_ylabel("km north")
fig.suptitle("Sentinel-1 RTC, relative orbit 21 (descending), Antakya")
fig.savefig(OUT / "01-pre-post-logratio.png")
plt.close(fig)

fig, ax = plt.subplots(1, 2, figsize=(11, 5), sharey=True)
for a, key, title in [(ax[0], "z_control", "control: 29 Jan vs 10 earlier dates (no event)"),
                      (ax[1], "z", "event: 10 Feb vs 11 pre-event dates")]:
    z = np.where(built, change[(21, "vh")][key], np.nan)
    a.imshow(np.where(built, np.nan, 1), cmap="Greys", vmin=0, vmax=4, extent=EXT)
    im = a.imshow(z, cmap="RdBu_r", vmin=-4, vmax=4, extent=EXT)
    a.scatter(*km, s=2, c="k", lw=0, label="destroyed building (OSM)")
    a.set_title(title)
    a.set_xlabel("km east")
ax[0].set_ylabel("km north")
ax[1].legend(loc="lower right", markerscale=4, fontsize=8)
fig.colorbar(im, ax=ax, shrink=0.8, label="z-score, VH (built-up pixels only)")
fig.savefig(OUT / "02-zscore-control-vs-event.png")
plt.close(fig)

# per-date exceedance: each pre-event date scored against the other pre-event dates (leave one out)
h, w = H // CELL, W // CELL
cnt_px = np.kron(count, np.ones((CELL, CELL)))
cnt_px = np.pad(cnt_px, ((0, H - cnt_px.shape[0]), (0, W - cnt_px.shape[1])))
in_sel = np.kron(sel, np.ones((CELL, CELL))).astype(bool)
in_sel = np.pad(in_sel, ((0, H - in_sel.shape[0]), (0, W - in_sel.shape[1])))
grp_dmg, grp_ok = built & in_sel & (cnt_px >= 3), built & in_sel & (cnt_px == 0)
fig, ax = plt.subplots(1, 3, figsize=(12, 3.4), sharey=True, layout="constrained")
for a, orbit in zip(ax, ORBITS):
    d = s1[orbit]
    db = 10 * np.log10(d["vv"])
    pre = np.flatnonzero(np.asarray(d["dates"] < EVENT))
    i_post = int(np.flatnonzero(np.asarray(d["dates"] > EVENT))[0])
    series = {"dmg": [], "ok": []}
    for t in list(pre) + [i_post]:
        others = [i for i in pre if i != t]
        z = (db[t] - db[others].mean(0)) / db[others].std(0, ddof=1)
        series["dmg"].append(np.nanmean(np.abs(z[grp_dmg]) > 2))
        series["ok"].append(np.nanmean(np.abs(z[grp_ok]) > 2))
    dates = d["dates"][list(pre) + [i_post]]
    a.plot(dates, series["dmg"], "o-", c="#d7191c", ms=3, label=">= 3 destroyed buildings in cell")
    a.plot(dates, series["ok"], "o-", c="#2c7bb6", ms=3, label="no destroyed building mapped")
    a.axvline(EVENT, c="k", ls="--", lw=0.8)
    a.set_title(f"orbit {orbit} ({ORBITS[orbit]}), VV")
    a.xaxis.set_major_locator(mdates.MonthLocator(bymonth=[10, 12, 2]))
    a.xaxis.set_major_formatter(mdates.DateFormatter("%b %Y"))
ax[0].set_ylabel("share of built pixels with |z| > 2")
ax[0].legend(fontsize=7, loc="upper left")
fig.savefig(OUT / "03-exceedance-time-series.png")
plt.close(fig)

fig, ax = plt.subplots(1, 2, figsize=(11, 5), sharey=True)
gext = (0, w * CELL / 100, 0, h * CELL / 100)
im0 = ax[0].imshow(np.where(sel, z_index, np.nan), cmap="YlOrRd", vmin=0.1, vmax=0.45, extent=gext)
ax[0].set_title("SAR damage index (200 m cells >= 50 % built-up)")
fig.colorbar(im0, ax=ax[0], shrink=0.75, label="share of built pixels with |z| > 2")
im1 = ax[1].imshow(np.where(sel, count, np.nan), cmap="Purples", vmin=0, vmax=15, extent=gext)
ax[1].set_title("reference: destroyed buildings per cell (OSM)")
fig.colorbar(im1, ax=ax[1], shrink=0.75, label="buildings", extend="max")
for a in ax:
    a.set_facecolor("#e6e6e6")
    a.set_xlabel("km east")
ax[0].set_ylabel("km north")
fig.savefig(OUT / "04-damage-index-vs-reference.png")
plt.close(fig)


def roc(score, lab):
    order = np.argsort(-score)
    l = lab[order]
    return np.r_[0, np.cumsum(~l) / (~l).sum()], np.r_[0, np.cumsum(l) / l.sum()]


fig, ax = plt.subplots(figsize=(4.6, 4.3))
for score, name, col in [(z_index[sel], "share |z|>2 (pre-event stack)", "#d7191c"),
                         (lr_index[sel], "|log ratio| (one pre date)", "#2c7bb6")]:
    fpr, tpr = roc(score, label)
    ax.plot(fpr, tpr, c=col, label=f"{name}: AUC {auc(score, label):.2f}")
ax.plot([0, 1], [0, 1], c="0.6", ls=":", label="random")
ax.plot(fp / (fp + tn), tp / (tp + fn), "ko", label="Otsu threshold")
ax.set_xlabel("false positive rate (cells with no mapped damage)")
ax.set_ylabel("true positive rate (cells with destroyed buildings)")
ax.legend(fontsize=7, loc="lower right")
fig.savefig(OUT / "05-roc.png")
plt.close(fig)

cat_grid = np.full((h, w), np.nan)
p_full = np.zeros((h, w), bool)
p_full[sel] = pred
l_full = count >= 1
cat_grid[sel & p_full & l_full] = 3
cat_grid[sel & p_full & ~l_full] = 2
cat_grid[sel & ~p_full & l_full] = 1
cat_grid[sel & ~p_full & ~l_full] = 0
cols4 = ["#f0f0f0", "#2c7bb6", "#fdae61", "#d7191c"]
fig, ax = plt.subplots(figsize=(6, 5.5))
ax.set_facecolor("#ffffff")
ax.imshow(cat_grid, cmap=ListedColormap(cols4), vmin=-0.5, vmax=3.5, extent=gext, interpolation="nearest")
ax.legend(handles=[Patch(color=c, label=f"{l} ({v})") for c, l, v in
                   zip(cols4[::-1], ["both: damaged", "SAR only", "reference only", "neither"], [tp, fp, fn, tn])],
          loc="lower right", fontsize=8)
ax.set_title(f"agreement at the Otsu threshold (kappa {metrics['kappa']:.2f})")
ax.set_xlabel("km east")
ax.set_ylabel("km north")
fig.savefig(OUT / "06-agreement-map.png")
plt.close(fig)

metrics.update({"auc_zindex": round(auc(z_index[sel], label), 3), "auc_logratio": round(auc(lr_index[sel], label), 3),
                "spearman_rho": round(float(rho), 3), "reference_buildings_in_area": int(len(ref))})
(OUT / "metrics.json").write_text(json.dumps(metrics, indent=2))
print("saved", ", ".join(sorted(p.name for p in OUT.glob("*.png"))))
print(f"total run time {time.perf_counter() - t0:.1f} s")