"""InSAR ground subsidence in Python: one year of Sentinel-1 interferograms over Mexico City. Lists the COMET LiCSAR interferograms for Sentinel-1 frame 143D_07132_131313, streams a 53 x 44 km window of the unwrapped phase and coherence for every short-baseline pair from September 2025 to September 2026, converts phase to line-of-sight (LOS) displacement in millimetres, masks low coherence, references everything to a stable point, inverts the network into a displacement time series (small-baseline least squares, the core of MintPy and LiCSBAS), fits a velocity, converts it to vertical subsidence and draws a map, a profile and time series. It also writes the stack in MintPy's HDF5 layout (ifgramStack.h5) for mintpy_crosscheck.py. Run with: python insar_subsidence.py """ import platform import re from datetime import datetime from importlib.metadata import version from pathlib import Path import h5py 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 from matplotlib.colors import TwoSlopeNorm import requests from rasterio.windows import from_bounds from scipy import ndimage 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", 120) print("=== versions ===") print(f"python {platform.python_version()} | rasterio {rasterio.__version__} | numpy {np.__version__} | " f"scipy {version('scipy')} | h5py {h5py.__version__} | matplotlib {matplotlib.__version__}") def section(title: str) -> None: print(f"\n=== {title} ===") FRAME = "143D_07132_131313" # Sentinel-1 relative orbit 143, descending BASE = f"https://gws-access.jasmin.ac.uk/public/nceo_geohazards/LiCSAR_products.public/143/{FRAME}" BBOX = (-99.30, 19.20, -98.80, 19.60) # lon/lat window: Mexico City, Nezahualcoyotl and Chalco START, END = "20250911", "20260913" # one year of 6-day Sentinel-1 acquisitions MAX_DAYS = 18 # keep only short temporal baselines: 6, 12 and 18 days WAVELENGTH = 0.0554658 # Sentinel-1 C-band wavelength in metres RAD2MM = -WAVELENGTH / (4 * np.pi) * 1000 # LiCSAR convention: positive = towards the satellite REF = (-99.185, 19.325) # stable reference: Pedregal lava field, Ciudad Universitaria SHOWCASE = "20260103_20260304" # a 60-day dry-season pair to show wrapped vs unwrapped phase POINTS = {"Chalco": (-98.900, 19.262), "Tláhuac": (-99.004, 19.286), "Nezahualcóyotl": (-99.010, 19.410), "Airport (AICM)": (-99.072, 19.436), "Centro Histórico": (-99.133, 19.433), "Reference (CU)": REF} PROFILE = ((-99.20, 19.43), (-98.88, 19.25)) # Chapultepec to Chalco, through the old lakebed rasterio.Env(GDAL_DISABLE_READDIR_ON_OPEN="EMPTY_DIR", GDAL_HTTP_MAX_RETRY=4, GDAL_HTTP_RETRY_DELAY=2).__enter__() 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), src.window_transform(win) def days(pair: str) -> int: return (datetime.strptime(pair[9:], "%Y%m%d") - datetime.strptime(pair[:8], "%Y%m%d")).days # --- Step 1: list the interferograms --------------------------------------------------- html = requests.get(f"{BASE}/interferograms/", timeout=60).text all_pairs = sorted(set(re.findall(r'href="(\d{8}_\d{8})/"', html))) pairs = [p for p in all_pairs if p[:8] >= START and p[9:] <= END and days(p) <= MAX_DAYS] epochs = sorted({p[:8] for p in pairs} | {p[9:] for p in pairs}) section("LiCSAR interferograms") print(f"frame {FRAME}: {len(all_pairs):,} interferograms online, {all_pairs[0][:8]} to {all_pairs[-1][9:]}") print(f"selected: {len(pairs)} pairs of <= {MAX_DAYS} days between {len(epochs)} epochs, {epochs[0]} to {epochs[-1]}") print("temporal baselines (days):", pd.Series([days(p) for p in pairs]).value_counts().sort_index().to_dict()) meta = dict(line.split("=", 1) for line in requests.get(f"{BASE}/metadata/metadata.txt", timeout=60).text.split() if "=" in line) print(f"heading {float(meta['heading']):.1f}°, mean incidence {float(meta['avg_incidence_angle']):.1f}°") # --- Step 2: stream the window of every interferogram --------------------------------------- cache = DATA / f"stack_{START}_{END}_{MAX_DAYS}d.npz" if cache.exists(): z = np.load(cache) unw, coh, transform = z["unw"], z["coh"], rasterio.Affine(*z["transform"]) else: unw, coh = [], [] for i, p in enumerate(pairs): u, transform = read_window(f"{BASE}/interferograms/{p}/{p}.geo.unw.tif") c, _ = read_window(f"{BASE}/interferograms/{p}/{p}.geo.cc.tif") unw.append(u) coh.append(c) if i % 40 == 0: print(f" read {i + 1}/{len(pairs)}: {p}") unw, coh = np.array(unw, "float32"), np.array(coh, "uint8") np.savez_compressed(cache, unw=unw, coh=coh, transform=np.array(transform)[:6]) coh = coh.astype("float32") / 255 # LiCSAR stores coherence as 0-255 unw[unw == 0] = np.nan # 0 = not unwrapped (masked by LiCSAR) up = read_window(f"{BASE}/metadata/{FRAME}.geo.U.tif")[0] # up component of the LOS unit vector hgt = read_window(f"{BASE}/metadata/{FRAME}.geo.hgt.tif")[0].astype("float32") n_ifg, H, W = unw.shape PIXEL_KM2 = transform.a * 111.32 * np.cos(np.radians(19.4)) * -transform.e * 110.6 section("stack") print(f"{n_ifg} interferograms x {H} x {W} pixels ({transform.a * 3600:.1f} arcsec, about " f"{transform.a * 111_320 * np.cos(np.radians(19.4)):.0f} m x {-transform.e * 110_600:.0f} m)") print(f"unwrapped in every pair: {np.isfinite(unw).all(0).mean():.1%} of pixels; LOS up component " f"{np.nanmin(up):.3f}..{np.nanmax(up):.3f} (local incidence {np.degrees(np.arccos(np.nanmean(up))):.1f}°)") # --- Step 3: one interferogram: wrapped vs unwrapped ------------------------------------------- inv = ~transform rc, rr = (int(v) for v in inv * REF) ref_win = (slice(rr - 2, rr + 3), slice(rc - 2, rc + 3)) # 5 x 5 pixel reference area wrapped, _ = read_window(f"{BASE}/interferograms/{SHOWCASE}/{SHOWCASE}.geo.diff_pha.tif") show_unw, _ = read_window(f"{BASE}/interferograms/{SHOWCASE}/{SHOWCASE}.geo.unw.tif") show_coh, _ = read_window(f"{BASE}/interferograms/{SHOWCASE}/{SHOWCASE}.geo.cc.tif") show_unw[show_unw == 0] = np.nan show_coh = show_coh / 255 show_los = (show_unw - np.nanmedian(show_unw[ref_win])) * RAD2MM section(f"single pair {SHOWCASE} ({days(SHOWCASE)} days)") print(f"wrapped phase range {np.nanmin(wrapped):.2f}..{np.nanmax(wrapped):.2f} rad; " f"unwrapped {np.nanmin(show_unw):.1f}..{np.nanmax(show_unw):.1f} rad") print(f"one fringe (2π) = {abs(2 * np.pi * RAD2MM):.1f} mm of LOS motion") print(f"LOS change relative to CU: 1st percentile {np.nanpercentile(show_los, 1):.0f} mm, " f"99th percentile {np.nanpercentile(show_los, 99):.0f} mm") print(f"median coherence {np.nanmedian(show_coh):.2f}; pixels with coherence < 0.3: {(show_coh < 0.3).mean():.1%}") # --- Step 4: coherence mask and reference point ---------------------------------------------- mean_coh = coh.mean(0) mask = (mean_coh >= 0.35) & np.isfinite(unw).all(0) mask = ndimage.binary_opening(mask, iterations=1) section("coherence mask") print(f"mean coherence: median {np.median(mean_coh):.2f}; reference area {mean_coh[ref_win].mean():.2f}") print(f"kept {mask.sum():,} of {mask.size:,} pixels ({mask.mean():.1%}) with mean coherence >= 0.35 " f"and unwrapped in all {n_ifg} pairs") unw_ref = unw - np.nanmedian(unw[:, ref_win[0], ref_win[1]].reshape(n_ifg, -1), axis=1)[:, None, None] # --- Step 5: small-baseline (SBAS) inversion ----------------------------------------------- dates = pd.to_datetime(epochs) t_years = np.array((dates - dates[0]).days) / 365.25 idx = {e: i for i, e in enumerate(epochs)} A = np.zeros((n_ifg, len(epochs) - 1)) # unknowns: phase at each epoch after the first for k, p in enumerate(pairs): i, j = idx[p[:8]], idx[p[9:]] A[k, j - 1] = 1 if i > 0: A[k, i - 1] = -1 section("SBAS inversion") print(f"design matrix {A.shape}, rank {np.linalg.matrix_rank(A)} (needs {len(epochs) - 1}: one connected network)") obs = unw_ref[:, mask] # n_ifg x n_pixels x, *_ = np.linalg.lstsq(A, obs, rcond=None) resid = obs - A @ x ts = np.full((len(epochs), H, W), np.nan, "float32") ts[0, mask] = 0 ts[1:, mask] = x * RAD2MM # cumulative LOS displacement in mm rms_ifg = np.sqrt(np.mean(resid ** 2, axis=1)) * abs(RAD2MM) print(f"residual RMS per interferogram: median {np.median(rms_ifg):.2f} mm, max {rms_ifg.max():.2f} mm " f"({pairs[int(rms_ifg.argmax())]})") print(f"pixels with a residual > π in any pair (likely unwrapping errors): {(np.abs(resid) > np.pi).any(0).mean():.2%}") # --- Step 6: velocity, vertical conversion and outputs ------------------------------------ G = np.column_stack([np.ones_like(t_years), t_years]) coef, *_ = np.linalg.lstsq(G, ts[:, mask], rcond=None) v_los = np.full((H, W), np.nan, "float32") v_los[mask] = coef[1] fit_resid = ts[:, mask] - G @ coef v_std = np.full((H, W), np.nan, "float32") v_std[mask] = np.sqrt((fit_resid ** 2).sum(0) / (len(epochs) - 2) / ((t_years - t_years.mean()) ** 2).sum()) v_up = v_los / up # assumes the motion is purely vertical section("velocity") print(f"LOS velocity: 1st percentile {np.nanpercentile(v_los, 1):.0f} mm/yr, median {np.nanmedian(v_los):.0f} mm/yr, " f"99th percentile {np.nanpercentile(v_los, 99):.0f} mm/yr") print(f"vertical (LOS / up component): fastest {np.nanmin(v_up):.0f} mm/yr; " f"area subsiding faster than 100 mm/yr: {(v_up < -100).sum() * PIXEL_KM2:.0f} km², " f"faster than 200 mm/yr: {(v_up < -200).sum() * PIXEL_KM2:.0f} km²") print(f"median 1-sigma velocity uncertainty (scatter about the line): {np.nanmedian(v_std):.1f} mm/yr") def window_median(arr: np.ndarray, lon: float, lat: float, half: int = 1): c, r = (int(v) for v in inv * (lon, lat)) return np.nanmedian(arr[..., r - half:r + half + 1, c - half:c + half + 1].reshape(*arr.shape[:-2], -1), -1) table = pd.DataFrame([{"place": name, "height_m": round(float(window_median(hgt, *xy))), "coherence": round(float(window_median(mean_coh, *xy)), 2), "los_mm_yr": round(float(window_median(v_los, *xy)), 1), "vertical_mm_yr": round(float(window_median(v_up, *xy)), 1), "los_total_mm": round(float(window_median(ts[-1], *xy)), 1)} for name, xy in POINTS.items()]) print(table.to_string(index=False)) series = pd.DataFrame({name: window_median(ts, *xy) for name, xy in POINTS.items()}, index=dates) series.round(1).to_csv(OUT / "los_timeseries_points.csv", index_label="date") with rasterio.open(OUT / "vertical_velocity_mm_yr.tif", "w", driver="GTiff", height=H, width=W, count=1, dtype="float32", crs="EPSG:4326", transform=transform, nodata=np.nan, compress="deflate") as dst: dst.write(v_up, 1) # The same referenced stack in MintPy's ifgramStack.h5 layout, for mintpy_crosscheck.py with h5py.File(DATA / "ifgramStack.h5", "w") as f: f.create_dataset("date", data=np.array([[p[:8].encode(), p[9:].encode()] for p in pairs])) f.create_dataset("dropIfgram", data=np.ones(n_ifg, bool)) f.create_dataset("bperp", data=np.zeros(n_ifg, "float32")) f.create_dataset("unwrapPhase", data=np.nan_to_num(unw_ref).astype("float32")) f.create_dataset("coherence", data=coh) f.create_dataset("connectComponent", data=np.isfinite(unw).astype("int16")) f.attrs.update({"FILE_TYPE": "ifgramStack", "LENGTH": H, "WIDTH": W, "X_FIRST": transform.c, "Y_FIRST": transform.f, "X_STEP": transform.a, "Y_STEP": transform.e, "X_UNIT": "degrees", "Y_UNIT": "degrees", "WAVELENGTH": WAVELENGTH, "REF_Y": rr, "REF_X": rc, "UNIT": "radian", "REF_LAT": REF[1], "REF_LON": REF[0], "PLATFORM": "Sen", "ORBIT_DIRECTION": "DESCENDING", "HEADING": float(meta["heading"]), "INCIDENCE_ANGLE": float(meta["avg_incidence_angle"]), "STARTING_RANGE": float(meta["centre_range_m"]), "CENTER_LINE_UTC": 44762, "ALOOKS": 1, "RLOOKS": 1, "RANGE_PIXEL_SIZE": 100.0, "AZIMUTH_PIXEL_SIZE": 100.0}) np.save(DATA / "mask.npy", mask) np.save(DATA / "v_los_numpy.npy", v_los) # --- Step 7: figures ----------------------------------------------------------------------------- extent = [transform.c, transform.c + W * transform.a, transform.f + H * transform.e, transform.f] halo = [pe.withStroke(linewidth=2.2, foreground="white")] hill = ndimage.gaussian_filter(np.gradient(hgt, axis=1) - np.gradient(hgt, axis=0), 1) def label_points(ax, names=None): for name, (lon, lat) in POINTS.items(): if names and name not in names: continue ax.plot(lon, lat, "o", ms=4, mfc="white", mec="black", mew=0.8) ax.annotate(name, (lon, lat), xytext=(5, 4), textcoords="offset points", fontsize=7.5, path_effects=halo) def geo_axes(ax): ax.set_xlabel("longitude (°)") ax.set_ylabel("latitude (°)") ax.set_aspect(1 / np.cos(np.radians(19.4))) pretty = f"{datetime.strptime(SHOWCASE[:8], '%Y%m%d'):%d %b} – {datetime.strptime(SHOWCASE[9:], '%Y%m%d'):%d %b %Y}" fig, axes = plt.subplots(1, 3, figsize=(13, 5), constrained_layout=True) im0 = axes[0].imshow(wrapped, cmap="twilight", vmin=-np.pi, vmax=np.pi, extent=extent, interpolation="nearest") axes[0].set_title(f"Wrapped phase: one colour cycle\n= {abs(2 * np.pi * RAD2MM):.1f} mm of LOS motion") fig.colorbar(im0, ax=axes[0], shrink=0.8, location="bottom", label="phase (radians)") im1 = axes[1].imshow(show_los, cmap="RdBu", vmin=-60, vmax=60, extent=extent, interpolation="nearest") axes[1].set_title("Unwrapped, as LOS displacement\n(red = away from the satellite)") fig.colorbar(im1, ax=axes[1], shrink=0.8, location="bottom", label="LOS displacement (mm)") im2 = axes[2].imshow(show_coh, cmap="Greys_r", vmin=0, vmax=1, extent=extent, interpolation="nearest") axes[2].set_title("Coherence\n(0 = noise, 1 = perfectly stable)") fig.colorbar(im2, ax=axes[2], shrink=0.8, location="bottom", label="coherence") for ax in axes: geo_axes(ax) label_points(ax, ["Chalco", "Airport (AICM)", "Reference (CU)"]) fig.suptitle(f"One Sentinel-1 interferogram over Mexico City, {pretty} ({days(SHOWCASE)} days)", fontsize=10) fig.savefig(OUT / "01-wrapped-unwrapped-coherence.png") plt.close(fig) fig, ax = plt.subplots(figsize=(9, 3.2)) arc = np.linspace(0, np.pi, 30) for p in pairs: t0, t1 = mdates.date2num(pd.Timestamp(p[:8])), mdates.date2num(pd.Timestamp(p[9:])) ax.plot(t0 + (t1 - t0) * (1 - np.cos(arc)) / 2, days(p) * np.sin(arc), lw=0.9, alpha=0.8, color="#1d6fa3" if days(p) <= 7 else "#2b8a3e" if days(p) <= 13 else "#e8590c") ax.plot(dates, np.zeros(len(dates)), "o", color="black", ms=2.5) for label, color in (("6-day pairs", "#1d6fa3"), ("12-day", "#2b8a3e"), ("18-day", "#e8590c")): ax.plot([], [], color=color, label=label) ax.legend(fontsize=7.5, loc="upper right", ncols=3) ax.set_ylim(-1, 24) ax.set_ylabel("pair length (days)") ax.set_title(f"Small-baseline network: {n_ifg} interferograms linking {len(epochs)} Sentinel-1 acquisitions") ax.xaxis.set_major_formatter(mdates.DateFormatter("%b %Y")) ax.grid(alpha=0.3) fig.savefig(OUT / "02-sbas-network.png") plt.close(fig) fig, axes = plt.subplots(1, 2, figsize=(11, 4.9), constrained_layout=True) im = axes[0].imshow(mean_coh, cmap="Greys_r", vmin=0, vmax=1, extent=extent) axes[0].set_title(f"Mean coherence of {n_ifg} short pairs") fig.colorbar(im, ax=axes[0], shrink=0.85, label="mean coherence") axes[1].imshow(hill, cmap="Greys", vmin=-40, vmax=40, extent=extent) axes[1].imshow(np.ma.masked_where(mask, np.ones_like(mean_coh)), cmap="Oranges", vmin=0, vmax=1.4, alpha=0.8, extent=extent) axes[1].set_title(f"Kept {mask.mean():.0%} of pixels (coherence ≥ 0.35); orange = masked") for ax in axes: geo_axes(ax) fig.savefig(OUT / "03-coherence-mask.png") plt.close(fig) (lon0, lat0), (lon1, lat1) = PROFILE fig, ax = plt.subplots(figsize=(8.6, 7.6), constrained_layout=True) ax.imshow(hill, cmap="Greys", vmin=-40, vmax=40, extent=extent) im = ax.imshow(v_up, cmap="Spectral", norm=TwoSlopeNorm(0, -350, 50), extent=extent, interpolation="nearest") fig.colorbar(im, ax=ax, shrink=0.8, label="vertical velocity (mm/yr); negative = subsidence") ax.plot([lon0, lon1], [lat0, lat1], "--", color="black", lw=1.2) ax.annotate("A", (lon0, lat0), xytext=(-12, 2), textcoords="offset points", fontweight="bold", path_effects=halo) ax.annotate("B", (lon1, lat1), xytext=(5, -10), textcoords="offset points", fontweight="bold", path_effects=halo) label_points(ax) ax.set_title(f"Mexico City vertical ground velocity, Sentinel-1, {dates[0]:%b %Y} – {dates[-1]:%b %Y}") geo_axes(ax) fig.savefig(OUT / "04-vertical-velocity-map.png") plt.close(fig) n = 400 plon, plat = np.linspace(lon0, lon1, n), np.linspace(lat0, lat1, n) pc, pr = inv * (plon, plat) dist_km = np.hypot((plon - lon0) * 111.32 * np.cos(np.radians(19.4)), (plat - lat0) * 110.6) prof_v = ndimage.map_coordinates(np.where(mask, v_up, -9999), [pr - 0.5, pc - 0.5], order=0, cval=-9999) prof_v = np.where(prof_v < -9000, np.nan, prof_v) prof_h = ndimage.map_coordinates(hgt, [pr - 0.5, pc - 0.5], order=1) section("profile A-B") print(f"length {dist_km[-1]:.1f} km; valid samples {np.isfinite(prof_v).mean():.0%}; " f"fastest {np.nanmin(prof_v):.0f} mm/yr at {dist_km[np.nanargmin(prof_v)]:.1f} km; " f"terrain {prof_h.min():.0f}-{prof_h.max():.0f} m") fig, ax = plt.subplots(figsize=(9, 3.8)) ax.plot(dist_km, prof_v, ".", ms=2.5, color="#c92a2a", alpha=0.6) ax.plot(dist_km, pd.Series(prof_v).rolling(9, center=True, min_periods=3).median(), color="#c92a2a", lw=1.6, label="vertical velocity") ax.axhline(0, color="black", lw=0.6) ax.set_xlabel("distance from A (km)") ax.set_ylabel("vertical velocity (mm/yr)", color="#c92a2a") ax2 = ax.twinx() ax2.plot(dist_km, prof_h, color="#495057", lw=1) ax2.set_ylabel("terrain height, DEM (m)", color="#495057") ax.set_title("Profile A–B, Chapultepec to Chalco: subsidence follows the old lakebed clays") ax.grid(alpha=0.3) fig.savefig(OUT / "05-profile.png") plt.close(fig) fig, ax = plt.subplots(figsize=(9, 4.3)) for (name, col), color in zip(series.items(), ["#c92a2a", "#862e9c", "#e8590c", "#1d6fa3", "#2b8a3e", "#495057"]): ax.plot(series.index, col, "o-", ms=2.5, lw=1, color=color, label=f"{name}: {table.set_index('place').loc[name, 'los_mm_yr']:+.0f} mm/yr LOS") ax.set_ylabel("cumulative LOS displacement (mm)") ax.set_title("Displacement time series from the SBAS inversion (relative to the first date and to CU)") ax.legend(fontsize=7.5, loc="lower left") ax.grid(alpha=0.3) ax.xaxis.set_major_formatter(mdates.DateFormatter("%b %Y")) fig.savefig(OUT / "06-timeseries-points.png") plt.close(fig) section("outputs") print("wrote vertical_velocity_mm_yr.tif, los_timeseries_points.csv, data/ifgramStack.h5 and six figures")