Source code: InSAR time series in Python: Svartsengi inflation cycles validated with GNSS

"""Volcano deformation with InSAR, part 3: an SBAS time series of the Svartsengi inflation cycles,
validated against GNSS.

Streams every short-baseline (<= 24 day) COMET LiCSAR Sentinel-1 interferogram on ascending frame
016A_02504_162118 from October 2023 to July 2026, multilooks them, checks the network and the
phase-closure (unwrapping) errors, inverts the stack pixel by pixel with a small-baseline (SBAS)
least-squares solver, and compares the LOS time series at the SENG GNSS station with the GNSS
east/north/up time series (Nevada Geodetic Laboratory) projected into the radar line of sight.
It also maps the coherence loss over new lava flows.
Run with:  python insar_gnss_timeseries.py   (first run downloads about 780 files; reruns use data/)
"""
import io
import os
import platform
import re
import time
import warnings
from concurrent.futures import ThreadPoolExecutor
from datetime import datetime
from importlib.metadata import version
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 rasterio
import requests
from cmcrameri import cm
from rasterio.windows import from_bounds

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,
                     "figure.facecolor": "white", "savefig.facecolor": "white"})
pd.set_option("display.width", 120)
warnings.filterwarnings("ignore", category=RuntimeWarning)  # all-NaN blocks over the sea

print("=== versions ===")
print(f"python {platform.python_version()} | rasterio {rasterio.__version__} | numpy {np.__version__} | "
      f"pandas {pd.__version__} | matplotlib {matplotlib.__version__} | cmcrameri {version('cmcrameri')}")


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


LICSAR = "https://gws-access.jasmin.ac.uk/public/nceo_geohazards/LiCSAR_products.public"
FRAME = "016A_02504_162118"                    # Sentinel-1 ascending track 16 over Reykjanes
BASE = f"{LICSAR}/16/{FRAME}"
START, END, MAX_DAYS = "20240201", "20260730", 24   # short-baseline network
BBOX = (-22.80, 63.72, -22.05, 64.05)          # same window as parts 1 and 2
ML = (2, 4)                                    # multilook 2 rows x 4 columns -> about 220 m x 195 m
SENG = (-22.4338, 63.8822)                     # GNSS station at the Svartsengi power station
REF_STATION = "NYLA"                           # GNSS station 18 km NW, inside the InSAR reference area
REF_KM = 15                                    # InSAR reference = land pixels > 15 km from SENG
MIN_COH = 0.25                                 # multilooked coherence needed to use a pixel in a pair
WAVELENGTH = 0.0554658                         # Sentinel-1 C-band wavelength in metres
RAD2MM = -WAVELENGTH / (4 * np.pi) * 1000      # LiCSAR convention: positive = towards the satellite
NGL = "https://geodesy.unr.edu/gps_timeseries/IGS20/tenv3/IGS20"
# Onset dates of the dyke intrusion and the nine Sundhnukur eruptions (Icelandic Met Office)
EVENTS = ["2023-11-10", "2023-12-18", "2024-01-14", "2024-02-08", "2024-03-16", "2024-05-29",
          "2024-08-22", "2024-11-20", "2025-04-01", "2025-07-16"]
ENDED = {"2025-07-16": "2025-08-05"}           # the 9th eruption lasted almost three weeks (IMO)
LAVA_EVENT = "2024-05-29"                      # eruption used to show coherence loss over new lava

# GDAL settings as environment variables so that they also apply in the download threads
os.environ.update(GDAL_DISABLE_READDIR_ON_OPEN="EMPTY_DIR", GDAL_HTTP_MAX_RETRY="4", GDAL_HTTP_RETRY_DELAY="2",
                  CPL_VSIL_CURL_ALLOWED_EXTENSIONS=".tif")


def days(pair: str) -> int:
    return (datetime.strptime(pair[9:], "%Y%m%d") - datetime.strptime(pair[:8], "%Y%m%d")).days


def read_window(url: str):
    """Read only the part of a remote LiCSAR GeoTIFF that covers BBOX (HTTP range requests)."""
    for attempt in range(4):
        try:
            with rasterio.open(url) as src:
                win = from_bounds(*BBOX, src.transform).round_offsets().round_lengths()
                return src.read(1, window=win).astype("float32"), src.window_transform(win)
        except rasterio.errors.RasterioIOError:
            if attempt == 3:
                raise
            time.sleep(5 * (attempt + 1))


def multilook(a: np.ndarray, valid: np.ndarray) -> np.ndarray:
    """Average ML[0] x ML[1] blocks, keeping a block only if at least half of its pixels are valid."""
    r, c = ML
    H, W = a.shape[0] // r * r, a.shape[1] // c * c
    a = np.where(valid, a, np.nan)[:H, :W].reshape(H // r, r, W // c, c)
    n = np.isfinite(a).sum((1, 3))
    return np.where(n >= r * c / 2, np.nanmean(a, (1, 3)), np.nan).astype("float32")


def read_pair(pair: str):
    """Multilooked unwrapped phase and coherence for one pair, cached in data/pairs/ (None if missing)."""
    f = DATA / "pairs" / f"{pair}.npz"
    if not f.exists():
        url = f"{BASE}/interferograms/{pair}/{pair}.geo"
        try:
            unw, t = read_window(f"{url}.unw.tif")
            coh, _ = read_window(f"{url}.cc.tif")
        except rasterio.errors.RasterioIOError:
            return None                                 # listed, but the files are not (yet) online
        land = unw != 0                                 # 0 = sea or not unwrapped by LiCSAR
        np.savez_compressed(f, unw=multilook(unw, land), coh=multilook(coh / 255, land), t=np.array(t)[:6])
    z = np.load(f)
    return z["unw"], z["coh"], rasterio.Affine(*z["t"])

def deramp(u: np.ndarray, use: np.ndarray, X: np.ndarray, Y: np.ndarray) -> np.ndarray:
    """Fit a + b*x + c*y to the far-field pixels and subtract it (orbit, ionosphere, long-wave troposphere)."""
    A = np.column_stack([np.ones(use.sum()), X[use], Y[use]])
    coef = np.linalg.lstsq(A, u[use], rcond=None)[0]
    return u - (coef[0] + coef[1] * X + coef[2] * Y)


# --- Step 1: download and multilook the short-baseline stack (cached) -------------------------
section(f"LiCSAR stack on {FRAME}")
html = requests.get(f"{BASE}/interferograms/", timeout=60).text
online = sorted(set(re.findall(r'href="(\d{8}_\d{8})/"', html)))
wanted = [p for p in online if p[:8] >= START and p[9:] <= END and days(p) <= MAX_DAYS]
(DATA / "pairs").mkdir(exist_ok=True)
t0 = time.time()
with ThreadPoolExecutor(6) as pool:                    # a few parallel HTTP readers
    res = dict(zip(wanted, pool.map(read_pair, wanted)))
pairs = [p for p in wanted if res[p] is not None]
print(f"{len(online):,} interferograms online; {len(wanted)} of <= {MAX_DAYS} days between {START} and {END}; "
      f"{len(pairs)} have unwrapped products ({time.time() - t0:.0f} s to read or load)")
unw = np.array([res[p][0] for p in pairs])
coh = np.array([res[p][1] for p in pairs])
t = res[pairs[0]][2]
transform = rasterio.Affine(t.a * ML[1], 0, t.c, 0, t.e * ML[0], t.f)

# --- Step 2: the network and its design matrix ------------------------------------------------
section("small-baseline network")
epochs = sorted({p[:8] for p in pairs} | {p[9:] for p in pairs})
idx = {e: i for i, e in enumerate(epochs)}
events = [e for e in EVENTS if epochs[0] < e.replace("-", "") < epochs[-1]]   # eruptions inside the series
A = np.zeros((len(pairs), len(epochs) - 1))         # unknowns: phase at each epoch after the first
for k, p in enumerate(pairs):
    A[k, idx[p[9:]] - 1] = 1
    if idx[p[:8]] > 0:
        A[k, idx[p[:8]] - 1] = -1
spans = pd.Series([days(p) for p in pairs]).value_counts().sort_index()
gaps = np.diff(pd.to_datetime(epochs).values).astype("timedelta64[D]").astype(int)
print(f"{len(pairs)} pairs between {epochs[0]} and {epochs[-1]} -> {len(epochs)} epochs")
print("pairs by time span (days): " + ", ".join(f"{d}: {n}" for d, n in spans.items()))
print(f"design matrix {A.shape}, rank {np.linalg.matrix_rank(A)} (needs {len(epochs) - 1}: one connected network)")
print(f"epoch spacing: median {np.median(gaps):.0f} days, longest {gaps.max()} days "
      f"({epochs[int(gaps.argmax())]} to {epochs[int(gaps.argmax()) + 1]})")
n_ifg, H, W = unw.shape
lon = transform.c + transform.a * (np.arange(W) + 0.5)
lat = transform.f + transform.e * (np.arange(H) + 0.5)
LON, LAT = np.meshgrid(lon, lat)
X = (LON - SENG[0]) * 111.32 * np.cos(np.radians(SENG[1]))   # km east of SENG
Y = (LAT - SENG[1]) * 110.6                                   # km north of SENG
land = np.isfinite(unw).any(0)
far = land & (np.hypot(X, Y) > REF_KM)
sr, sc = int((SENG[1] - transform.f) / transform.e), int((SENG[0] - transform.c) / transform.a)
section("multilooked stack")
print(f"{n_ifg} interferograms x {H} x {W} pixels ({transform.a:.3f} x {-transform.e:.3f} deg, about "
      f"{transform.a * 111.32 * np.cos(np.radians(63.9)) * 1000:.0f} m x {-transform.e * 110600:.0f} m)")
med_coh = np.array([np.nanmedian(c[land]) for c in coh])
print(f"median coherence per pair: {np.median(med_coh):.2f} (range {med_coh.min():.2f}-{med_coh.max():.2f}); "
      f"lowest {pairs[int(med_coh.argmin())]}")
month = pd.to_datetime([p[:8] for p in pairs]).month
print(f"median coherence in Dec-Mar pairs {np.median(med_coh[np.isin(month, [12, 1, 2, 3])]):.2f} "
      f"vs Jun-Sep pairs {np.median(med_coh[np.isin(month, [6, 7, 8, 9])]):.2f} (snow)")

# Reference every interferogram to the far field and remove a plane fitted there.
use_obs = np.isfinite(unw) & (coh >= MIN_COH)
for k in range(n_ifg):
    f = far & use_obs[k]
    unw[k] = deramp(unw[k], f, X, Y)
    unw[k] -= np.nanmedian(unw[k][f])
print(f"reference: {far.sum():,} land pixels > {REF_KM} km from SENG; plane removed from each pair")

fig, ax = plt.subplots(figsize=(10, 3.4))
t1 = pd.to_datetime([p[:8] for p in pairs])
t2 = pd.to_datetime([p[9:] for p in pairs])
for a, b, c, d in zip(t1, t2, med_coh, [days(p) for p in pairs]):
    ax.plot([a, b], [c, c], color=cm.batlow((d - 6) / 18), lw=1.2, alpha=0.9)
for e in events:
    ax.axvline(pd.Timestamp(e), color="0.35", lw=0.8, ls=":")
ax.set_ylabel("median coherence of the pair")
ax.set_title(f"Short-baseline network: {n_ifg} interferograms, {len(epochs)} Sentinel-1 epochs "
             f"(dotted lines = dyke intrusion and eruptions)")
sm = plt.cm.ScalarMappable(cmap=cm.batlow, norm=plt.Normalize(6, 24))
fig.colorbar(sm, ax=ax, label="time span (days)", pad=0.01)
ax.xaxis.set_major_formatter(mdates.DateFormatter("%b %Y"))
ax.grid(alpha=0.3)
fig.savefig(OUT / "01-network-coherence.png")
plt.close(fig)

# --- Step 3: phase closure: a test for unwrapping errors --------------------------------------
# For three pairs i->j, j->k and i->k the unwrapped phases must satisfy phi_ij + phi_jk - phi_ik = 0
# (apart from small non-zero closure from soil moisture and filtering). A multiple of 2*pi means at
# least one of the three was unwrapped wrongly. The plane removal is linear, so it closes too.
section("phase closure (unwrapping errors)")
pid = {p: k for k, p in enumerate(pairs)}
tri = [(pid[f"{a}_{b}"], pid[f"{b}_{c}"], pid[f"{a}_{c}"]) for a in epochs for b in epochs if b > a
       if f"{a}_{b}" in pid for c in epochs if c > b and f"{b}_{c}" in pid and f"{a}_{c}" in pid]
tri = np.array(tri)
clos = unw[tri[:, 0]] + unw[tri[:, 1]] - unw[tri[:, 2]]
ok3 = use_obs[tri[:, 0]] & use_obs[tri[:, 1]] & use_obs[tri[:, 2]]
bad = ok3 & (np.abs(clos) > np.pi)
frac_tri = bad.sum((1, 2)) / ok3.sum((1, 2))
frac_pix = np.where(ok3.sum(0) > 0, bad.sum(0) / np.maximum(ok3.sum(0), 1), np.nan)
frac_pix[~land] = np.nan
print(f"{len(tri)} closed triplets; |closure| > pi in {bad.sum() / ok3.sum():.2%} of valid pixel-triplets")
print(f"median closure magnitude {np.nanmedian(np.abs(clos[ok3])):.2f} rad "
      f"({np.nanmedian(np.abs(clos[ok3])) * abs(RAD2MM):.1f} mm)")
worst = np.argsort(frac_tri)[::-1][:3]
for w in worst:
    a, b, c = pairs[tri[w, 0]][:8], pairs[tri[w, 1]][:8], pairs[tri[w, 2]][9:]
    print(f"  worst triplet {a}-{b}-{c}: {frac_tri[w]:.1%} of pixels have a 2*pi closure error")
print(f"pixels with closure errors in > 10% of their triplets: {np.nanmean(frac_pix[land] > 0.1):.1%} of land")

# --- Step 4: SBAS inversion, pixel by pixel ---------------------------------------------------
# For each pixel, minimise sum_k w_k (phi_k - (x_j - x_i))^2 + lam^2 * |D x|^2, where w_k = 1 if the
# pixel is coherent and unwrapped in pair k and 0 otherwise. D is a second-difference operator; the
# tiny lam only bridges epochs that a pixel's own network leaves disconnected (rank deficiency).
section("SBAS inversion")
m = len(epochs) - 1
D = np.diff(np.eye(m + 1), 2, axis=0)[:, 1:]         # second difference of [0, x1, ..., xm]
LAM = 0.05
R = LAM ** 2 * D.T @ D
frac_valid = use_obs.mean(0)
solve = land & (frac_valid >= 0.5)
rows, cols = np.nonzero(solve)
obs = np.nan_to_num(unw[:, rows, cols])                 # n_ifg x n_pixels
w = use_obs[:, rows, cols].astype("float64")
AA = (A[:, :, None] * A[:, None, :]).reshape(n_ifg, -1)
x = np.zeros((m, len(rows)))
for s in range(0, len(rows), 1500):
    sl = slice(s, s + 1500)
    N = (w[:, sl].T @ AA).reshape(-1, m, m) + R
    b = (w[:, sl] * obs[:, sl]).T @ A
    x[:, sl] = np.linalg.solve(N, b[..., None])[..., 0].T
resid = (obs - A @ x) * w
n_obs = w.sum(0)
rms_pix = np.sqrt((resid ** 2).sum(0) / n_obs) * abs(RAD2MM)
ts = np.full((len(epochs), H, W), np.nan, "float32")
ts[0, rows, cols] = 0
ts[1:, rows, cols] = x * RAD2MM                          # cumulative LOS displacement in mm
print(f"solved {len(rows):,} of {land.sum():,} land pixels (coherent in >= 50% of pairs); "
      f"median {np.median(n_obs):.0f} pairs per pixel")
print(f"residual RMS per pixel: median {np.median(rms_pix):.1f} mm, 95th percentile {np.percentile(rms_pix, 95):.1f} mm")
print(f"pixels with a residual > pi in any pair: {(np.abs(resid) > np.pi).any(0).mean():.1%}")
inc = np.diff(ts[:, far & solve], axis=0)                 # far-field change from one epoch to the next
noise = 1.4826 * np.nanmedian(np.abs(inc - np.nanmedian(inc, 1, keepdims=True)), 1)   # robust std (MAD)
print(f"far-field scatter of epoch-to-epoch changes (atmosphere + noise, robust std): median {np.median(noise):.1f} mm, "
      f"max {noise.max():.1f} mm ({epochs[int(noise.argmax())]} to {epochs[int(noise.argmax()) + 1]})")

# --- Step 5: time-series maps ----------------------------------------------------------------
dt = pd.to_datetime(epochs)


def epoch_before(day: str) -> int:
    return int(np.nonzero(dt < pd.Timestamp(day))[0][-1])


def epoch_after(day: str) -> int:
    return int(np.nonzero(dt >= pd.Timestamp(day))[0][0])


panels = [(epoch_after("2024-11-20") + 1, epoch_before("2025-04-01"), "inflation between eruptions 7 and 8"),
          (epoch_after("2025-08-05"), len(epochs) - 1, "inflation after the 9th eruption"),
          (epoch_before("2025-04-01"), epoch_after("2025-04-01"), "across the 1 April 2025 eruption")]
fig, axes = plt.subplots(1, 3, figsize=(13, 4), sharey=True, layout="constrained")
extent = (X[0, 0], X[0, -1], Y[-1, 0], Y[0, 0])
section("cumulative LOS change at SENG")
for ax, (i, j, label) in zip(axes, panels):
    dmap = ts[j] - ts[i]
    lim = 150 if j - i > 3 else 80
    im = ax.imshow(dmap, cmap=cm.vik, vmin=-lim, vmax=lim, extent=extent, interpolation="nearest")
    ax.plot(0, 0, "k^", ms=6, mfc="none")
    ax.set_title(f"{label}\n{epochs[i]}\u2013{epochs[j]}")
    ax.set_xlabel("x, km east of SENG")
    fig.colorbar(im, ax=ax, shrink=0.8, label="LOS change (mm)")
    print(f"{label:40s} {epochs[i]}-{epochs[j]}: SENG {np.nanmedian(dmap[sr - 1:sr + 2, sc - 1:sc + 2]):+.0f} mm, "
          f"1st/99th percentile {np.nanpercentile(dmap, 1):+.0f}/{np.nanpercentile(dmap, 99):+.0f} mm")
axes[0].set_ylabel("y, km north of SENG")
fig.savefig(OUT / "02-cumulative-los-maps.png")
plt.close(fig)

# --- Step 6: GNSS -> LOS and the comparison at SENG --------------------------------------------
section(f"GNSS: SENG relative to {REF_STATION} (Nevada Geodetic Laboratory, IGS20)")


def load_ngl(station: str) -> pd.DataFrame:
    f = DATA / f"{station}.tenv3"
    if not f.exists():
        f.write_text(requests.get(f"{NGL}/{station}.tenv3", timeout=60).text)
    g = pd.read_csv(f, sep=r"\s+")
    g.columns = [c.strip("_").split("(")[0] for c in g.columns]
    g["date"] = pd.to_datetime(g["YYMMMDD"], format="%y%b%d")
    g = g.set_index("date")
    return pd.DataFrame({"E": (g["e0"] + g["east"]) * 1000, "N": (g["n0"] + g["north"]) * 1000,
                         "U": (g["u0"] + g["up"] + g["ant"]) * 1000})


seng, ref = load_ngl("SENG"), load_ngl(REF_STATION)
rel = (seng - ref).dropna().loc[dt[0] - pd.Timedelta(days=3):dt[-1] + pd.Timedelta(days=3)]
vec = {}
for k in "ENU":
    a, t = read_window(f"{BASE}/metadata/{FRAME}.geo.{k}.tif")
    vec[k] = float(a[int((SENG[1] - t.f) / t.e), int((SENG[0] - t.c) / t.a)])
print(f"{len(rel)} common daily solutions; LOS unit vector at SENG: E {vec['E']:+.3f} N {vec['N']:+.3f} U {vec['U']:+.3f}")
rel["LOS"] = vec["E"] * rel["E"] + vec["N"] * rel["N"] + vec["U"] * rel["U"]
g_ep = np.array([rel.loc[d - pd.Timedelta(days=1):d + pd.Timedelta(days=1), "LOS"].mean() for d in dt])
i_ep = np.array([np.nanmedian(ts[i, sr - 1:sr + 2, sc - 1:sc + 2]) for i in range(len(epochs))])
both = np.isfinite(g_ep) & np.isfinite(i_ep)
k0 = int(np.argmax(both))                                  # first epoch with both InSAR and GNSS
g_ep, i_ep = g_ep - g_ep[k0], i_ep - i_ep[k0]               # both relative to that epoch
diff = i_ep[both] - g_ep[both]
r = np.corrcoef(i_ep[both], g_ep[both])[0, 1]
print(f"InSAR at SENG: {np.isfinite(i_ep).sum()} epochs; GNSS on {both.sum()} of them")
print(f"total LOS change {epochs[k0]}-{epochs[-1]}: InSAR {i_ep[-1]:+.0f} mm, GNSS {g_ep[both][-1]:+.0f} mm")
print(f"raw InSAR - GNSS: RMS {np.sqrt(np.mean(diff ** 2)):.0f} mm, max {np.abs(diff).max():.0f} mm; r = {r:.3f}")

# Why they drift apart: the interferograms that span an eruption are often not unwrapped at SENG.
use_seng = use_obs[:, sr - 1:sr + 2, sc - 1:sc + 2].any((1, 2))
print("co-eruptive steps (last epoch before -> first epoch after each eruption onset):")
for e in events:
    i, j = epoch_before(e), epoch_after(e)
    span = [k for k, p in enumerate(pairs) if p[:8] < e.replace("-", "") <= p[9:]]
    print(f"  {e}: InSAR {i_ep[j] - i_ep[i]:+5.0f} mm, GNSS {g_ep[j] - g_ep[i]:+5.0f} mm; "
          f"{use_seng[span].sum()} of {len(span)} spanning pairs unwrapped and coherent at SENG")

# Validate the inflation phases: restart both series at the first epoch after each eruption.
section("InSAR vs GNSS between eruptions")
starts = [0] + [epoch_after(ENDED.get(e, e)) for e in events]
ends = [epoch_before(e) for e in events] + [len(epochs) - 1]


def by_phase(i_series: np.ndarray):
    """Restart InSAR and GNSS at zero at the start of each inflation phase and compare them."""
    seg = np.full(len(epochs), -1)
    i_seg, g_seg = np.full(len(epochs), np.nan), np.full(len(epochs), np.nan)
    out = []
    for n, (a0, a1) in enumerate(zip(starts, ends)):
        k = np.arange(a0, a1 + 1)
        k = k[both[k] & np.isfinite(i_series[k])]
        if len(k) < 3:
            continue
        seg[k] = n
        i_seg[k], g_seg[k] = i_series[k] - i_series[k[0]], g_ep[k] - g_ep[k[0]]
        yrs = (dt[k] - dt[k[0]]).days / 365.25
        out.append({"from": f"{dt[k[0]]:%Y-%m-%d}", "to": f"{dt[k[-1]]:%Y-%m-%d}", "epochs": len(k),
                    "InSAR_mm": i_seg[k[-1]], "GNSS_mm": g_seg[k[-1]],
                    "InSAR_mm_per_yr": np.polyfit(yrs, i_seg[k], 1)[0],
                    "GNSS_mm_per_yr": np.polyfit(yrs, g_seg[k], 1)[0]})
    ok = np.isfinite(i_seg) & np.isfinite(g_seg)
    d = i_seg[ok] - g_seg[ok]
    return pd.DataFrame(out), seg, i_seg, g_seg, np.sqrt(np.mean(d ** 2)), np.corrcoef(i_seg[ok], g_seg[ok])[0, 1]


inter, seg, i_seg, g_seg, rms2, r2 = by_phase(i_ep)
print(inter.round(0).to_string(index=False))
print(f"{np.isfinite(i_seg).sum()} epochs in {len(inter)} inflation phases: InSAR - GNSS RMS {rms2:.1f} mm, r = {r2:.3f}")
tot_i, tot_g = inter["InSAR_mm"].sum(), inter["GNSS_mm"].sum()
print(f"summed inflation-phase LOS uplift: InSAR {tot_i:+.0f} mm, GNSS {tot_g:+.0f} mm ({tot_i / tot_g - 1:+.0%})")

# Short-interval interferograms can carry a small systematic phase bias that adds up in an SBAS
# network. Test it: re-solve the SENG pixel without the shortest pairs and compare with GNSS.
section("SENG re-solved as one averaged pixel, with and without the shortest pairs")
obs_seng = np.array([np.nanmedian(u[sr - 1:sr + 2, sc - 1:sc + 2]) for u in unw])
for min_days in (0, 12, 18):
    keep = np.array([days(p) >= min_days for p in pairs]) & np.isfinite(obs_seng) & use_seng
    Ak = A[keep]
    rank = np.linalg.matrix_rank(Ak)
    xk = np.linalg.solve(Ak.T @ Ak + R, Ak.T @ obs_seng[keep])
    ik = np.r_[0, xk] * RAD2MM
    ik = ik - ik[k0]
    tab, *_, rms_k, r_k = by_phase(ik)
    print(f"pairs >= {min_days:2d} days: {keep.sum():3d} pairs, rank {rank}/{m}; summed inflation InSAR "
          f"{tab['InSAR_mm'].sum():+.0f} vs GNSS {tab['GNSS_mm'].sum():+.0f} mm "
          f"({tab['InSAR_mm'].sum() / tab['GNSS_mm'].sum() - 1:+.0%}); last phase {tab['InSAR_mm_per_yr'].iloc[-1]:+.0f} "
          f"vs {tab['GNSS_mm_per_yr'].iloc[-1]:+.0f} mm/yr; RMS {rms_k:.1f} mm")
pd.DataFrame({"date": dt, "insar_los_mm": i_ep, "gnss_los_mm": g_ep, "phase": seg, "insar_phase_mm": i_seg,
              "gnss_phase_mm": g_seg}).to_csv(
    OUT / "seng_insar_gnss_los.csv", index=False, float_format="%.2f")

fig, (ax, ax1, ax2) = plt.subplots(3, 1, figsize=(10, 7.6), sharex=True, height_ratios=[3, 3, 1.2],
                                   layout="constrained")
ax.plot(rel.index, rel["LOS"] - rel["LOS"].loc[dt[k0] - pd.Timedelta(days=1):dt[k0] + pd.Timedelta(days=1)].mean(),
        color="0.55", lw=0.8, label=f"GNSS SENG \u2212 {REF_STATION}, daily, projected into LOS")
ax.plot(dt, i_ep, "o", ms=3, color="#c0392b", label="InSAR SBAS at SENG (median of 3 x 3 pixels)")
ax.set_title(f"(a) Whole series: errors at eruptions accumulate (RMS difference {np.sqrt(np.mean(diff ** 2)):.0f} mm)")
for n in sorted(set(seg[seg >= 0])):
    k = np.nonzero(seg == n)[0]
    ax1.plot(dt[k], g_seg[k], "-", color="0.55", lw=1.6, label="GNSS" if n == seg[seg >= 0].min() else None)
    ax1.plot(dt[k], i_seg[k], "o", ms=3, color="#c0392b", label="InSAR" if n == seg[seg >= 0].min() else None)
ax1.set_title(f"(b) Each inflation phase restarted at zero: RMS difference {rms2:.0f} mm, r = {r2:.2f}")
for a in (ax, ax1):
    for e in events:
        a.axvline(pd.Timestamp(e), color="0.35", lw=0.8, ls=":")
    a.set_ylabel("LOS displacement (mm)")
    a.legend(loc="upper left", fontsize=8)
    a.grid(alpha=0.3)
ax2.bar(dt, i_seg - g_seg, width=5, color="#2c7fb8")
ax2.axhline(0, color="0.3", lw=0.8)
ax2.set_ylabel("InSAR \u2212 GNSS\nin (b), mm")
ax2.grid(alpha=0.3)
ax2.xaxis.set_major_formatter(mdates.DateFormatter("%b %Y"))
fig.savefig(OUT / "03-seng-insar-vs-gnss.png")
plt.close(fig)

# --- Step 7: coherence loss over new lava -----------------------------------------------------
section(f"coherence across the {LAVA_EVENT} eruption")
before = [p for p in pairs if p[9:] < LAVA_EVENT.replace("-", "")]
across = [p for p in pairs if p[:8] < LAVA_EVENT.replace("-", "") < p[9:]]
p_across = min(across, key=days)
recent = [p for p in before if days(p) == days(p_across) and (dt[idx[p_across[:8]]] - pd.Timestamp(p[:8])).days <= 36]
p_before = max(recent, key=lambda p: med_coh[pairs.index(p)])     # the most coherent recent pair of equal span
full = {}
for p in (p_before, p_across):
    f = DATA / f"cc_{p}.npz"
    if not f.exists():
        a, t = read_window(f"{BASE}/interferograms/{p}/{p}.geo.cc.tif")
        np.savez_compressed(f, a=a.astype("uint8"), t=np.array(t)[:6])
    full[p] = np.load(f)["a"] / 255
    full[p][full[p] == 0] = np.nan
zt = rasterio.Affine(*np.load(DATA / f"cc_{p_across}.npz")["t"])
fx = (zt.c + zt.a * (np.arange(full[p_across].shape[1]) + 0.5) - SENG[0]) * 111.32 * np.cos(np.radians(SENG[1]))
fy = (zt.f + zt.e * (np.arange(full[p_across].shape[0]) + 0.5) - SENG[1]) * 110.6
zoom = (slice(np.searchsorted(-fy, -4), np.searchsorted(-fy, 9)), slice(np.searchsorted(fx, -4), np.searchsorted(fx, 9)))
cb, ca = full[p_before][zoom], full[p_across][zoom]
lost = (cb > 0.5) & (ca < 0.2)
pix_km2 = zt.a * 111.32 * np.cos(np.radians(SENG[1])) * -zt.e * 110.6
print(f"before: {p_before} median coherence {np.nanmedian(cb):.2f}; across: {p_across} {np.nanmedian(ca):.2f} "
      f"(zoom 13 x 13 km around Sundhnukur)")
print(f"pixels coherent before (> 0.5) but incoherent across the eruption (< 0.2): {lost.sum() * pix_km2:.1f} km2")
fig, axes = plt.subplots(1, 3, figsize=(13, 4.2), sharey=True, layout="constrained")
zext = (fx[zoom[1]][0], fx[zoom[1]][-1], fy[zoom[0]][-1], fy[zoom[0]][0])
fmt = lambda p: f"{pd.Timestamp(p[:8]):%d %b}\u2013{pd.Timestamp(p[9:]):%d %b %Y}"
for ax, arr, title in [(axes[0], cb, f"(a) before: {fmt(p_before)}"), (axes[1], ca, f"(b) across: {fmt(p_across)}")]:
    im = ax.imshow(arr, cmap=cm.batlow, vmin=0, vmax=1, extent=zext, interpolation="nearest")
    ax.set_title(title)
fig.colorbar(im, ax=axes[:2], shrink=0.8, label="coherence")
im = axes[2].imshow(ca - cb, cmap=cm.vik_r, vmin=-0.8, vmax=0.8, extent=zext, interpolation="nearest")
axes[2].set_title("(c) change in coherence, (b) \u2212 (a)")
fig.colorbar(im, ax=axes[2], shrink=0.8, label="coherence change")
for ax in axes:
    ax.plot(0, 0, "w^", ms=7, mec="k")
    ax.set_xlabel("x, km east of SENG")
axes[0].set_ylabel("y, km north of SENG")
fig.savefig(OUT / "04-coherence-loss-lava.png")
plt.close(fig)

# --- Step 8: closure errors and usable pairs --------------------------------------------------
fig, (ax, ax2) = plt.subplots(1, 2, figsize=(12, 4), layout="constrained", width_ratios=[1.3, 1])
tmid = pd.to_datetime([pairs[k][:8] for k in tri[:, 0]])
ax.scatter(tmid, frac_tri * 100, s=8, c="#c0392b")
for e in events:
    ax.axvline(pd.Timestamp(e), color="0.35", lw=0.8, ls=":")
ax.set_ylabel("pixels with a 2\u03c0 closure error (%)")
ax.set_title(f"(a) Phase closure of {len(tri)} triplets (by first date)")
ax.xaxis.set_major_locator(mdates.MonthLocator((1, 7)))
ax.xaxis.set_major_formatter(mdates.DateFormatter("%b %Y"))
ax.grid(alpha=0.3)
im = ax2.imshow(np.where(land, frac_valid, np.nan), cmap=cm.batlow, vmin=0, vmax=1, extent=extent,
                interpolation="nearest")
ax2.plot(0, 0, "w^", ms=7, mec="k")
ax2.set_title(f"(b) share of the {n_ifg} pairs usable (coherence \u2265 {MIN_COH})")
ax2.set_xlabel("x, km east of SENG")
ax2.set_ylabel("y, km north of SENG")
fig.colorbar(im, ax=ax2, shrink=0.85, label="fraction of pairs")
fig.savefig(OUT / "05-closure-and-usable-pairs.png")
plt.close(fig)
print(f"land pixels usable in < 50% of pairs: {np.mean(frac_valid[land] < 0.5):.1%}; at SENG {frac_valid[sr, sc]:.0%}")
print("\nfigures and seng_insar_gnss_los.csv written to output/")