"""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")