Source code: InSAR ground subsidence in Python: one year of Mexico City sinking, from Sentinel-1

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

"""Cross-check the NumPy SBAS inversion in insar_subsidence.py with MintPy.

insar_subsidence.py writes the referenced LiCSAR stack to data/ifgramStack.h5 in MintPy's layout.
Here MintPy inverts the same network (ifgram_inversion), fits the velocity (timeseries2velocity)
and we compare the result pixel by pixel with our NumPy version.
Run with:  python mintpy_crosscheck.py   (after insar_subsidence.py)
"""
import contextlib
import os
from pathlib import Path

import h5py
import mintpy
import numpy as np
from mintpy.cli import ifgram_inversion, timeseries2velocity

HERE = Path(__file__).parent
DATA = HERE / "data"
os.chdir(DATA)
print(f"MintPy {mintpy.__version__}")

# Unweighted least squares: the same estimator as np.linalg.lstsq. MintPy's progress log goes to mintpy.log.
with open("mintpy.log", "w") as log, contextlib.redirect_stdout(log):
    ifgram_inversion.main(["ifgramStack.h5", "--weight-func", "no", "--update"])
    timeseries2velocity.main(["timeseries.h5", "-o", "velocity.h5", "--update"])

with h5py.File("timeseries.h5") as f:
    ts = f["timeseries"][:]                    # metres, positive = towards the satellite
    dates = [d.decode() for d in f["date"][:]]
with h5py.File("velocity.h5") as f:
    v_mintpy = f["velocity"][:] * 1000         # m/yr -> mm/yr
mask = np.load("mask.npy")
v_numpy = np.load("v_los_numpy.npy")
diff = v_mintpy[mask] - v_numpy[mask]
print(f"\nMintPy time series: {ts.shape[0]} dates, {dates[0]} to {dates[-1]}")
print(f"LOS velocity, MintPy: 1st percentile {np.percentile(v_mintpy[mask], 1):.0f} mm/yr, "
      f"median {np.median(v_mintpy[mask]):.0f} mm/yr")
print(f"MintPy minus NumPy over {mask.sum():,} pixels: mean {diff.mean():+.2f} mm/yr, "
      f"max |difference| {np.abs(diff).max():.2f} mm/yr, correlation {np.corrcoef(v_mintpy[mask], v_numpy[mask])[0, 1]:.4f}")