Source code: Sentinel-1 soil moisture in Python: change detection checked against COSMOS-UK

"""Sentinel-1 soil moisture in Python: TU Wien-style change detection checked against COSMOS-UK.

Streams two years (2023-2024) of Sentinel-1 RTC (gamma0 VV, 10 m) over four COSMOS-UK soil moisture
stations in lowland England (two arable, two grassland) from the Microsoft Planetary Computer (no account
or key), builds per-orbit dry and wet references from the backscatter time series, converts backscatter to
relative surface soil moisture (0-100 %), and compares it with the cosmic-ray (COSMOS) and 10 cm probe
soil moisture and rainfall from the UKCEH COSMOS-UK API (Open Government Licence). Finishes with 100 m
relative soil moisture maps around Wimpole, Cambridgeshire, for a dry and a wet date.
Run with:  python soil_moisture.py
"""
import platform
import time
from collections import defaultdict
from concurrent.futures import ThreadPoolExecutor
from pathlib import Path

import matplotlib

matplotlib.use("Agg")
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 pyproj import Transformer
from rasterio.enums import Resampling
from rasterio.transform import from_origin
from rasterio.vrt import WarpedVRT
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"pandas {pd.__version__} | rasterio {rasterio.__version__} | pystac-client {pystac_client.__version__}")


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


SITES = {  # COSMOS-UK station code: (name, land cover, lat, lon)
    "WIMPL": ("Wimpole", "arable", 52.1322, -0.0444),
    "ROTHD": ("Rothamsted", "arable", 51.8136, -0.3778),
    "CARDT": ("Cardington", "grassland", 52.1056, -0.4249),
    "WADDN": ("Waddesdon", "grassland", 51.8389, -0.9478),
}
START, END = "2023-01-01", "2024-12-31"
BNG = "EPSG:27700"                       # British National Grid: one metric CRS for every orbit and site
to_bng = Transformer.from_crs(4326, BNG, always_xy=True)

# --- Step 1: in-situ soil moisture and rain from the COSMOS-UK API --------------------
section("Step 1: COSMOS-UK daily soil moisture and rainfall")


def cosmos_daily(site: str) -> pd.DataFrame:
    """Daily COSMOS volumetric water content (%), 10 cm TDT probe (%) and rain (mm) as a DataFrame."""
    cache = DATA / f"cosmos_{site}.csv"
    if not cache.exists():
        r = requests.get(f"https://cosmos-api.ceh.ac.uk/collections/1D/locations/{site}", timeout=120,
                         params={"datetime": f"{START}T00:00:00Z/{END}T00:00:00Z",
                                 "parameter-name": "cosmos_vwc,tdt1_vwc,precip"})
        r.raise_for_status()
        cov = r.json()["coverages"][0]
        df = pd.DataFrame({p: cov["ranges"][p]["values"] for p in ("cosmos_vwc", "tdt1_vwc", "precip")},
                          index=pd.to_datetime(cov["domain"]["axes"]["t"]["values"]).date)
        df.index.name = "date"
        df.to_csv(cache)
    return pd.read_csv(cache, index_col="date", parse_dates=True).astype(float)


insitu = {s: cosmos_daily(s) for s in SITES}
rows = []
for s, df in insitu.items():
    rows.append({"site": s, "name": SITES[s][0], "cover": SITES[s][1], "days": len(df),
                 "cosmos_vwc_min": df.cosmos_vwc.min(), "cosmos_vwc_max": df.cosmos_vwc.max(),
                 "tdt_10cm_min": df.tdt1_vwc.min(), "tdt_10cm_max": df.tdt1_vwc.max(),
                 "rain_mm": round(df.precip.sum())})
print(pd.DataFrame(rows).round(1).to_string(index=False))

# --- Step 2: Sentinel-1 backscatter time series over each station ---------------------
section("Step 2: Sentinel-1 RTC VV time series over each station footprint")
catalog = pystac_client.Client.open("https://planetarycomputer.microsoft.com/api/stac/v1",
                                    modifier=planetary_computer.sign_inplace)
HALF = 200                               # 400 m x 400 m box: roughly the COSMOS footprint


def read_mean(href: str, transform, width: int, height: int, resampling=Resampling.nearest) -> np.ndarray:
    with rasterio.open(href) as src, WarpedVRT(src, crs=BNG, transform=transform, width=width, height=height,
                                               nodata=0, resampling=resampling) as vrt:
        a = vrt.read(1).astype("float32")
    a[a <= 0] = np.nan
    return a


def site_series(site: str) -> pd.DataFrame:
    """Footprint-mean linear gamma0 VV per acquisition (slices of one pass are merged)."""
    cache = DATA / f"s1_{site}.csv"
    if not cache.exists():
        _, _, lat, lon = SITES[site]
        x, y = to_bng.transform(lon, lat)
        tr = from_origin(x - HALF, y + HALF, 10, 10)
        items = list(catalog.search(collections=["sentinel-1-rtc"], intersects={"type": "Point", "coordinates": [lon, lat]},
                                    datetime=f"{START}/{END}").items())
        with ThreadPoolExecutor(8) as ex:
            chips = list(ex.map(lambda i: read_mean(i.assets["vv"].href, tr, 2 * HALF // 10, 2 * HALF // 10), items))
        recs = defaultdict(list)
        for item, a in zip(items, chips):
            if np.isfinite(a).mean() > 0.9:
                recs[(item.datetime.strftime("%Y-%m-%d %H:%M"), item.properties["sat:relative_orbit"],
                      item.properties["sat:orbit_state"], item.properties["platform"])].append(float(np.nanmean(a)))
        df = pd.DataFrame([{"time": k[0], "orbit": k[1], "pass": k[2], "platform": k[3], "vv": np.mean(v)}
                           for k, v in recs.items()]).sort_values("time")
        df.to_csv(cache, index=False)
    df = pd.read_csv(cache, parse_dates=["time"])
    df["date"] = df.time.dt.normalize()
    df["vv_db"] = 10 * np.log10(df.vv)
    return df


t0 = time.perf_counter()
s1 = {s: site_series(s) for s in SITES}
print(f"read in {time.perf_counter() - t0:.1f} s")
for s, df in s1.items():
    per_orbit = df.groupby(["orbit", "pass"]).vv_db.agg(["count", "mean", "std", "min", "max"]).round(2)
    print(f"\n{s} ({SITES[s][0]}, {SITES[s][1]}): {len(df)} acquisitions")
    print(per_orbit.to_string())

# --- Step 3: change detection: dry and wet references per orbit -----------------------
section("Step 3: dry/wet references and relative surface soil moisture")
LO, HI = 5, 95                           # percentiles of the backscatter record used as dry and wet references


def change_detection(df: pd.DataFrame, by_orbit: bool = True) -> pd.Series:
    """Relative soil moisture (%) = (sigma0 - dry) / (wet - dry), references taken per orbit or pooled."""
    groups = df.groupby("orbit").vv_db if by_orbit else df.assign(all=1).groupby("all").vv_db
    dry = groups.transform(lambda v: np.percentile(v, LO))
    wet = groups.transform(lambda v: np.percentile(v, HI))
    return (100 * (df.vv_db - dry) / (wet - dry)).clip(0, 100)


rows = []
for s, df in s1.items():
    df["ssm"] = change_detection(df)
    df["ssm_pooled"] = change_detection(df, by_orbit=False)
    ref = df.groupby("orbit").vv_db.agg(dry=lambda v: np.percentile(v, LO), wet=lambda v: np.percentile(v, HI))
    for orbit, r in ref.iterrows():
        rows.append({"site": s, "orbit": orbit, "dry_dB": round(r.dry, 2), "wet_dB": round(r.wet, 2),
                     "sensitivity_dB": round(r.wet - r.dry, 2)})
print(pd.DataFrame(rows).to_string(index=False))

# --- Step 4: compare with COSMOS-UK ---------------------------------------------------
section("Step 4: agreement with COSMOS-UK (Pearson r, same-day match)")
matched = {}
rows = []
for s, df in s1.items():
    m = df.merge(insitu[s], left_on="date", right_index=True, how="inner").dropna(subset=["cosmos_vwc"])
    m["rain_3d"] = [insitu[s].precip.loc[d - pd.Timedelta(days=2):d].sum() for d in m.date]
    matched[s] = m
    winter = m.date.dt.month.isin([10, 11, 12, 1, 2, 3])
    rows.append({"site": s, "cover": SITES[s][1], "n": len(m),
                 "r_cosmos": stats.pearsonr(m.ssm, m.cosmos_vwc)[0],
                 "r_tdt10cm": stats.pearsonr(m.dropna(subset=["tdt1_vwc"]).ssm, m.dropna(subset=["tdt1_vwc"]).tdt1_vwc)[0],
                 "r_pooled_orbits": stats.pearsonr(m.ssm_pooled, m.cosmos_vwc)[0],
                 "r_raw_dB": stats.pearsonr(m.vv_db, m.cosmos_vwc)[0],
                 "r_Oct-Mar": stats.pearsonr(m[winter].ssm, m[winter].cosmos_vwc)[0],
                 "r_Apr-Sep": stats.pearsonr(m[~winter].ssm, m[~winter].cosmos_vwc)[0],
                 "rmse_vwc_%": np.nan})
    # an absolute estimate: rescale relative SSM to the station's own observed range (needs in-situ data!)
    lo, hi = np.percentile(m.cosmos_vwc, [LO, HI])
    est = lo + m.ssm / 100 * (hi - lo)
    rows[-1]["rmse_vwc_%"] = float(np.sqrt(np.mean((est - m.cosmos_vwc) ** 2)))
agree = pd.DataFrame(rows)
print(agree.round(2).to_string(index=False))
print("r_pooled_orbits: one dry/wet reference for all orbits (no incidence-angle handling)")
print("rmse_vwc_%: relative SSM rescaled to the station's 5th-95th percentile VWC range")

rows = []
for s, m in matched.items():
    wet_after_rain = m[m.rain_3d >= 10]
    dry_spell = m[m.rain_3d == 0]
    rows.append({"site": s, "n_rain>=10mm_3d": len(wet_after_rain), "ssm_after_rain": wet_after_rain.ssm.mean(),
                 "n_no_rain_3d": len(dry_spell), "ssm_no_rain": dry_spell.ssm.mean(),
                 "cosmos_after_rain": wet_after_rain.cosmos_vwc.mean(), "cosmos_no_rain": dry_spell.cosmos_vwc.mean()})
print("\nmean relative SSM (%) and COSMOS VWC (%) by rain in the previous 3 days:")
print(pd.DataFrame(rows).round(1).to_string(index=False))

# --- Step 5: from the surface skin to the root zone: the soil water index -------------
section("Step 5: soil water index (exponential filter) vs COSMOS")


def soil_water_index(times: pd.Series, ssm: pd.Series, t_days: float) -> np.ndarray:
    """Recursive exponential filter (Albergel et al. 2008): a smoothed, lagged, deeper-layer proxy."""
    t = (times - times.iloc[0]).dt.total_seconds().to_numpy() / 86400
    swi, gain = np.empty(len(ssm)), 1.0
    swi[0] = ssm.iloc[0]
    for n in range(1, len(ssm)):
        gain = gain / (gain + np.exp(-(t[n] - t[n - 1]) / t_days))
        swi[n] = swi[n - 1] + gain * (ssm.iloc[n] - swi[n - 1])
    return swi


T_VALUES = (1, 5, 10, 20, 40, 60)
rows = []
for s, m in matched.items():
    row = {"site": s, "cover": SITES[s][1], "r_SSM": stats.pearsonr(m.ssm, m.cosmos_vwc)[0]}
    for T in T_VALUES:
        m[f"swi_{T}"] = soil_water_index(m.time, m.ssm, T)
        row[f"r_SWI_T{T}"] = stats.pearsonr(m[f"swi_{T}"], m.cosmos_vwc)[0]
    rows.append(row)
swi_table = pd.DataFrame(rows)
print(swi_table.round(2).to_string(index=False))
T_BEST = int(swi_table[[f"r_SWI_T{T}" for T in T_VALUES]].mean().idxmax().split("T")[-1])
print(f"best characteristic time length on average: T = {T_BEST} days")

# --- Step 6: relative soil moisture maps around Wimpole -------------------------------
section("Step 6: 100 m relative soil moisture maps around Wimpole (orbit 81, descending)")
MAP_ORBIT, MAP_RES = 81, 100
cx, cy = to_bng.transform(SITES["WIMPL"][3], SITES["WIMPL"][2])
MAP_HALF = 6000
map_tr = from_origin(round(cx, -2) - MAP_HALF, round(cy, -2) + MAP_HALF, MAP_RES, MAP_RES)
MW = 2 * MAP_HALF // MAP_RES
map_cache = DATA / f"map_orbit{MAP_ORBIT}.npz"
if not map_cache.exists():
    lon0, lat0 = SITES["WIMPL"][3], SITES["WIMPL"][2]
    items = list(catalog.search(collections=["sentinel-1-rtc"], bbox=(lon0 - 0.1, lat0 - 0.06, lon0 + 0.1, lat0 + 0.06),
                                datetime=f"{START}/{END}", query={"sat:relative_orbit": {"eq": MAP_ORBIT}}).items())
    with ThreadPoolExecutor(8) as ex:   # GDAL "average" resampling = 10 x 10 multilook in linear units
        chips = list(ex.map(lambda i: read_mean(i.assets["vv"].href, map_tr, MW, MW, Resampling.average), items))
    by_date = defaultdict(list)
    for item, a in zip(items, chips):
        by_date[item.datetime.strftime("%Y-%m-%d")].append(a)
    dates = sorted(by_date)
    np.savez_compressed(map_cache, dates=np.array(dates),
                        vv=np.stack([np.nanmax(np.stack(by_date[d]), 0) for d in dates]))
    wc_tiles = []                        # the map straddles 0 degrees longitude, so it needs two 3 x 3 degree tiles
    for wc_item in catalog.search(collections=["esa-worldcover"], bbox=(lon0 - 0.1, lat0 - 0.06, lon0 + 0.1, lat0 + 0.06),
                                  query={"esa_worldcover:product_version": {"eq": "2.0.0"}}).items():
        with rasterio.open(wc_item.assets["map"].href) as src, WarpedVRT(src, crs=BNG, transform=map_tr, width=MW,
                                                                        height=MW, resampling=Resampling.mode) as vrt:
            wc_tiles.append(vrt.read(1))
    np.save(DATA / "worldcover_map.npy", np.max(wc_tiles, axis=0))
z = np.load(map_cache)
map_dates, cube = pd.to_datetime(z["dates"]), 10 * np.log10(z["vv"])
wc = np.load(DATA / "worldcover_map.npy")
soil = np.isin(wc, [30, 40])            # WorldCover grassland and cropland: where the model is meant to work
dry_ref, wet_ref = np.nanpercentile(cube, LO, axis=0), np.nanpercentile(cube, HI, axis=0)
ssm_cube = np.clip(100 * (cube - dry_ref) / (wet_ref - dry_ref), 0, 100)
cos_w = insitu["WIMPL"].cosmos_vwc
on_dates = cos_w.reindex(map_dates.normalize()).to_numpy()
rain5 = insitu["WIMPL"].precip.rolling(5, min_periods=1).sum().reindex(map_dates.normalize()).to_numpy()
dry_i = int(np.nanargmin(np.where(rain5 == 0, on_dates, np.nan)))   # driest COSMOS day after 5 rain-free days
winter = np.isin(map_dates.month, [10, 11, 12, 1, 2, 3])
wet_i = int(np.nanargmax(np.where(winter, on_dates, np.nan)))            # wettest day with bare or low crops
veg_i = int(np.nanargmax(np.where(np.isin(map_dates.month, [4, 5, 6]), on_dates, np.nan)))  # wet, crops growing
print(f"{len(map_dates)} dates on a {MW} x {MW} grid of {MAP_RES} m ({np.isfinite(cube).mean():.1%} valid); "
      f"grassland+cropland = {soil.mean():.0%} of pixels")
print(f"median sensitivity (wet - dry) on cropland {np.nanmedian((wet_ref - dry_ref)[wc == 40]):.2f} dB, "
      f"grassland {np.nanmedian((wet_ref - dry_ref)[wc == 30]):.2f} dB, tree cover {np.nanmedian((wet_ref - dry_ref)[wc == 10]):.2f} dB, "
      f"built-up {np.nanmedian((wet_ref - dry_ref)[wc == 50]):.2f} dB")
rows = []
for name, i in (("dry", dry_i), ("wet (winter)", wet_i), ("wet (crops growing)", veg_i)):
    rows.append({"case": name, "date": f"{map_dates[i]:%d %b %Y}", "COSMOS_vwc_%": on_dates[i],
                 "rain_prev_5d_mm": rain5[i],
                 "median_SSM_cropland_%": float(np.nanmedian(ssm_cube[i][wc == 40])),
                 "median_SSM_grassland_%": float(np.nanmedian(ssm_cube[i][wc == 30]))})
print(pd.DataFrame(rows).round(1).to_string(index=False))
daily_median = pd.Series(np.nanmedian(ssm_cube[:, soil], axis=1), index=map_dates.normalize())
both = pd.concat([daily_median.rename("map_median"), cos_w], axis=1, join="inner").dropna()
print(f"area-median SSM (grass+crop) vs Wimpole COSMOS over {len(both)} dates: "
      f"r = {stats.pearsonr(both.map_median, both.cosmos_vwc)[0]:.2f}")

# --- Step 7: figures ------------------------------------------------------------------
section("Step 7: figures")
ORBIT_COLOURS = {30: "#1a9641", 59: "#d7191c", 81: "#2c7bb6", 132: "#fdae61", 154: "#7b3294"}

df, ins = s1["WIMPL"], insitu["WIMPL"]
fig, ax = plt.subplots(2, 1, figsize=(10, 5.6), sharex=True, height_ratios=[1.3, 1])
for orbit, g in df.groupby("orbit"):
    c = ORBIT_COLOURS[orbit]
    ax[0].plot(g.time, g.vv_db, "o-", ms=2.5, lw=0.6, c=c, label=f"orbit {orbit} ({g['pass'].iloc[0]})")
    ax[0].axhline(np.percentile(g.vv_db, LO), c=c, ls="--", lw=0.8)
    ax[0].axhline(np.percentile(g.vv_db, HI), c=c, ls=":", lw=0.8)
ax[0].set_ylabel("gamma0 VV (dB)")
ax[0].set_title("Wimpole (arable): backscatter per orbit; dashed = dry reference (5th pct), dotted = wet (95th pct)")
ax[0].legend(fontsize=7, ncol=3, loc="lower left")
ax[1].bar(ins.index, ins.precip, width=1, color="#abd9e9", label="rain (mm/day)")
ax[1].set_ylabel("rain (mm/day)")
a2 = ax[1].twinx()
a2.plot(ins.index, ins.cosmos_vwc, c="k", lw=1, label="COSMOS VWC")
a2.plot(ins.index, ins.tdt1_vwc, c="#d7191c", lw=0.8, alpha=0.8, label="probe at 10 cm")
a2.set_ylabel("soil moisture (% vol)")
a2.legend(fontsize=7, loc="lower left")
ax[1].xaxis.set_major_formatter(mdates.DateFormatter("%b %y"))
fig.savefig(OUT / "01-backscatter-orbits-wimpole.png")
plt.close(fig)

fig, ax = plt.subplots(len(SITES), 1, figsize=(10, 9), sharex=True)
for a, (s, m) in zip(ax, matched.items()):
    ins = insitu[s]
    for orbit, g in m.groupby("orbit"):
        a.plot(g.time, g.ssm, "o", ms=2.2, c=ORBIT_COLOURS[orbit], alpha=0.8)
    a.plot(m.time, m[f"swi_{T_BEST}"], c="#7b3294", lw=1.2, label=f"SWI (T = {T_BEST} d)")
    a.set_ylim(-5, 105)
    a.set_ylabel("relative SSM (%)")
    a2 = a.twinx()
    a2.plot(ins.index, ins.cosmos_vwc, c="k", lw=1, label="COSMOS VWC")
    a2.set_ylabel("VWC (%)")
    r = agree.set_index("site").loc[s]
    a.set_title(f"{s} {SITES[s][0]} ({SITES[s][1]}): r(SSM, COSMOS) = {r.r_cosmos:.2f}, "
                f"r(SWI, COSMOS) = {swi_table.set_index('site').loc[s, f'r_SWI_T{T_BEST}']:.2f}", fontsize=9)
ax[0].legend(loc="lower left", fontsize=7)
a2.legend(loc="lower right", fontsize=7)
ax[-1].xaxis.set_major_formatter(mdates.DateFormatter("%b %y"))
fig.savefig(OUT / "02-ssm-vs-cosmos-timeseries.png")
plt.close(fig)

fig, ax = plt.subplots(1, len(SITES), figsize=(12, 3.3), sharey=True)
for a, (s, m) in zip(ax, matched.items()):
    winter = m.date.dt.month.isin([10, 11, 12, 1, 2, 3])
    a.scatter(m[winter].cosmos_vwc, m[winter].ssm, s=8, c="#2c7bb6", label="Oct-Mar")
    a.scatter(m[~winter].cosmos_vwc, m[~winter].ssm, s=8, c="#fdae61", label="Apr-Sep")
    r = agree.set_index("site").loc[s]
    a.set_title(f"{s} ({SITES[s][1]})\nr Oct-Mar {r['r_Oct-Mar']:.2f} | Apr-Sep {r['r_Apr-Sep']:.2f}", fontsize=8)
    a.set_xlabel("COSMOS VWC (%)")
ax[0].set_ylabel("Sentinel-1 relative SSM (%)")
ax[0].legend(fontsize=7)
fig.savefig(OUT / "03-scatter-by-season.png")
plt.close(fig)

fig, axs = plt.subplots(2, 2, figsize=(10.5, 9.4), layout="constrained")
ax = axs.ravel()
ext = (0, 2 * MAP_HALF / 1000, 0, 2 * MAP_HALF / 1000)
for a, i in zip(ax[:3], (dry_i, wet_i, veg_i)):
    a.imshow(np.where(soil, np.nan, 1), cmap="Greys", vmin=0, vmax=3, extent=ext)
    im = a.imshow(np.where(soil, ssm_cube[i], np.nan), cmap="RdYlBu", vmin=0, vmax=100, extent=ext)
    a.plot(MAP_HALF / 1000, MAP_HALF / 1000, "k^", ms=7, mec="w")
    a.set_title(f"{map_dates[i]:%d %b %Y}\nCOSMOS {on_dates[i]:.0f} % vol, 5-day rain {rain5[i]:.0f} mm")
    fig.colorbar(im, ax=a, shrink=0.85, label="relative SSM (%), grass and crop")
im = ax[3].imshow(np.where(np.isin(wc, [10, 30, 40, 50]), wet_ref - dry_ref, np.nan), cmap="viridis", vmin=1, vmax=10, extent=ext)
ax[3].set_title("sensitivity\nwet minus dry reference (dB)")
fig.colorbar(im, ax=ax[3], shrink=0.85, label="dB")
for a in ax:
    a.set_xlabel("km")
    a.set_ylabel("km")
fig.savefig(OUT / "04-ssm-maps-wimpole.png")
plt.close(fig)

fig, ax = plt.subplots(figsize=(7, 3.4))
cols = ["r_raw_dB", "r_pooled_orbits", "r_cosmos", "r_tdt10cm"]
labels = ["raw dB", "SSM, pooled refs", "SSM, per-orbit refs", "SSM vs 10 cm probe"]
x = np.arange(len(SITES))
for k, (c, lab, colr) in enumerate(zip(cols + ["swi"], labels + [f"SWI T={T_BEST} d"],
                                       ["#bababa", "#abd9e9", "#2c7bb6", "#fdae61", "#d7191c"])):
    vals = agree[c] if c != "swi" else swi_table[f"r_SWI_T{T_BEST}"]
    ax.bar(x + (k - 2) * 0.16, vals, 0.16, color=colr, label=lab)
ax.set_xticks(x, [f"{s}\n{SITES[s][1]}" for s in SITES])
ax.set_ylabel("Pearson r with COSMOS VWC")
ax.axhline(0, c="k", lw=0.6)
ax.legend(fontsize=7, ncol=3, loc="upper left")
ax.set_ylim(0, 1)
fig.savefig(OUT / "05-correlation-summary.png")
plt.close(fig)
print("saved", ", ".join(sorted(p.name for p in OUT.glob("*.png"))))
print(f"total run time {time.perf_counter() - t0:.1f} s")