Source code: InSAR volcano deformation in Python: mapping magma inflation at Svartsengi, Iceland, with Sentinel-1

"""Volcano deformation with InSAR, part 1: mapping magma inflation at Svartsengi, Iceland.

Streams COMET LiCSAR Sentinel-1 interferograms (ascending frame 016A and descending frame 155D)
over the Reykjanes Peninsula, shows wrapped phase, unwrapped phase and coherence for one pair,
counts the fringes, adds consecutive pairs into about 48 days of line-of-sight (LOS) displacement
on each track, decomposes the two LOS maps into vertical and east-west motion using the LOS unit
vectors, and checks the result against the SENG GNSS station (Nevada Geodetic Laboratory).
It also writes the LOS maps and unit vectors to data/ and GeoTIFFs to output/ for part 2.
Run with:  python insar_volcano_deformation.py
"""
import io
import platform
import warnings
import re
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.patheffects as pe
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 rows 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"
TRACKS = {  # Sentinel-1 frame, and two consecutive 24-day pairs that span about 48 days
    "ascending": ("016A_02504_162118", ["20250815_20250908", "20250908_20251002"]),
    "descending": ("155D_02484_191814", ["20250813_20250906", "20250906_20250930"]),
}
BBOX = (-22.80, 63.72, -22.05, 64.05)    # lon/lat window: Reykjanes Peninsula around Svartsengi
SENG = (-22.4338, 63.8822)               # GNSS station SENG at the Svartsengi power station
REF_KM = 15                              # reference = median of land pixels > 15 km from SENG
MIN_COH = 0.3                            # ignore pixels with coherence below this in any pair
WAVELENGTH = 0.0554658                   # Sentinel-1 C-band wavelength in metres
RAD2MM = -WAVELENGTH / (4 * np.pi) * 1000  # LiCSAR convention: positive = towards the satellite
SHOWCASE = "20250908_20260605"           # ascending 270-day pair used to explain fringes
PROFILE_LAT, PROFILE_LON = 63.868, (-22.75, -22.40)  # from the far field in the west to the LOS peak
NGL = "https://geodesy.unr.edu/gps_timeseries/IGS20/tenv3/IGS20"

rasterio.Env(GDAL_DISABLE_READDIR_ON_OPEN="EMPTY_DIR", GDAL_HTTP_MAX_RETRY=4, GDAL_HTTP_RETRY_DELAY=2).__enter__()


def base(frame: str) -> str:
    return f"{LICSAR}/{int(frame[:3])}/{frame}"


def read_window(url: str):
    """Read only the rows of a remote LiCSAR GeoTIFF that cover BBOX (HTTP range requests)."""
    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)


def cached(name: str, url: str):
    """Read a window once and keep it in data/ so that reruns are instant."""
    f = DATA / f"{name}.npz"
    if not f.exists():
        a, t = read_window(url)
        np.savez_compressed(f, a=a, t=np.array(t)[:6])
    z = np.load(f)
    return z["a"], rasterio.Affine(*z["t"])


def deramp(u: np.ndarray, use: np.ndarray, LON: np.ndarray, LAT: np.ndarray):
    """Fit a + b*lon + c*lat to the far-field pixels and subtract it (orbit/ionosphere/long-wave tropo)."""
    A = np.column_stack([np.ones(use.sum()), LON[use] - LON.mean(), LAT[use] - LAT.mean()])
    coef = np.linalg.lstsq(A, u[use], rcond=None)[0]
    return u - (coef[0] + coef[1] * (LON - LON.mean()) + coef[2] * (LAT - LAT.mean())), coef


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


# --- Step 1: the two Sentinel-1 frames ---------------------------------------------------------
section("LiCSAR frames over Svartsengi")
meta = {}
for track, (frame, pairs) in TRACKS.items():
    html = requests.get(f"{base(frame)}/interferograms/", timeout=60).text
    online = sorted(set(re.findall(r'href="(\d{8}_\d{8})/"', html)))
    m = dict(tok.split("=", 1) for tok in requests.get(f"{base(frame)}/metadata/metadata.txt", timeout=60)
             .text.split() if "=" in tok)
    meta[track] = m
    print(f"{track:10s} {frame}: {len(online):,} interferograms online ({online[0][:8]} to {online[-1][9:]}), "
          f"heading {float(m['heading']):.1f} deg, mean incidence {float(m['avg_incidence_angle']):.1f} deg, "
          f"overpass {':'.join(m['center_time'].split(':')[:2]):>5s} UTC")
    assert all(p in online for p in pairs), f"missing pair on {frame}"
    print(f"           using {' + '.join(pairs)} = {sum(days(p) for p in pairs)} days")

# --- Step 2: one interferogram: wrapped phase, unwrapped phase and coherence ---------------------
frame = TRACKS["ascending"][0]
url = f"{base(frame)}/interferograms/{SHOWCASE}/{SHOWCASE}"
wrapped, transform = cached(f"asc_{SHOWCASE}_pha", f"{url}.geo.diff_pha.tif")
unw, _ = cached(f"asc_{SHOWCASE}_unw", f"{url}.geo.unw.tif")
coh, _ = cached(f"asc_{SHOWCASE}_cc", f"{url}.geo.cc.tif")
hgt, _ = cached("asc_hgt", f"{base(frame)}/metadata/{frame}.geo.hgt.tif")
H, W = unw.shape
land = unw != 0                                          # 0 = sea or not unwrapped
unw[~land], wrapped[~land] = np.nan, np.nan
coh = np.where(land, coh / 255, np.nan)                  # LiCSAR stores coherence as 0-255
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)
dist_km = np.hypot((LON - SENG[0]) * 111.32 * np.cos(np.radians(SENG[1])), (LAT - SENG[1]) * 110.6)
far = land & (dist_km > REF_KM)
inv = ~transform
sc, sr = (int(v) for v in inv * SENG)
section("showcase interferogram")
print(f"pair {SHOWCASE} ({days(SHOWCASE)} days), window {H} x {W} pixels of {transform.a:.3f} deg "
      f"(about {transform.a * 111_320 * np.cos(np.radians(63.9)):.0f} m x {-transform.e * 110_600:.0f} m)")
print(f"land pixels {land.mean():.0%}; median coherence on land {np.nanmedian(coh):.2f}; "
      f"far-field reference pixels (> {REF_KM} km from SENG): {far.sum():,}")
unw_ref = unw - np.nanmedian(unw[far])
at_seng = np.nanmedian(unw_ref[sr - 1:sr + 2, sc - 1:sc + 2])
print(f"unwrapped phase at SENG relative to the far field: {at_seng:.2f} rad = {at_seng / (2 * np.pi):.2f} "
      f"fringes = {at_seng * RAD2MM:+.1f} mm LOS (one fringe = {abs(2 * np.pi * RAD2MM):.1f} mm)")

# A west-east profile from the far field to the centre of the fringes: count the 2*pi jumps
pr = int((~transform * (0, PROFILE_LAT))[1])
row_w = np.angle(np.nanmean(np.exp(1j * wrapped[pr - 3:pr + 4]), axis=0))  # circular mean of 7 rows
row_u = np.nanmedian(unw_ref[pr - 3:pr + 4], axis=0)
ok = (lon > PROFILE_LON[0]) & (lon < PROFILE_LON[1]) & np.isfinite(row_w) & np.isfinite(row_u)
ok &= lon <= lon[ok][np.argmax(np.abs(row_u[ok]))]      # stop at the peak of the LOS signal
jumps = np.abs(np.diff(row_w[ok])) > np.pi
cycles = (row_u[ok][-1] - row_u[ok][0]) / (2 * np.pi)
print(f"profile along {PROFILE_LAT} N from {lon[ok][0]:.3f} to the LOS peak at {lon[ok][-1]:.3f} E: "
      f"{jumps.sum()} wrapped-phase jumps counted; unwrapped phase changes by {cycles:.2f} cycles "
      f"= {abs(cycles * 2 * np.pi * RAD2MM):.0f} mm")

extent = (transform.c, transform.c + transform.a * W, transform.f + transform.e * H, transform.f)
gy, gx = np.gradient(np.where(np.isfinite(hgt), hgt, 0), 110.6 * -transform.e * 1000,
                     111_320 * np.cos(np.radians(63.9)) * transform.a)
hill = np.where(land, -(gx - gy) / np.sqrt(1 + gx ** 2 + gy ** 2), np.nan)  # simple NW-lit hillshade
asp = 1 / np.cos(np.radians(63.9))


def mark(ax, label=True):
    ax.plot(*SENG, "^", ms=6, mfc="white", mec="black", mew=1)
    if label:
        ax.annotate("SENG", SENG, xytext=(5, 5), textcoords="offset points", fontsize=8,
                    path_effects=[pe.withStroke(linewidth=2.5, foreground="white")])
    ax.set_aspect(asp)
    ax.set_facecolor("#dfe7ee")                          # sea / no data


fig, axes = plt.subplots(1, 3, figsize=(13, 3.9), sharey=True)
for ax, arr, cmap, vmin, vmax, lab in [
        (axes[0], wrapped, cm.romaO, -np.pi, np.pi, "wrapped phase (rad)"),
        (axes[1], unw_ref * RAD2MM, cm.vik, -200, 200, "unwrapped LOS displacement (mm)"),
        (axes[2], coh, cm.batlow, 0, 1, "coherence")]:
    im = ax.imshow(arr, cmap=cmap, vmin=vmin, vmax=vmax, extent=extent, interpolation="nearest")
    fig.colorbar(im, ax=ax, shrink=0.8, label=lab, orientation="horizontal", pad=0.17)
    mark(ax)
axes[0].set_title(f"(a) wrapped interferogram {SHOWCASE}")
axes[1].set_title("(b) unwrapped (red = towards satellite)")
axes[2].set_title("(c) coherence")
axes[0].set_ylabel("latitude")
for ax in axes:
    ax.set_xlabel("longitude")
fig.savefig(OUT / "01-wrapped-unwrapped-coherence.png")
plt.close(fig)

fig, ax = plt.subplots(figsize=(8, 3.6))
ax.plot(lon[ok], row_w[ok], ".", ms=2.5, color="#7a5195", label="wrapped phase")
ax.plot(lon[ok], row_u[ok], "-", lw=2, color="#ef5675", label="unwrapped phase (far field = 0)")
for x in lon[ok][1:][jumps]:
    ax.axvline(x, color="0.6", lw=0.6, ls=":")
ax.annotate(f"{abs(cycles):.1f} cycles\n= {abs(cycles * 2 * np.pi * RAD2MM):.0f} mm LOS", (lon[ok][-1], row_u[ok][-1]),
            xytext=(-60, 25), textcoords="offset points", fontsize=8, arrowprops=dict(arrowstyle="->"))
ax.set_xlabel("longitude along latitude %.3f N" % PROFILE_LAT)
ax.set_ylabel("phase (rad)")
sec = ax.secondary_yaxis("right", functions=(lambda r: r * RAD2MM, lambda mm: mm / RAD2MM))
sec.set_ylabel("LOS displacement (mm)")
ax.set_title(f"Counting fringes: {jumps.sum()} jumps of 2\u03c0 from the far field to the peak ({SHOWCASE})")
ax.legend(loc="lower left", fontsize=8)
ax.grid(alpha=0.3)
fig.savefig(OUT / "02-fringe-profile.png")
plt.close(fig)

# --- Step 3: add consecutive pairs into ~48 days of LOS displacement per track -----------------
section("cumulative LOS displacement (~48 days)")
los, vec, dates = {}, {}, {}
for track, (frame, pairs) in TRACKS.items():
    tag = "asc" if track == "ascending" else "desc"
    total, good = np.zeros((H, W), "float32"), np.ones((H, W), bool)
    for p in pairs:
        u, t = cached(f"{tag}_{p}_unw", f"{base(frame)}/interferograms/{p}/{p}.geo.unw.tif")
        c, _ = cached(f"{tag}_{p}_cc", f"{base(frame)}/interferograms/{p}/{p}.geo.cc.tif")
        assert u.shape == (H, W) and np.allclose(np.array(t)[:6], np.array(transform)[:6]), "grids differ"
        ok_p = (u != 0) & (c / 255 >= MIN_COH) & land
        good &= ok_p
        u = np.where(ok_p, u, np.nan)
        u, coef = deramp(u, far & ok_p, LON, LAT)        # re-reference AND remove a plane, every pair
        total += np.nan_to_num(u) * RAD2MM
        print(f"{track:10s} {p}: median coherence {np.median(c[land]) / 255:.2f}, "
              f"usable {ok_p.sum() / land.sum():.0%} of land, ramp removed {abs(coef[1] * RAD2MM) * 0.75:.0f} mm "
              f"W-E and {abs(coef[2] * RAD2MM) * 0.33:.0f} mm S-N across the window")
    los[track] = np.where(good, total, np.nan)
    vec[track] = np.stack([cached(f"{tag}_{k}", f"{base(frame)}/metadata/{frame}.geo.{k}.tif")[0]
                           for k in "ENU"])
    dates[track] = (pairs[0][:8], pairs[-1][9:])
    s = np.nanmedian(los[track][sr - 1:sr + 2, sc - 1:sc + 2])
    print(f"{track:10s} {dates[track][0]}-{dates[track][1]}: max {np.nanmax(los[track]):+.1f} mm, "
          f"SENG {s:+.1f} mm, far-field scatter (std) {np.nanstd(los[track][far]):.1f} mm, "
          f"usable {np.isfinite(los[track]).sum() / land.sum():.0%} of land")

section("LOS unit vectors at SENG (ground -> satellite)")
for track in TRACKS:
    e, n, u = vec[track][:, sr, sc]
    print(f"{track:10s} E {e:+.3f}  N {n:+.3f}  U {u:+.3f}  -> incidence {np.degrees(np.arccos(u)):.1f} deg, "
          f"look azimuth {np.degrees(np.arctan2(-e, -n)) % 360:.1f} deg")

fig, axes = plt.subplots(1, 2, figsize=(11, 4.2), sharey=True)
for ax, track in zip(axes, TRACKS):
    ax.imshow(hill, cmap="Greys", vmin=-1.2, vmax=0.6, extent=extent, interpolation="bilinear")
    im = ax.imshow(los[track], cmap=cm.vik, vmin=-60, vmax=60, extent=extent,
                   interpolation="nearest", alpha=0.85)
    mark(ax)
    e, n, _ = vec[track][:, sr, sc]
    ax.annotate("", xy=(-22.12 - 0.06 * e / np.hypot(e, n), 63.76 - 0.03 * n / np.hypot(e, n)),
                xytext=(-22.12, 63.76), arrowprops=dict(arrowstyle="-|>", color="black", lw=1.5))
    ax.text(-22.12, 63.775, "look", ha="center", fontsize=7)
    ax.set_title(f"{track} ({TRACKS[track][0][:4]}) {dates[track][0]}\u2013{dates[track][1]}")
    ax.set_xlabel("longitude")
axes[0].set_ylabel("latitude")
fig.colorbar(im, ax=axes, shrink=0.85, label="LOS displacement (mm, + = towards satellite)")
fig.savefig(OUT / "03-los-ascending-descending.png")
plt.close(fig)

# --- Step 4: decompose two LOS maps into east and up -------------------------------------------
# d_asc = E_a*e + N_a*n + U_a*u and d_desc = E_d*e + N_d*n + U_d*u. Two equations, three unknowns:
# north is dropped because both near-polar orbits are almost blind to it (N is about +/-0.1).
section("asc + desc decomposition")
Ea, Na, Ua = vec["ascending"]
Ed, Nd, Ud = vec["descending"]
det = Ea * Ud - Ed * Ua
d_a, d_d = los["ascending"], los["descending"]
east = (Ud * d_a - Ua * d_d) / det
up = (Ea * d_d - Ed * d_a) / det
G = np.array([[Ea[sr, sc], Ua[sr, sc]], [Ed[sr, sc], Ud[sr, sc]]])
cov = np.linalg.inv(G.T @ G)
print(f"design matrix at SENG [[E_a U_a] [E_d U_d]] = {[[round(float(v), 3) for v in row] for row in G]}, condition number "
      f"{np.linalg.cond(G):.2f}")
print(f"noise amplification for 1 mm LOS noise: east {np.sqrt(cov[0, 0]):.2f} mm, up {np.sqrt(cov[1, 1]):.2f} mm")
iy, ix = np.unravel_index(np.nanargmax(up), up.shape)
print(f"maximum uplift {np.nanmax(up):.1f} mm at {lon[ix]:.3f} E, {lat[iy]:.3f} N; "
      f"east-west range {np.nanpercentile(east, 1):+.1f} to {np.nanpercentile(east, 99):+.1f} mm (1st-99th pct)")
up_s = np.nanmedian(up[sr - 1:sr + 2, sc - 1:sc + 2])
east_s = np.nanmedian(east[sr - 1:sr + 2, sc - 1:sc + 2])
print(f"at SENG: up {up_s:+.1f} mm, east {east_s:+.1f} mm")
px_km2 = transform.a * 111.32 * np.cos(np.radians(63.9)) * -transform.e * 110.6
print(f"area with > 20 mm uplift: {(up > 20).sum() * px_km2:.0f} km2")

fig, axes = plt.subplots(1, 2, figsize=(11, 4.2), sharey=True)
for ax, arr, lim, title in [(axes[0], up, 80, "vertical (+ = up)"),
                            (axes[1], east, 30, "east\u2013west (+ = east)")]:
    ax.imshow(hill, cmap="Greys", vmin=-1.2, vmax=0.6, extent=extent, interpolation="bilinear")
    im = ax.imshow(arr, cmap=cm.vik, vmin=-lim, vmax=lim, extent=extent, interpolation="nearest",
                   alpha=0.85)
    if ax is axes[0]:
        ax.contour(lon, lat, np.where(np.isfinite(up), up, 0), levels=[20, 40, 60], colors="black",
                   linewidths=0.6)
    fig.colorbar(im, ax=ax, shrink=0.8, label="displacement (mm)", orientation="horizontal", pad=0.17)
    mark(ax)
    ax.set_title(f"{title}, Aug\u2013Oct 2025")
    ax.set_xlabel("longitude")
axes[0].set_ylabel("latitude")
fig.savefig(OUT / "04-vertical-east-decomposition.png")
plt.close(fig)

# --- Step 5: compare with the SENG GNSS station ------------------------------------------------
section("GNSS check: SENG (Nevada Geodetic Laboratory, IGS20)")
txt = requests.get(f"{NGL}/SENG.tenv3", timeout=60).text
g = pd.read_csv(io.StringIO(txt), 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")
g["E"] = (g["e0"] + g["east"]) * 1000                    # mm; e0/n0/u0 change when the reference does
g["N"] = (g["n0"] + g["north"]) * 1000
g["U"] = (g["u0"] + g["up"] + g["ant"]) * 1000           # add the antenna height back (0 / 0.5 m steps)
print(f"{len(g):,} daily solutions {g.index[0]:%Y-%m-%d} to {g.index[-1]:%Y-%m-%d}; "
      f"median sigma E/N/U {g['sig_e'].median() * 1000:.1f}/{g['sig_n'].median() * 1000:.1f}/"
      f"{g['sig_u'].median() * 1000:.1f} mm")


def gnss_at(day: str) -> pd.Series:
    t = pd.Timestamp(day)
    return g.loc[t - pd.Timedelta(days=1):t + pd.Timedelta(days=1), ["E", "N", "U"]].mean()


rows = []
for track in TRACKS:
    d0, d1 = dates[track]
    dg = gnss_at(d1) - gnss_at(d0)
    e, n, u = vec[track][:, sr, sc]
    g_los = e * dg["E"] + n * dg["N"] + u * dg["U"]
    s = np.nanmedian(los[track][sr - 1:sr + 2, sc - 1:sc + 2])
    rows.append((f"{track[:3] if track == 'ascending' else track[:4]}. LOS", float(s), float(g_los)))
    print(f"{track:10s} {d0}-{d1}: GNSS dE {dg['E']:+.1f} dN {dg['N']:+.1f} dU {dg['U']:+.1f} mm -> "
          f"LOS {g_los:+.1f} mm; InSAR {s:+.1f} mm; difference {s - g_los:+.1f} mm")
dg = gnss_at("20251001") - gnss_at("20250814")
rows += [("vertical", float(up_s), float(dg["U"])), ("east", float(east_s), float(dg["E"]))]
dg_sc = gnss_at(SHOWCASE[9:]) - gnss_at(SHOWCASE[:8])
e, n, u = vec["ascending"][:, sr, sc]
g_sc = e * dg_sc["E"] + n * dg_sc["N"] + u * dg_sc["U"]
print(f"showcase pair {SHOWCASE}: GNSS predicts {g_sc:+.1f} mm LOS = {g_sc / abs(2 * np.pi * RAD2MM):.2f} fringes "
      f"at SENG; the unwrapped interferogram gives {at_seng * RAD2MM:+.1f} mm = {abs(at_seng) / (2 * np.pi):.2f}")
print(f"decomposed vs GNSS (14 Aug-1 Oct 2025): up {up_s:+.1f} vs {dg['U']:+.1f} mm, "
      f"east {east_s:+.1f} vs {dg['E']:+.1f} mm")
cmp_df = pd.DataFrame(rows, columns=["component", "InSAR_mm", "GNSS_mm"]).set_index("component")
print(cmp_df.astype("float64").round(1).to_string())

fig, (ax, bx) = plt.subplots(1, 2, figsize=(12, 3.8), gridspec_kw={"width_ratios": [2.3, 1]})
s = g.loc["2023-06-01":, "U"]
ax.plot(s.index, s - s.iloc[0], ".", ms=1.5, color="#003f5c")
ax.axvspan(pd.Timestamp("2025-08-14"), pd.Timestamp("2025-10-01"), color="#ffa600", alpha=0.35,
           label="InSAR window")
ax.set_ylabel("SENG vertical (mm, from Jun 2023)")
ax.xaxis.set_major_locator(mdates.MonthLocator((1, 7)))
ax.xaxis.set_major_formatter(mdates.DateFormatter("%b\n%Y"))
ax.set_title("SENG GNSS: inflation, then sudden drops at dyke intrusions and eruptions")
ax.legend(loc="lower left", fontsize=8)
ax.grid(alpha=0.3)
x = np.arange(len(cmp_df))
bx.bar(x - 0.2, cmp_df["InSAR_mm"], 0.4, color="#ef5675", label="InSAR")
bx.bar(x + 0.2, cmp_df["GNSS_mm"], 0.4, color="#003f5c", label="GNSS")
bx.set_xticks(x, cmp_df.index, rotation=20)
bx.axhline(0, color="black", lw=0.6)
bx.set_ylabel("displacement at SENG (mm)")
bx.set_title("Aug\u2013Oct 2025 at SENG")
bx.legend(fontsize=8)
bx.grid(alpha=0.3, axis="y")
fig.savefig(OUT / "05-gnss-comparison.png")
plt.close(fig)

# --- Step 6: save the products for part 2 (source modelling) ------------------------------------
section("outputs")
np.savez_compressed(DATA / "svartsengi_los_2025.npz", lon=lon, lat=lat, asc=los["ascending"],
                    desc=los["descending"], vec_asc=vec["ascending"], vec_desc=vec["descending"],
                    up=up, east=east, transform=np.array(transform)[:6])
profile = dict(driver="GTiff", height=H, width=W, count=1, dtype="float32", crs="EPSG:4326",
               transform=transform, nodata=np.nan, compress="deflate")
for name, arr in [("los_ascending_mm", los["ascending"]), ("los_descending_mm", los["descending"]),
                  ("vertical_mm", up), ("east_mm", east)]:
    with rasterio.open(OUT / f"svartsengi_{name}.tif", "w", **profile) as dst:
        dst.write(arr.astype("float32"), 1)
print("wrote data/svartsengi_los_2025.npz and output/svartsengi_{los_ascending,los_descending,vertical,east}_mm.tif")