API

The command line wraps these modules; everything can be called from Python.

Running

One frame's end-to-end validation run (plan T21-T25).

run_validation(config) streams the stack once: the calibrated velocity (DISP - CAL, chained) and the station series; the raw velocity, mean coherence and raw station series are independent of the calibration and are cached per stack. It then compares with MIDAS by station class, computes the double-difference sills (raw, calibrated, baseline) with paired bootstrap CIs, applies the gate and writes the outputs (report.write_outputs).

Compatibility with the trade study (e2e_validation/validate_stack.py): the raw cache uses the same key and file names, and the outputs the same names, so trade-study validation/ folders serve as baselines.

RunConfig

Bases: BaseModel

One validation run of one candidate on one frame.

Source code in src/disp2vlm_validation/pipeline.py

class RunConfig(BaseModel):
    """One validation run of one candidate on one frame."""

    label: str
    frame: FrameInputs
    candidate_dir: Path = Field(description="DISP-CAL products of the candidate")
    out_dir: Path
    baseline_dir: Path | None = Field(
        default=None, description="validation output folder of the baseline run"
    )
    gate: Path | None = None
    golden: GoldenCheck | None = None
    window: int = 7
    min_window_pixels: int = 10
    min_epochs: int = 5
    coh_threshold: float = 0.5
    min_common_dates: int = 30
    sill_min_km: float = 50.0
    near_max_km: float = 10.0
    n_boot: int = 1000
    seed: int = 0
    workers: int = 6
    allow_partial: bool = False
    min_span_years: float = 2.0

    model_config = {"extra": "forbid"}

FrameInputs

Bases: BaseModel

Everything about a frame that does not depend on the candidate.

Source code in src/disp2vlm_validation/pipeline.py

class FrameInputs(BaseModel):
    """Everything about a frame that does not depend on the candidate."""

    frame: str
    disp_dir: Path
    water_mask: Path = Field(description="GeoTIFF on the DISP grid, 1 = land")
    los_file: Path = Field(description="3-band LOS ENU GeoTIFF on the DISP grid")
    cache_dir: Path
    ref_pixel: tuple[int, int] | None = Field(
        default=None,
        description="raw velocity reference; default: the last DISP product's",
    )
    exclude_stations: list[str] = Field(default_factory=list)
    held_out_fraction: float = 0.0
    defo_geojson: Path | None = None
    event_geojson: Path | None = None
    pixel_m: float = 30.0

    model_config = {"extra": "forbid"}

GoldenCheck

Bases: BaseModel

Golden-epoch check: the candidate's product for tag vs the golden.

Source code in src/disp2vlm_validation/pipeline.py

class GoldenCheck(BaseModel):
    """Golden-epoch check: the candidate's product for `tag` vs the golden."""

    golden_product: Path
    tag: str = Field(description="<ref>_<sec> dates, e.g. 20220111_20220722")
    executable: str = "cal-disp"

    model_config = {"extra": "forbid"}

load_run_config(path)

RunConfig from YAML.

Source code in src/disp2vlm_validation/pipeline.py

def load_run_config(path: Path | str) -> RunConfig:
    """`RunConfig` from YAML."""
    from .config_io import read_config

    return RunConfig(**read_config(path))

run_validation(cfg, *, use_cache=True)

Validate one candidate on one frame; return the metrics (also written).

Source code in src/disp2vlm_validation/pipeline.py

def run_validation(  # noqa: PLR0912, PLR0915
    cfg: RunConfig, *, use_cache: bool = True
) -> dict[str, Any]:
    """Validate one candidate on one frame; return the metrics (also written)."""
    from .report import write_outputs

    fr = cfg.frame
    out = Path(cfg.out_dir)
    out.mkdir(parents=True, exist_ok=True)
    epochs = list_epochs(
        fr.disp_dir, cfg.candidate_dir, allow_partial=cfg.allow_partial
    )
    run_hash = config_hash(cfg, epochs)
    metrics_path = out / "metrics.json"
    if use_cache and metrics_path.exists():
        old = json.loads(metrics_path.read_text())
        if old.get("config_hash") == run_hash:
            logger.info("cached: %s (config hash %s)", metrics_path, run_hash)
            return old
    with rasterio.open(fr.los_file) as ds:
        profile, los = ds.profile, ds.read().astype(np.float32)
    shape = los.shape[1:]
    land = _read_tif(fr.water_mask).astype(bool)
    if land.shape != shape:
        msg = f"water mask {land.shape} is not on the LOS grid {shape}"
        raise ValueError(msg)
    all_epochs = list_epochs(fr.disp_dir)
    ref_pixel: tuple[int, int] = (
        (int(fr.ref_pixel[0]), int(fr.ref_pixel[1]))
        if fr.ref_pixel
        else product_reference_point(all_epochs[-1].disp)
    )
    logger.info(
        "%s: %d epochs, %d reference-date blocks; raw reference pixel %s",
        fr.frame,
        len(epochs),
        len({e.ref for e in epochs}),
        ref_pixel,
    )

    # GNSS (frozen per frame)
    from rasterio.warp import transform_bounds

    bbox = transform_bounds(
        profile["crs"],
        "EPSG:4326",
        *rasterio.transform.array_bounds(*shape, profile["transform"]),
    )
    midas_df, daily = load_or_fetch_gnss(
        fr.cache_dir, bbox, min_span_years=cfg.min_span_years
    )
    st_all = station_pixels(midas_df, profile["crs"], profile["transform"], los)
    split = split_held_out(
        st_all.index,
        exclude=fr.exclude_stations,
        held_out_fraction=fr.held_out_fraction,
    )
    st = st_all.loc[split.validation].copy()
    sampler = StationSampler(st, shape, land, cfg.window, cfg.min_window_pixels)

    # raw part: cached per stack (independent of the candidate)
    partial = len(epochs) != len(all_epochs)
    cache = (
        out / "_partial_raw"
        if partial
        else Path(fr.cache_dir) / f"stack_{raw_cache_key(all_epochs, cfg, ref_pixel)}"
    )
    cache.mkdir(parents=True, exist_ok=True)
    raw_tif = cache / "raw_velocity_mm_yr.tif"
    coh_tif = cache / "mean_temporal_coherence.tif"
    raw_csv = cache / "raw_station_series_mm.csv"
    need_raw = partial or not all(p.exists() for p in (raw_tif, coh_tif, raw_csv))
    if not need_raw:
        raw_ser_cached = pd.read_csv(raw_csv, index_col=0, parse_dates=True)
        need_raw = not set(st.index) <= set(raw_ser_cached.columns)

    first_ref = pd.Timestamp(min(e.ref for e in epochs)).normalize()
    dates = pd.DatetimeIndex(
        [first_ref] + [pd.Timestamp(e.sec).normalize() for e in epochs]
    )
    t0 = float(np.mean([e.t for e in epochs]))
    fit_cal = VelocityFit(shape, t0)
    fit_raw = VelocityFit(shape, t0) if need_raw else None
    coh_sum = np.zeros(shape) if need_raw else None
    ser_cal = np.zeros((len(dates), len(st)))
    ser_raw = np.zeros((len(dates), len(st))) if need_raw else None
    logger.info(
        "streaming epochs (raw part %s)",
        "computed" if need_raw else f"cached in {cache}",
    )
    for i, (e, (cal_disp, raw, coh)) in enumerate(
        stream_epochs(epochs, ref_pixel if need_raw else None, cfg.workers), 1
    ):
        assert cal_disp is not None
        ser_cal[i] = sampler(fit_cal.add(e.ref, e.t, cal_disp, land)) * 1000.0
        if need_raw:
            assert fit_raw is not None and raw is not None and coh is not None
            ser_raw[i] = sampler(fit_raw.add(e.ref, e.t, raw, land)) * 1000.0  # type: ignore[index]
            np.add(coh_sum, coh, out=coh_sum, where=np.isfinite(coh))
        if i % 25 == 0 or i == len(epochs):
            logger.info("  %d/%d epochs", i, len(epochs))
    if need_raw:
        assert fit_raw is not None and coh_sum is not None
        _write_tif(coh_tif, coh_sum / len(epochs), profile, "mean temporal coherence")
        _write_tif(
            raw_tif,
            np.where(land, fit_raw.velocity_mm_yr(cfg.min_epochs), np.nan),
            profile,
            f"raw LOS velocity (mm/yr), referenced to pixel {ref_pixel}",
        )
        pd.DataFrame(ser_raw, index=dates, columns=st.index).to_csv(raw_csv)
        del fit_raw
    raw_ser = pd.read_csv(raw_csv, index_col=0, parse_dates=True).reindex(index=dates)[
        st.index
    ]
    cal_ser = pd.DataFrame(ser_cal, index=dates, columns=st.index)
    cal_ser.to_csv(out / SERIES_CSV)
    keep = land & (_read_tif(coh_tif) >= cfg.coh_threshold)
    vel_raw = np.where(keep, _read_tif(raw_tif), np.nan).astype(np.float32)
    vel_cal = np.where(keep, fit_cal.velocity_mm_yr(cfg.min_epochs), np.nan).astype(
        np.float32
    )
    del fit_cal
    _write_tif(out / VELOCITY_TIF, vel_cal, profile, "calibrated LOS velocity (mm/yr)")

    # variants: candidate, raw, baseline
    variants = {"calibrated": vel_cal, "raw": vel_raw}
    series = {"calibrated": cal_ser, "raw": raw_ser}
    if cfg.baseline_dir is not None:
        b = Path(cfg.baseline_dir)
        if not ((b / VELOCITY_TIF).exists() and (b / SERIES_CSV).exists()):
            msg = f"baseline {b} has no {VELOCITY_TIF} / {SERIES_CSV}"
            raise FileNotFoundError(msg)
        variants["baseline"] = _read_tif(b / VELOCITY_TIF)
        series["baseline"] = pd.read_csv(
            b / SERIES_CSV, index_col=0, parse_dates=True
        ).reindex(index=dates, columns=st.index)
    for name, vel in variants.items():
        st[f"pix::{name}"] = vel[st["row"].to_numpy(), st["col"].to_numpy()].astype(
            float
        )
        st[f"med::{name}"] = sampler(vel)

    st["class"] = classify_stations(
        st,
        defo_areas=load_polygons(fr.defo_geojson),
        events=load_polygons(fr.event_geojson),
        land=land,
        pixel_m=fr.pixel_m,
    )
    common = st[np.all(np.isfinite(st[[f"pix::{n}" for n in variants]]), axis=1)]
    velocity = {
        name: velocity_stats(
            common[f"pix::{name}"].to_numpy(), common["gps_los_mm_yr"].to_numpy()
        )
        for name in variants
    }
    # the velocity rules see only what the calibration is responsible for
    gate_cfg = load_gate(cfg.gate)
    excluded_classes = set(gate_cfg.velocity_excluded_classes)
    gated = common[~common["class"].isin(excluded_classes)]
    velocity_gate = {
        name: velocity_stats(
            gated[f"pix::{name}"].to_numpy(), gated["gps_los_mm_yr"].to_numpy()
        )
        for name in variants
    }
    by_class = metrics_by_class(
        common["pix::calibrated"].to_numpy(),
        common["gps_los_mm_yr"].to_numpy(),
        common["class"],
        sigma_insar=None,
        sigma_gnss=None,
    )
    # raw is tied to one reference pixel: remove its constant offset (median
    # residual on the gated stations) before comparing class RMSEs with it
    raw_offset = float(
        np.nanmedian(gated["pix::raw"].to_numpy() - gated["gps_los_mm_yr"].to_numpy())
    )
    by_class_raw = metrics_by_class(
        common["pix::raw"].to_numpy() - raw_offset,
        common["gps_los_mm_yr"].to_numpy(),
        common["class"],
    )

    gps_ser = gnss_los_series(daily, st, dates)
    coords = {s: (float(st.at[s, "lon"]), float(st.at[s, "lat"])) for s in st.index}
    dd = double_differences(
        series,
        gps_ser,
        coords,
        min_common_dates=cfg.min_common_dates,
        sill_min_km=cfg.sill_min_km,
        near_max_km=cfg.near_max_km,
        n_boot=cfg.n_boot,
        seed=cfg.seed,
    )

    golden_ok, golden_out = None, "not run"
    if cfg.golden is not None:
        match = [e for e in epochs if e.tag == cfg.golden.tag]
        if match and match[0].cal is not None:
            golden_ok, golden_out = golden_epoch_check(
                cfg.golden.golden_product,
                match[0].cal,
                executable=cfg.golden.executable,
            )
        else:
            golden_out = f"no candidate product for {cfg.golden.tag}"

    has_base = "baseline" in series
    gate_inputs = GateInputs(
        sill_candidate=dd.sills["calibrated"].sill,
        sill_raw=dd.sills["raw"].sill,
        velocity_candidate=velocity_gate["calibrated"],
        sill_baseline=dd.sills["baseline"].sill if has_base else None,
        sill_change_vs_baseline_ci=(
            dd.change_pct("calibrated", "baseline")[1] if has_base else None
        ),
        velocity_baseline=velocity_gate["baseline"] if has_base else None,
        golden_epoch_ok=golden_ok,
        class_metrics=by_class.to_dict(orient="index"),
        class_metrics_raw=by_class_raw.to_dict(orient="index"),
    )
    result = evaluate(gate_inputs, gate_cfg)

    metrics: dict[str, Any] = {
        "label": cfg.label,
        "frame": fr.frame,
        "config_hash": run_hash,
        "package_version": version("disp2vlm-validation"),
        "n_epochs": len(epochs),
        "partial": partial,
        "reference_pixel": list(ref_pixel),
        "n_gnss_stations": len(st),
        "excluded_stations": split.excluded,
        "held_out_stations": split.held_out,
        "velocity_los": velocity,
        "velocity_los_gate": {
            **velocity_gate,
            "excluded_classes": sorted(excluded_classes),
            "n_excluded": int(len(common) - len(gated)),
        },
        "velocity_by_class": by_class.reset_index().to_dict(orient="records"),
        "velocity_by_class_raw": {
            "offset_removed_mm_yr": raw_offset,
            "classes": by_class_raw.reset_index().to_dict(orient="records"),
        },
        "double_difference": {
            name: {k: v for k, v in vars(s).items() if k != "draws"}
            for name, s in dd.sills.items()
        },
        "sill_change_vs_raw": dict(
            zip(("pct", "ci95"), dd.change_pct("calibrated", "raw"), strict=True)
        ),
        "sill_change_vs_baseline": (
            dict(
                zip(
                    ("pct", "ci95"),
                    dd.change_pct("calibrated", "baseline"),
                    strict=True,
                )
            )
            if has_base
            else None
        ),
        "golden_epoch_check": golden_out,
        "verdict": result.verdict,
        "verdict_reasons": result.reasons,
        "notes": result.notes,
        "deltas": result.deltas,
        "sill_min_km": cfg.sill_min_km,
    }
    write_outputs(metrics, st, dd, out)
    metrics_path.write_text(json.dumps(metrics, indent=1, default=str))
    logger.info("%s %s: %s", fr.frame, cfg.label, result.verdict)
    return metrics

VLM validation against independent GNSS (Venti plan T57).

A VLM run (opera_vlm) writes one product per target pair (T0, Tk), all sharing T0, with vertical, vertical_std, east, mode. For every product and every GNSS station inside it:

• pair residuals: product (3 x 3 pixel median at the station) minus the GNSS displacement over the same span (median of the daily solutions within half_window_days of T0 and of Tk), vertical and east;
• station rates: the slope of the product's vertical series against time minus the slope of the GNSS up series at the same epochs; the gate is on these (owner target, 2026-10-11: 1-2 mm/yr);
• sigma realism: z = residual / vertical_std per span bin, coverage of |z| < 1.96 and its robust spread.

Only independent GNSS (UNR daily, MIDAS), frozen per frame in the same cache the DISP-CAL validation uses. Stations in deforming areas or event footprints are reported but not gated (decision 0002).

VlmRunConfig

Bases: BaseModel

One VLM validation run (vlm.yaml).

Source code in src/disp2vlm_validation/vlm.py

class VlmRunConfig(BaseModel):
    """One VLM validation run (``vlm.yaml``)."""

    label: str
    frame: str
    products_dir: Path = Field(description="Folder of the run's VLM products")
    cache_dir: Path = Field(description="Frozen GNSS cache of the frame")
    out_dir: Path
    defo_geojson: Path | None = None
    event_geojson: Path | None = None
    exclude_stations: list[str] = Field(default_factory=list)
    half_window_days: int = 3
    min_pairs: int = 5
    min_span_years: float = 2.5
    min_gnss_span_years: float = 2.0

    model_config = {"extra": "forbid"}

run_vlm(cfg)

Validate one VLM run; write metrics_vlm.json and the residual tables.

Source code in src/disp2vlm_validation/vlm.py

def run_vlm(cfg: VlmRunConfig) -> dict[str, Any]:
    """Validate one VLM run; write ``metrics_vlm.json`` and the residual tables."""
    from shapely import contains_xy, union_all

    products = list_vlm_products(cfg.products_dir)
    with h5py.File(products[0].path, "r") as f:
        x, y = f["x"][()], f["y"][()]
        epsg = int(np.atleast_1d(f.attrs["epsg"])[0])
    from pyproj import Transformer

    lon, lat = Transformer.from_crs(epsg, "EPSG:4326", always_xy=True).transform(
        [x.min(), x.max(), x.min(), x.max()], [y.min(), y.min(), y.max(), y.max()]
    )
    bbox = (min(lon), min(lat), max(lon), max(lat))
    midas, daily = load_or_fetch_gnss(
        cfg.cache_dir, bbox, min_span_years=cfg.min_gnss_span_years
    )
    stations = midas.drop(index=[s for s in cfg.exclude_stations if s in midas.index])
    res = pair_residuals(products, stations, daily, cfg.half_window_days)
    rates = station_rates(res, cfg.min_pairs, cfg.min_span_years)

    polys = load_polygons(cfg.defo_geojson) + load_polygons(cfg.event_geojson)
    st = stations.reindex(rates.index)
    rates["in_area"] = (
        contains_xy(union_all(polys), st["lon"].to_numpy(), st["lat"].to_numpy())
        if polys
        else False
    )
    gated = rates[(rates["mode"] == 1) & ~rates["in_area"]]
    # r2 of a constant series (e.g. no east motion) is undefined, not an error
    with np.errstate(invalid="ignore", divide="ignore"):
        return _metrics(cfg, products, res, rates, gated)

evaluate_vlm(stats, gate=VLM_GATE)

FAIL reasons of the VLM gate on the gated station-rate statistics.

Source code in src/disp2vlm_validation/vlm.py

def evaluate_vlm(stats: dict[str, Any], gate: VlmGate = VLM_GATE) -> list[str]:
    """FAIL reasons of the VLM gate on the gated station-rate statistics."""
    reasons = []
    n = stats.get("n", 0)
    if n < gate.min_stations:
        return [f"{n} gated stations with rates (< {gate.min_stations}): not assessed"]
    if abs(stats["median"]) > gate.rate_abs_median_max_mm_yr:
        reasons.append(
            f"|median| vertical rate {abs(stats['median']):.2f} mm/yr >"
            f" {gate.rate_abs_median_max_mm_yr:g}"
        )
    if stats["nmad"] > gate.rate_nmad_max_mm_yr:
        reasons.append(
            f"vertical rate NMAD {stats['nmad']:.2f} mm/yr >"
            f" {gate.rate_nmad_max_mm_yr:g}"
        )
    return reasons

A small synthetic frame to try the validator on, offline (demo command).

The frame is 40 x 60 pixels of 2 km (EPSG:32615) with a vertical look. Truth is a 10 mm/yr east-west velocity gradient. Each DISP pair adds a random per-date ramp, the long-wavelength error a calibration should remove, plus 0.5 mm of noise; the stack has two reference blocks (21 dates, 12 days apart).

Three DISP-CAL candidates differ only in how much of the ramp they remove: perfect (all of it), half and none. Sixteen GNSS stations carry the true velocity; their MIDAS rates and daily series are written to the frozen GNSS cache, so nothing is downloaded.

write_demo also writes ready-to-run configs in <root>/configs: perfect.yaml (no baseline), half.yaml and none.yaml (both against perfect), and gate.yaml (a per-class target).

demo_config(root, candidate, baseline=None, **overrides)

Run config (as a dict) for one demo candidate.

The demo stack is short (21 dates over 240 days) and small (120 km), so min_common_dates and sill_min_km are lowered from their full-frame defaults; n_boot is lowered to keep it quick. Paths are absolute, so the configs run from any directory.

Source code in src/disp2vlm_validation/demo.py

def demo_config(
    root: Path | str, candidate: str, baseline: str | None = None, **overrides
) -> dict:
    """Run config (as a dict) for one demo candidate.

    The demo stack is short (21 dates over 240 days) and small (120 km), so
    `min_common_dates` and `sill_min_km` are lowered from their full-frame
    defaults; `n_boot` is lowered to keep it quick. Paths are absolute, so the
    configs run from any directory.
    """
    root = Path(root).resolve()
    cfg = {
        "label": candidate,
        "frame": {
            "frame": "F99999",
            "disp_dir": str(root / "disp"),
            "water_mask": str(root / "water_mask.tif"),
            "los_file": str(root / "los.tif"),
            "cache_dir": str(root / "cache"),
            "ref_pixel": list(REF_PIXEL),
            "pixel_m": PIX,
        },
        "candidate_dir": str(root / candidate),
        "out_dir": str(root / "runs" / candidate),
        "baseline_dir": str(root / "runs" / baseline) if baseline else None,
        "min_common_dates": 15,
        "sill_min_km": 30.0,
        "n_boot": 200,
        "workers": 0,
    }
    cfg.update(overrides)
    return cfg

make_demo_frame(root, seed=1)

Write the synthetic frame under root.

Layout: disp/ (DISP pairs), perfect/, half/, none/ (each a products/<ref>_<sec>/ tree of DISP-CAL products, as cal-disp writes it), water_mask.tif, los.tif and cache/ (the GNSS cache).

Returns

Path root.

Source code in src/disp2vlm_validation/demo.py

def make_demo_frame(root: Path | str, seed: int = 1) -> Path:
    """Write the synthetic frame under `root`.

    Layout: ``disp/`` (DISP pairs), ``perfect/``, ``half/``, ``none/`` (each a
    ``products/<ref>_<sec>/`` tree of DISP-CAL products, as cal-disp writes
    it), ``water_mask.tif``, ``los.tif`` and ``cache/`` (the GNSS cache).

    Returns
    -------
    Path
        `root`.

    """
    from pyproj import Transformer

    root = Path(root)
    root.mkdir(parents=True, exist_ok=True)
    rng = np.random.default_rng(seed)
    _, xx = np.mgrid[0:NY, 0:NX]
    x_km = (xx - NX / 2) * PIX / 1000.0
    vel = (-5.0 + 10.0 * xx / NX) / 1000.0  # m/yr, a 10 mm/yr gradient
    t0 = datetime(2020, 1, 1)
    dates = [t0 + timedelta(days=12 * i) for i in range(21)]
    ramp = {d: rng.normal(0, 0.2, 1)[0] * x_km / 1000.0 for d in dates}  # 0.2 mm/km
    ramp[dates[0]] = np.zeros_like(x_km)
    blocks = [(dates[0], dates[1:11]), (dates[10], dates[11:])]
    for ref, secs in blocks:
        for sec in secs:
            yrs = (sec - ref).days / 365.25
            dramp = ramp[sec] - ramp[ref]
            disp = vel * yrs + dramp + rng.normal(0, 0.0005, (NY, NX))
            _h5(
                root / "disp" / _name("L3_DISP-S1", ref, sec),
                displacement=disp,
                temporal_coherence=np.full((NY, NX), 0.9, np.float32),
            )
            tag = f"{ref:%Y%m%d}_{sec:%Y%m%d}"
            for cand, frac in CANDIDATES.items():
                d = root / cand / "products" / tag
                _h5(
                    d / _name("L4_DISP-CAL-S1", ref, sec),
                    calibration=(frac * dramp).astype(np.float32),
                )
                (d / ".done").touch()
    _tif(root / "water_mask.tif", np.ones((NY, NX)))
    _tif(
        root / "los.tif",
        np.stack([np.zeros((NY, NX)), np.zeros((NY, NX)), np.ones((NY, NX))]),
    )

    # 16 stations on a grid, away from the edges
    rows, cols = np.meshgrid([6, 15, 24, 33], [8, 22, 37, 51], indexing="ij")
    to_ll = Transformer.from_crs(CRS, "EPSG:4326", always_xy=True)
    midas, daily = {}, []
    for i, (r, c) in enumerate(zip(rows.ravel(), cols.ravel(), strict=True)):
        sid = f"G{i:03d}"
        lon, lat = to_ll.transform(X0 + (c + 0.5) * PIX, Y0 - (r + 0.5) * PIX)
        vu = vel[r, c] * 1000.0
        midas[sid] = {
            "lon": lon,
            "lat": lat,
            "ve_mm_yr": 0.0,
            "vn_mm_yr": 0.0,
            "vu_mm_yr": vu,
            "sigma_ve_mm_yr": 0.3,
            "sigma_vn_mm_yr": 0.3,
            "sigma_vu_mm_yr": 0.5,
        }
        for day in dates:
            up = vu / 1000.0 * (day - t0).days / 365.25 + rng.normal(0, 0.001)
            daily.append(
                {"station": sid, "date": day, "east": 0.0, "north": 0.0, "up": up}
            )
    gnss = root / "cache"
    gnss.mkdir(exist_ok=True)
    pd.DataFrame.from_dict(midas, orient="index").rename_axis("station").to_csv(
        gnss / "gnss_midas_enu.csv"
    )
    pd.DataFrame(daily).to_csv(gnss / "gnss_enu_daily.csv.gz", index=False)
    return root

write_demo(root, seed=1)

Write the demo frame and its run configs; return the config paths.

Source code in src/disp2vlm_validation/demo.py

def write_demo(root: Path | str, seed: int = 1) -> list[Path]:
    """Write the demo frame and its run configs; return the config paths."""
    root = make_demo_frame(Path(root).resolve(), seed)
    cfg_dir = root / "configs"
    cfg_dir.mkdir(exist_ok=True)
    gate = cfg_dir / "gate.yaml"
    gate.write_text(
        "# A per-class target for the demo: the stable class must stay within\n"
        "# 1 mm/yr of GNSS. Classes and their rules: docs/guide/metrics.md.\n"
        "per_class_targets:\n"
        "  stable: {max_abs_bias_mm_yr: 1.0, max_rmse_mm_yr: 1.0}\n"
    )
    paths = []
    for cand, base in (("perfect", None), ("half", "perfect"), ("none", "perfect")):
        cfg = demo_config(root, cand, base, gate=str(gate))
        path = cfg_dir / f"{cand}.yaml"
        path.write_text(yaml.safe_dump(cfg, sort_keys=False))
        paths.append(path)
    return paths

Gate and requirements

The CalVal release gate as data (Venti PRD section 2.7).

One place for the numbers a run is judged by, so the gate logic (plan T23), the reports (T24) and the traceability table all read the same values. A change here changes what passes; it needs a decision record.

GateThresholds dataclass

Release-gate thresholds; a run FAILS if any condition is violated.

Attributes

sill_must_drop : bool The calibrated double-difference sill must be below the raw one. sill_rise_max_frac : float Candidate vs baseline: a sill rise above this fraction that is also statistically significant fails. sill_min_pair_km : float Station pairs shorter than this are not used for the sill. velocity_rmse_rise_max_mm_yr, velocity_bias_rise_max_mm_yr : float Candidate vs baseline: velocity RMSE or |bias| vs MIDAS may not rise by more than this. bias_statistic : str Which bias the bias rules use: "median" of the station residuals (default since decision 0001) or "mean". A mean is moved by a few outliers, so fixing them on a frame with a bias of the other sign raised |mean| and failed the gate. target_abs_bias_mm_yr : float Release target per frame (reported, not a pass/fail condition), on the same bias statistic. bootstrap_ci : float Confidence level of the station-bootstrap intervals. velocity_excluded_classes : tuple of str Station classes left out of the velocity rules (decision 0002): the deforming areas and event footprints are restored, not calibrated, so they are reported per class instead of gating the calibration. Empty restores the all-station rule.

Source code in src/disp2vlm_validation/requirements.py

@dataclass(frozen=True)
class GateThresholds:
    """Release-gate thresholds; a run FAILS if any condition is violated.

    Attributes
    ----------
    sill_must_drop : bool
        The calibrated double-difference sill must be below the raw one.
    sill_rise_max_frac : float
        Candidate vs baseline: a sill rise above this fraction that is also
        statistically significant fails.
    sill_min_pair_km : float
        Station pairs shorter than this are not used for the sill.
    velocity_rmse_rise_max_mm_yr, velocity_bias_rise_max_mm_yr : float
        Candidate vs baseline: velocity RMSE or |bias| vs MIDAS may not rise
        by more than this.
    bias_statistic : str
        Which bias the bias rules use: ``"median"`` of the station residuals
        (default since decision 0001) or ``"mean"``. A mean is moved by a few
        outliers, so fixing them on a frame with a bias of the other sign
        raised |mean| and failed the gate.
    target_abs_bias_mm_yr : float
        Release target per frame (reported, not a pass/fail condition), on
        the same bias statistic.
    bootstrap_ci : float
        Confidence level of the station-bootstrap intervals.
    velocity_excluded_classes : tuple of str
        Station classes left out of the velocity rules (decision 0002): the
        deforming areas and event footprints are restored, not calibrated,
        so they are reported per class instead of gating the calibration.
        Empty restores the all-station rule.

    """

    sill_must_drop: bool = True
    sill_rise_max_frac: float = 0.05
    sill_min_pair_km: float = 50.0
    velocity_rmse_rise_max_mm_yr: float = 0.2
    velocity_bias_rise_max_mm_yr: float = 0.2
    bias_statistic: str = "median"
    target_abs_bias_mm_yr: float = 1.0
    bootstrap_ci: float = 0.95
    velocity_excluded_classes: tuple[str, ...] = ("defo_area", "near_event")
    independent_gnss_only: bool = field(default=True, init=False)

VlmGate dataclass

VLM gate on per-station vertical rates vs GNSS up (decision 0003, accepted).

Owner target (2026-10-11): VLM vertical rate within 1-2 mm/yr of GNSS. Gated stations: both looks (mode 1), outside deforming areas and event footprints (decision 0002). RMSE is reported, not gated: a handful of local sites (landslides, wells, coast) dominate it. Baseline: LA (F16940 + F18904, 2016-2023), 172 stations: median -0.49, NMAD 1.47, RMSE 2.77 (1.67 without 6 outliers) mm/yr.

Attributes

rate_abs_median_max_mm_yr, rate_nmad_max_mm_yr : float Limits on |median| and NMAD of (VLM rate - GNSS up rate). min_stations : int Fewer gated stations than this: not assessed (FAIL with that reason).

Source code in src/disp2vlm_validation/requirements.py

@dataclass(frozen=True)
class VlmGate:
    """VLM gate on per-station vertical rates vs GNSS up (decision 0003, accepted).

    Owner target (2026-10-11): VLM vertical rate within 1-2 mm/yr of GNSS.
    Gated stations: both looks (mode 1), outside deforming areas and event
    footprints (decision 0002). RMSE is reported, not gated: a handful of
    local sites (landslides, wells, coast) dominate it. Baseline: LA
    (F16940 + F18904, 2016-2023), 172 stations: median -0.49, NMAD 1.47,
    RMSE 2.77 (1.67 without 6 outliers) mm/yr.

    Attributes
    ----------
    rate_abs_median_max_mm_yr, rate_nmad_max_mm_yr : float
        Limits on |median| and NMAD of (VLM rate - GNSS up rate).
    min_stations : int
        Fewer gated stations than this: not assessed (FAIL with that reason).

    """

    rate_abs_median_max_mm_yr: float = 1.0
    rate_nmad_max_mm_yr: float = 2.0
    min_stations: int = 20

The release gate: configuration and evaluation (PRD section 2.7; plan T23).

A run FAILS when any of these holds:

1. the calibrated double-difference sill is not below the raw one;
2. versus the baseline, the sill rises more than sill_rise_max_pct and the rise is significant (the paired bootstrap CI of the change is above 0);
3. versus the baseline, the velocity RMSE or |bias| (vs MIDAS) rises more than velocity_rmse_rise_max_mm_yr / velocity_bias_rise_max_mm_yr; the bias is the median of the station residuals unless bias_statistic says mean (decision 0001);
4. the golden-epoch check fails.

The thresholds come from gate.yaml (per release); the defaults are disp2vlm_validation.requirements.GATE. |bias| <= target_abs_bias_mm_yr and per-class targets are reported, not gating.

GateConfig

Bases: BaseModel

Thresholds of one release (gate.yaml).

Source code in src/disp2vlm_validation/gate.py

class GateConfig(BaseModel):
    """Thresholds of one release (``gate.yaml``)."""

    release: str = Field(default="v0.5", description="Release the gate applies to")
    sill_min_km: float = GATE.sill_min_pair_km
    sill_rise_max_pct: float = 100.0 * GATE.sill_rise_max_frac
    velocity_rmse_rise_max_mm_yr: float = GATE.velocity_rmse_rise_max_mm_yr
    velocity_bias_rise_max_mm_yr: float = GATE.velocity_bias_rise_max_mm_yr
    bias_statistic: Literal["median", "mean"] = GATE.bias_statistic  # type: ignore[assignment]
    target_abs_bias_mm_yr: float = GATE.target_abs_bias_mm_yr
    significance_ci: float = GATE.bootstrap_ci
    velocity_excluded_classes: list[str] = Field(
        default_factory=lambda: list(GATE.velocity_excluded_classes),
        description="Station classes left out of the velocity rules (decision 0002)",
    )
    per_class_targets: dict[str, ClassTarget] = Field(default_factory=dict)

    model_config = {"extra": "forbid"}

ClassTarget

Bases: BaseModel

Reported per-class target (not gating).

Source code in src/disp2vlm_validation/gate.py

class ClassTarget(BaseModel):
    """Reported per-class target (not gating)."""

    max_abs_bias_mm_yr: float | None = None
    max_rmse_mm_yr: float | None = None
    rmse_must_drop_vs_raw: bool = Field(
        default=False,
        description=(
            "The class RMSE must be below raw's, raw with its constant offset"
            " removed (decision 0002): the calibration must improve on DISP there"
        ),
    )
    reference_rmse_mm_yr: float | None = Field(
        default=None,
        ge=0,
        description=(
            "With rmse_must_drop_vs_raw: the accuracy of the GNSS reference for this"
            " class (grid vs stations, from the TS-G1 study; the validator never"
            " reads the grid). The calibrated RMSE may then reach"
            " sqrt(raw^2 + reference^2): the calibration adds no error beyond its"
            " reference (decision 0002 amendment, 2026-10-13)"
        ),
    )
    gating: bool = Field(
        default=False,
        description=(
            "A miss FAILS the run instead of being a note: a frame's benchmark"
            " condition (decision 0002, e.g. Houston's subsiding stations)"
        ),
    )

    model_config = {"extra": "forbid"}

load_gate(path)

GateConfig from a YAML file; defaults when path is None.

Source code in src/disp2vlm_validation/gate.py

def load_gate(path: Path | str | None) -> GateConfig:
    """`GateConfig` from a YAML file; defaults when `path` is None."""
    if path is None:
        return GateConfig()
    from .config_io import read_config

    return GateConfig(**read_config(path))

evaluate(inputs, gate=None)

Apply the gate (see the module docstring).

Source code in src/disp2vlm_validation/gate.py

def evaluate(inputs: GateInputs, gate: GateConfig | None = None) -> GateResult:
    """Apply the gate (see the module docstring)."""
    g = gate or GateConfig()
    label = "median bias" if g.bias_statistic == "median" else "bias"
    reasons: list[str] = []
    notes: list[str] = []
    deltas: dict[str, Any] = {}

    if not inputs.sill_candidate < inputs.sill_raw:
        reasons.append(
            f"DD sill not reduced by calibration ({inputs.sill_candidate:.1f} vs raw"
            f" {inputs.sill_raw:.1f} mm²)"
        )
    deltas["sill_change_vs_raw_pct"] = 100.0 * (
        inputs.sill_candidate / inputs.sill_raw - 1.0
    )

    if inputs.sill_baseline is not None:
        d_sill = 100.0 * (inputs.sill_candidate / inputs.sill_baseline - 1.0)
        ci = inputs.sill_change_vs_baseline_ci
        significant = ci is None or not ci[0] <= 0.0
        deltas["sill_change_vs_baseline_pct"] = d_sill
        deltas["sill_change_vs_baseline_ci"] = ci
        if d_sill > g.sill_rise_max_pct and significant:
            ci_txt = f" [95% CI {ci[0]:+.1f}, {ci[1]:+.1f}]" if ci else ""
            reasons.append(
                f"DD sill up {d_sill:+.1f}%{ci_txt} vs baseline (limit"
                f" {g.sill_rise_max_pct:g}%)"
            )
    if inputs.velocity_baseline is not None:
        cur, ref = inputs.velocity_candidate, inputs.velocity_baseline
        d_rmse = cur["rmse"] - ref["rmse"]
        d_bias = abs(_bias(cur, g.bias_statistic)) - abs(_bias(ref, g.bias_statistic))
        deltas["velocity_rmse_change_mm_yr"] = d_rmse
        deltas["velocity_abs_bias_change_mm_yr"] = d_bias
        if d_rmse > g.velocity_rmse_rise_max_mm_yr:
            reasons.append(
                f"velocity RMSE up {d_rmse:+.2f} mm/yr vs baseline (limit"
                f" {g.velocity_rmse_rise_max_mm_yr:g})"
            )
        if d_bias > g.velocity_bias_rise_max_mm_yr:
            reasons.append(
                f"velocity |{label}| up {d_bias:+.2f} mm/yr vs baseline (limit"
                f" {g.velocity_bias_rise_max_mm_yr:g})"
            )
    if inputs.golden_epoch_ok is False:
        reasons.append("golden-epoch check failed")

    deltas["bias_statistic"] = g.bias_statistic
    key = "median" if g.bias_statistic == "median" else "bias"
    bias = inputs.velocity_candidate.get(key)
    if bias is not None and abs(bias) > g.target_abs_bias_mm_yr:
        notes.append(
            f"|{label}| {abs(bias):.2f} mm/yr above the release target"
            f" {g.target_abs_bias_mm_yr:g} mm/yr (reported, not gating)"
        )
    for cls, target in g.per_class_targets.items():
        m = inputs.class_metrics.get(cls)
        if not m or m.get("n", 0) == 0:
            continue
        misses = _class_target_misses(
            target, m, inputs.class_metrics_raw.get(cls), g.bias_statistic, label
        )
        tag = f"{cls} (benchmark condition)" if target.gating else cls
        (reasons if target.gating else notes).extend(f"{tag}: {x}" for x in misses)
    return GateResult(passed=not reasons, reasons=reasons, notes=notes, deltas=deltas)

golden_epoch_check(golden, candidate, *, executable='cal-disp', timeout_s=1800.0)

Compare a candidate pair product with the golden via cal-disp validate.

Returns (ok, output). The comparison and its tolerance are cal-disp's own (validate GOLDEN TEST), so this package does not re-implement it.

Source code in src/disp2vlm_validation/gate.py

def golden_epoch_check(
    golden: Path | str,
    candidate: Path | str,
    *,
    executable: str = "cal-disp",
    timeout_s: float = 1800.0,
) -> tuple[bool, str]:
    """Compare a candidate pair product with the golden via ``cal-disp validate``.

    Returns ``(ok, output)``. The comparison and its tolerance are cal-disp's
    own (``validate GOLDEN TEST``), so this package does not re-implement it.
    """
    exe = shutil.which(executable)
    if exe is None:
        msg = f"{executable!r} not found on PATH; install cal-disp or pass executable="
        raise FileNotFoundError(msg)
    proc = subprocess.run(
        [exe, "validate", str(golden), str(candidate)],
        capture_output=True,
        text=True,
        timeout=timeout_s,
        check=False,
    )
    return proc.returncode == 0, (proc.stdout + proc.stderr).strip()

Metrics

Validation metrics: velocity statistics and the double-difference sill.

Double differences (geepers.validation): for every station pair, the InSAR-minus-GNSS misfit of the relative LOS series (B - A) is reduced to an RMSE about its mean; gamma = rmse**2 / 2 is the misfit semivariance. Double differences cancel the datum, so raw and calibrated are directly comparable; calibration removes long-wavelength error, so the sill (median gamma of pairs at least sill_min_km apart) must drop.

The sill confidence interval is a station bootstrap: stations are resampled with replacement, each pair weighted by the product of its stations' counts, and the weighted median recomputed. All variants share the same resamples, so their differences are paired. Ported from the trade study (plan T21.1).

DDResult dataclass

Double-difference results of several variants on common station-dates.

Source code in src/disp2vlm_validation/metrics.py

@dataclass
class DDResult:
    """Double-difference results of several variants on common station-dates."""

    pairs: dict[str, pd.DataFrame]
    semivariograms: dict[str, pd.DataFrame]
    sills: dict[str, SillResult]
    sill_min_km: float

    def change_pct(self, a: str, b: str) -> tuple[float, tuple[float, float]]:
        """Sill change of `a` relative to `b` in %, with its paired 95% CI."""
        d = 100.0 * (self.sills[a].sill / self.sills[b].sill - 1.0)
        draws_a, draws_b = self.sills[a].draws, self.sills[b].draws
        if draws_a.size and draws_b.size:
            return d, ci95(100.0 * (draws_a / draws_b - 1.0))
        return d, (float("nan"), float("nan"))

change_pct(a, b)

Sill change of a relative to b in %, with its paired 95% CI.

Source code in src/disp2vlm_validation/metrics.py

def change_pct(self, a: str, b: str) -> tuple[float, tuple[float, float]]:
    """Sill change of `a` relative to `b` in %, with its paired 95% CI."""
    d = 100.0 * (self.sills[a].sill / self.sills[b].sill - 1.0)
    draws_a, draws_b = self.sills[a].draws, self.sills[b].draws
    if draws_a.size and draws_b.size:
        return d, ci95(100.0 * (draws_a / draws_b - 1.0))
    return d, (float("nan"), float("nan"))

SillResult dataclass

Sill of one variant with its station-bootstrap 95% interval.

Source code in src/disp2vlm_validation/metrics.py

@dataclass
class SillResult:
    """Sill of one variant with its station-bootstrap 95% interval."""

    sill: float
    ci_lo: float
    ci_hi: float
    n_pairs: int
    n_stations: int
    near_gamma: float
    median_pair_rmse: float
    draws: np.ndarray = field(repr=False, default_factory=lambda: np.empty(0))

bootstrap_sills(pairs, *, sill_min_km=50.0, near_max_km=10.0, n_boot=1000, seed=0)

Sill (median gamma of far pairs) per variant with paired bootstrap CIs.

pairs are pairwise_differential_rmse tables (station1, station2, distance_km, rmse) per variant.

Source code in src/disp2vlm_validation/metrics.py

def bootstrap_sills(
    pairs: Mapping[str, pd.DataFrame],
    *,
    sill_min_km: float = 50.0,
    near_max_km: float = 10.0,
    n_boot: int = 1000,
    seed: int = 0,
) -> dict[str, SillResult]:
    """Sill (median gamma of far pairs) per variant with paired bootstrap CIs.

    `pairs` are ``pairwise_differential_rmse`` tables (``station1, station2,
    distance_km, rmse``) per variant.
    """
    names: Sequence[str] = list(pairs)
    stations = sorted(
        set().union(
            *(set(df["station1"]) | set(df["station2"]) for df in pairs.values())
        )
    )
    idx = {s: i for i, s in enumerate(stations)}
    rng = np.random.default_rng(seed)
    out: dict[str, SillResult] = {}
    far = {}
    for n in names:
        df = pairs[n]
        gamma = 0.5 * df["rmse"].to_numpy() ** 2
        dist = df["distance_km"].to_numpy()
        f = dist >= sill_min_km
        far[n] = (
            df["station1"][f].map(idx).to_numpy(),
            df["station2"][f].map(idx).to_numpy(),
            gamma[f],
        )
        near = dist < near_max_km
        n_st = len(set(df["station1"]) | set(df["station2"]))
        out[n] = SillResult(
            sill=float(np.median(gamma[f])) if f.any() else float("nan"),
            ci_lo=float("nan"),
            ci_hi=float("nan"),
            n_pairs=len(df),
            n_stations=n_st,
            near_gamma=float(np.median(gamma[near])) if near.any() else float("nan"),
            median_pair_rmse=float(np.median(df["rmse"])) if len(df) else float("nan"),
        )
    if n_boot > 0 and stations:
        draws = {n: np.empty(n_boot) for n in names}
        for b in range(n_boot):
            counts = np.bincount(
                rng.integers(0, len(stations), len(stations)), minlength=len(stations)
            ).astype(float)
            for n in names:
                i, j, g = far[n]
                draws[n][b] = weighted_median(g, counts[i] * counts[j])
        for n in names:
            lo, hi = ci95(draws[n])
            out[n].ci_lo, out[n].ci_hi, out[n].draws = lo, hi, draws[n]
    return out

ci95(x)

2.5 and 97.5 percentiles.

Source code in src/disp2vlm_validation/metrics.py

def ci95(x: np.ndarray) -> tuple[float, float]:
    """2.5 and 97.5 percentiles."""
    return float(np.nanpercentile(x, 2.5)), float(np.nanpercentile(x, 97.5))

double_differences(insar, gps, coords, *, min_common_dates=30, sill_min_km=50.0, near_max_km=10.0, n_bins=12, n_boot=1000, seed=0)

Pairwise DD RMSE, semivariograms and bootstrap sills per variant.

Parameters

insar : mapping of name -> DataFrame dates x stations InSAR LOS series (mm) per variant (raw, calibrated, baseline...). gps : DataFrame dates x stations GNSS LOS series (mm). coords : mapping station -> (lon, lat).

Only station-dates valid in every variant and in GNSS are used, so the variants are compared on the same data.

Source code in src/disp2vlm_validation/metrics.py

def double_differences(
    insar: Mapping[str, pd.DataFrame],
    gps: pd.DataFrame,
    coords: Mapping[str, tuple[float, float]],
    *,
    min_common_dates: int = 30,
    sill_min_km: float = 50.0,
    near_max_km: float = 10.0,
    n_bins: int = 12,
    n_boot: int = 1000,
    seed: int = 0,
) -> DDResult:
    """Pairwise DD RMSE, semivariograms and bootstrap sills per variant.

    Parameters
    ----------
    insar : mapping of name -> DataFrame
        dates x stations InSAR LOS series (mm) per variant (raw, calibrated,
        baseline...).
    gps : DataFrame
        dates x stations GNSS LOS series (mm).
    coords : mapping
        station -> (lon, lat).

    Only station-dates valid in every variant and in GNSS are used, so the
    variants are compared on the same data.

    """
    from geepers.validation import misfit_semivariogram, pairwise_differential_rmse

    valid = gps.notna()
    for df in insar.values():
        valid &= df.notna()
    pairs, svs = {}, {}
    for name, df in insar.items():
        merged = {
            s: pd.DataFrame({"los_gps": gps[s], "los_insar": df[s]})[valid[s]]
            for s in gps.columns
            if valid[s].sum() >= min_common_dates
        }
        rmse_df = pairwise_differential_rmse(
            merged, {s: coords[s] for s in merged}, min_common_dates=min_common_dates
        )
        pairs[name] = rmse_df
        svs[name] = (
            misfit_semivariogram(rmse_df, n_bins=n_bins, log_bins=True)
            if len(rmse_df)
            else pd.DataFrame()
        )
    sills = bootstrap_sills(
        pairs,
        sill_min_km=sill_min_km,
        near_max_km=near_max_km,
        n_boot=n_boot,
        seed=seed,
    )
    return DDResult(
        pairs=pairs, semivariograms=svs, sills=sills, sill_min_km=sill_min_km
    )

velocity_stats(insar, gps)

Residual statistics of InSAR - GNSS (same units as the inputs).

Source code in src/disp2vlm_validation/metrics.py

def velocity_stats(insar: np.ndarray, gps: np.ndarray) -> dict[str, float]:
    """Residual statistics of InSAR - GNSS (same units as the inputs)."""
    insar = np.asarray(insar, float)
    gps = np.asarray(gps, float)
    ok = np.isfinite(insar) & np.isfinite(gps)
    res = insar[ok] - gps[ok]
    n = int(ok.sum())
    if n < 2:
        return {"n": n}
    return {
        "n": n,
        "bias": float(np.mean(res)),
        "median": float(np.median(res)),
        "rmse": float(np.sqrt(np.mean(res**2))),
        "std": float(np.std(res, ddof=1)),
        "nmad": float(1.4826 * np.median(np.abs(res - np.median(res)))),
        "r2": float(np.corrcoef(gps[ok], insar[ok])[0, 1] ** 2),
    }

weighted_median(x, w)

Weighted median (lower); NaN when no weight is positive.

Source code in src/disp2vlm_validation/metrics.py

def weighted_median(x: np.ndarray, w: np.ndarray) -> float:
    """Weighted median (lower); NaN when no weight is positive."""
    keep = w > 0
    x, w = x[keep], w[keep]
    if x.size == 0:
        return float("nan")
    o = np.argsort(x)
    cw = np.cumsum(w[o])
    return float(x[o][np.searchsorted(cw, 0.5 * cw[-1])])

Station classes and per-class metrics (PRD section 3.4; plan T22).

A station gets the first class that applies, in this order:

1. near_event: inside an event area of the event GeoJSON;
2. defo_area: inside a deforming area of the defo GeoJSON;
3. subsiding: GNSS vertical rate below subsiding_mm_yr;
4. coastal: within coast_km of the sea: a water body of the DISP water mask of at least sea_min_km2 (bays and tidal channels joined to it count; ponds, lakes and rivers do not; sea_min_km2=0 is any water);
5. stable: everything else.

Metrics per class compare InSAR and GNSS velocities and, when sigmas are given, report sigma realism: z = residual / sqrt(sigma_insar^2 + sigma_gnss^2) should have std near 1 and ~95% of |z| below 1.96 (R-E4, report only).

classify_stations(stations, *, defo_areas=(), events=(), land=None, pixel_m=30.0, subsiding_mm_yr=-3.0, coast_km=10.0, sea_min_km2=500.0)

Class of every station (see the module docstring for the order).

stations needs lon, lat, vu_mm_yr and, for the coastal class, row, col on the grid of land (True = land).

Source code in src/disp2vlm_validation/classes.py

def classify_stations(
    stations: pd.DataFrame,
    *,
    defo_areas: Sequence = (),
    events: Sequence = (),
    land: np.ndarray | None = None,
    pixel_m: float = 30.0,
    subsiding_mm_yr: float = -3.0,
    coast_km: float = 10.0,
    sea_min_km2: float = 500.0,
) -> pd.Series:
    """Class of every station (see the module docstring for the order).

    `stations` needs ``lon, lat, vu_mm_yr`` and, for the coastal class,
    ``row, col`` on the grid of `land` (True = land).
    """
    cls = pd.Series("stable", index=stations.index, dtype=object)
    coastal = (
        _distance_to_water_km(stations, land, pixel_m, sea_min_km2) <= coast_km
        if land is not None
        else np.zeros(len(stations), dtype=bool)
    )
    subsiding = stations["vu_mm_yr"].to_numpy() < subsiding_mm_yr
    # assign from the lowest priority up, so the first class in CLASS_ORDER wins
    cls[coastal] = "coastal"
    cls[subsiding] = "subsiding"
    cls[_inside(stations, defo_areas)] = "defo_area"
    cls[_inside(stations, events)] = "near_event"
    return cls

load_polygons(path)

Shapely geometries of a GeoJSON FeatureCollection (EPSG:4326).

Source code in src/disp2vlm_validation/classes.py

def load_polygons(path: Path | str | None) -> list:
    """Shapely geometries of a GeoJSON FeatureCollection (EPSG:4326)."""
    if path is None:
        return []
    from shapely.geometry import shape

    data = json.loads(Path(path).read_text())
    return [shape(f["geometry"]) for f in data.get("features", [])]

metrics_by_class(insar, gnss, classes, *, sigma_insar=None, sigma_gnss=None)

class, n, bias, median, rmse, corr, std_z, coverage95 per class and all.

bias is the mean residual, median the median one.

Source code in src/disp2vlm_validation/classes.py

def metrics_by_class(
    insar: np.ndarray,
    gnss: np.ndarray,
    classes: pd.Series,
    *,
    sigma_insar: np.ndarray | float | None = None,
    sigma_gnss: np.ndarray | float | None = None,
) -> pd.DataFrame:
    """``class, n, bias, median, rmse, corr, std_z, coverage95`` per class and ``all``.

    ``bias`` is the mean residual, ``median`` the median one.
    """
    insar = np.asarray(insar, float)
    gnss = np.asarray(gnss, float)
    ok = np.isfinite(insar) & np.isfinite(gnss)
    res = insar - gnss
    z = (
        normalized_residuals(res, sigma_insar, sigma_gnss)
        if sigma_insar is not None and sigma_gnss is not None
        else np.full(res.shape, np.nan)
    )
    labels = classes.to_numpy()
    rows = []
    for name in (*CLASS_ORDER, "all"):
        sel = ok & ((labels == name) if name != "all" else True)
        n = int(sel.sum())
        row = {"class": name, "n": n}
        if n >= 1:
            r = res[sel]
            row["bias"] = float(np.mean(r))
            row["median"] = float(np.median(r))
            row["rmse"] = float(np.sqrt(np.mean(r**2)))
        if n >= 3 and np.std(insar[sel]) > 0 and np.std(gnss[sel]) > 0:
            row["corr"] = float(np.corrcoef(insar[sel], gnss[sel])[0, 1])
        zs = z[sel]
        zs = zs[np.isfinite(zs)]
        if zs.size >= 3:
            row["std_z"] = float(np.std(zs, ddof=1))
            row["coverage95"] = float(np.mean(np.abs(zs) <= 1.96))
        rows.append(row)
    return pd.DataFrame(rows).set_index("class")

normalized_residuals(residual, sigma_insar, sigma_gnss)

residual / sqrt(sigma_insar^2 + sigma_gnss^2) (NaN where sigma is 0).

Source code in src/disp2vlm_validation/classes.py

def normalized_residuals(
    residual: np.ndarray,
    sigma_insar: np.ndarray | float,
    sigma_gnss: np.ndarray | float,
) -> np.ndarray:
    """``residual / sqrt(sigma_insar^2 + sigma_gnss^2)`` (NaN where sigma is 0)."""
    s = np.sqrt(
        np.asarray(sigma_insar, float) ** 2 + np.asarray(sigma_gnss, float) ** 2
    )
    with np.errstate(divide="ignore", invalid="ignore"):
        return np.where(s > 0, np.asarray(residual, float) / s, np.nan)

sea_mask(land, pixel_m, sea_min_km2)

Water bodies of at least sea_min_km2 (8-connected); all water if 0.

On the four validation frames the seas are 900-15 000 km^2 (Gulf of Mexico, Galveston Bay, Pacific, Atlantic and Long Island Sound) while the largest inland water body is ~300 km^2 (Oklahoma City has ~20 000 small ones), so 500 km^2 separates them.

Source code in src/disp2vlm_validation/classes.py

def sea_mask(land: np.ndarray, pixel_m: float, sea_min_km2: float) -> np.ndarray:
    """Water bodies of at least `sea_min_km2` (8-connected); all water if 0.

    On the four validation frames the seas are 900-15 000 km^2 (Gulf of
    Mexico, Galveston Bay, Pacific, Atlantic and Long Island Sound) while the
    largest inland water body is ~300 km^2 (Oklahoma City has ~20 000 small
    ones), so 500 km^2 separates them.
    """
    from scipy.ndimage import label

    water = ~np.asarray(land, dtype=bool)
    if sea_min_km2 <= 0 or not water.any():
        return water
    lab, n = label(water, structure=np.ones((3, 3), dtype=bool))
    area_km2 = np.bincount(lab.ravel(), minlength=n + 1) * (pixel_m / 1000.0) ** 2
    big = np.flatnonzero(area_km2 >= sea_min_km2)
    big = big[big > 0]
    return np.isin(lab, big)

Data

DISP stacks, their DISP-CAL products, and chained per-pixel velocities.

Ported from trade studies/e2e_validation/validate_stack.py (plan T21.2, T21.4) without hard-coded paths. A stack is the set of DISP-S1 products of a frame; each epoch may have a DISP-CAL product with the same reference and secondary dates. Calibrated displacement is DISP - calibration.

Epochs are chained across reference-date blocks as opera_utils.disp.create_rebased_displacement does: the last cumulative value of a block is the offset of the next one.

Epoch dataclass

One DISP product and, if found, its DISP-CAL product.

Source code in src/disp2vlm_validation/stack.py

@dataclass(frozen=True)
class Epoch:
    """One DISP product and, if found, its DISP-CAL product."""

    ref: datetime
    sec: datetime
    disp: Path
    cal: Path | None = None

    @property
    def tag(self) -> str:
        """``<ref>_<sec>`` dates, the key used to match DISP and DISP-CAL."""
        return f"{self.ref:%Y%m%d}_{self.sec:%Y%m%d}"

    @property
    def t(self) -> float:
        """Secondary date as a decimal year."""
        return decimal_year(self.sec)

t property

Secondary date as a decimal year.

tag property

<ref>_<sec> dates, the key used to match DISP and DISP-CAL.

VelocityFit

Streaming per-pixel OLS of chained cumulative displacement against time.

Buffers are preallocated and updated in place: fresh full-frame temporaries every epoch cost more in page faults than the arithmetic.

Source code in src/disp2vlm_validation/stack.py

class VelocityFit:
    """Streaming per-pixel OLS of chained cumulative displacement against time.

    Buffers are preallocated and updated in place: fresh full-frame
    temporaries every epoch cost more in page faults than the arithmetic.
    """

    def __init__(self, shape: tuple[int, int], t0: float) -> None:
        self.t0 = t0
        self.n, self.st, self.st2, self.sy, self.sty = (
            np.zeros(shape) for _ in range(5)
        )
        self.offset, self.cum, self.tmp = (np.zeros(shape) for _ in range(3))
        self.ok = np.zeros(shape, dtype=bool)
        self.block: object = None

    def add(
        self, block_key: object, t: float, disp_m: np.ndarray, mask: np.ndarray
    ) -> np.ndarray:
        """Add an epoch of block `block_key` (its reference date); return cumulative."""
        if block_key != self.block:
            if self.block is not None:
                np.copyto(
                    self.offset, self.cum
                )  # last cumulative of the previous block
            self.block = block_key
        np.add(self.offset, disp_m, out=self.cum)
        np.isfinite(self.cum, out=self.ok)
        np.logical_and(self.ok, mask, out=self.ok)
        tt = t - self.t0
        np.add(self.n, 1.0, out=self.n, where=self.ok)
        np.add(self.st, tt, out=self.st, where=self.ok)
        np.add(self.st2, tt * tt, out=self.st2, where=self.ok)
        np.add(self.sy, self.cum, out=self.sy, where=self.ok)
        np.multiply(self.cum, tt, out=self.tmp)
        np.add(self.sty, self.tmp, out=self.sty, where=self.ok)
        return self.cum

    def velocity_mm_yr(self, min_epochs: int = 5) -> np.ndarray:
        """OLS slope in mm/yr; NaN with fewer than `min_epochs` valid epochs."""
        den = self.n * self.st2 - self.st**2
        good = (self.n >= min_epochs) & (den > 0)
        v = np.full(self.n.shape, np.nan)
        v[good] = (self.n[good] * self.sty[good] - self.st[good] * self.sy[good]) / den[
            good
        ]
        return (v * 1000.0).astype(np.float32)

add(block_key, t, disp_m, mask)

Add an epoch of block block_key (its reference date); return cumulative.

Source code in src/disp2vlm_validation/stack.py

def add(
    self, block_key: object, t: float, disp_m: np.ndarray, mask: np.ndarray
) -> np.ndarray:
    """Add an epoch of block `block_key` (its reference date); return cumulative."""
    if block_key != self.block:
        if self.block is not None:
            np.copyto(
                self.offset, self.cum
            )  # last cumulative of the previous block
        self.block = block_key
    np.add(self.offset, disp_m, out=self.cum)
    np.isfinite(self.cum, out=self.ok)
    np.logical_and(self.ok, mask, out=self.ok)
    tt = t - self.t0
    np.add(self.n, 1.0, out=self.n, where=self.ok)
    np.add(self.st, tt, out=self.st, where=self.ok)
    np.add(self.st2, tt * tt, out=self.st2, where=self.ok)
    np.add(self.sy, self.cum, out=self.sy, where=self.ok)
    np.multiply(self.cum, tt, out=self.tmp)
    np.add(self.sty, self.tmp, out=self.sty, where=self.ok)
    return self.cum

velocity_mm_yr(min_epochs=5)

OLS slope in mm/yr; NaN with fewer than min_epochs valid epochs.

Source code in src/disp2vlm_validation/stack.py

def velocity_mm_yr(self, min_epochs: int = 5) -> np.ndarray:
    """OLS slope in mm/yr; NaN with fewer than `min_epochs` valid epochs."""
    den = self.n * self.st2 - self.st**2
    good = (self.n >= min_epochs) & (den > 0)
    v = np.full(self.n.shape, np.nan)
    v[good] = (self.n[good] * self.sty[good] - self.st[good] * self.sy[good]) / den[
        good
    ]
    return (v * 1000.0).astype(np.float32)

apply_calibration(epoch)

Calibrated displacement DISP - calibration of one epoch (plan T21.4).

Source code in src/disp2vlm_validation/stack.py

def apply_calibration(epoch: Epoch) -> np.ndarray:
    """Calibrated displacement ``DISP - calibration`` of one epoch (plan T21.4)."""
    if epoch.cal is None:
        msg = f"epoch {epoch.tag} has no DISP-CAL product"
        raise ValueError(msg)
    calibrated, _, _ = read_epoch(epoch)
    assert calibrated is not None
    return calibrated

decimal_year(t)

Decimal year with 365.25-day years from 2000-01-01.

Source code in src/disp2vlm_validation/stack.py

def decimal_year(t: datetime) -> float:
    """Decimal year with 365.25-day years from 2000-01-01."""
    return 2000.0 + (t - _EPOCH0).total_seconds() / 86400.0 / DAYS_PER_YEAR

list_epochs(disp_dir, cal_dir=None, *, allow_partial=False)

Epochs of a stack sorted by secondary date, each with its CAL product.

Raises

FileNotFoundError If disp_dir has no DISP product, or a CAL product is missing and allow_partial is False (then those epochs are skipped instead).

Source code in src/disp2vlm_validation/stack.py

def list_epochs(
    disp_dir: Path | str,
    cal_dir: Path | str | None = None,
    *,
    allow_partial: bool = False,
) -> list[Epoch]:
    """Epochs of a stack sorted by secondary date, each with its CAL product.

    Raises
    ------
    FileNotFoundError
        If `disp_dir` has no DISP product, or a CAL product is missing and
        `allow_partial` is False (then those epochs are skipped instead).

    """
    disp_files = sorted(Path(disp_dir).glob(DISP_GLOB))
    if not disp_files:
        msg = f"no {DISP_GLOB} in {disp_dir}"
        raise FileNotFoundError(msg)
    cal_index = _index_cal_products(Path(cal_dir)) if cal_dir is not None else {}
    epochs: list[Epoch] = []
    missing: list[str] = []
    for f in disp_files:
        ref, sec = parse_dates(f.name)
        e = Epoch(ref=ref, sec=sec, disp=f)
        if cal_dir is not None:
            cal = cal_index.get(e.tag)
            if cal is None:
                missing.append(e.tag)
                continue
            e = Epoch(ref=ref, sec=sec, disp=f, cal=cal)
        epochs.append(e)
    if missing and not allow_partial:
        msg = (
            f"{len(missing)} epochs have no DISP-CAL product in {cal_dir}"
            f" (first: {missing[:3]}); pass allow_partial to skip them"
        )
        raise FileNotFoundError(msg)
    epochs.sort(key=lambda e: e.sec)
    return epochs

parse_dates(name)

Return the reference and secondary datetimes of a DISP or DISP-CAL file.

Source code in src/disp2vlm_validation/stack.py

def parse_dates(name: str) -> tuple[datetime, datetime]:
    """Return the reference and secondary datetimes of a DISP or DISP-CAL file."""
    m = PAIR_DATES.search(name)
    if m is None:
        msg = f"no _<ref>Z_<sec>Z_ dates in {name!r}"
        raise ValueError(msg)
    ref, sec = m.groups()
    return (
        datetime.strptime(ref, "%Y%m%dT%H%M%S"),
        datetime.strptime(sec, "%Y%m%dT%H%M%S"),
    )

read_epoch(epoch, ref_pixel=None)

Return (calibrated, raw, coherence) for one epoch (float32, metres).

calibrated is DISP - calibration (None without a CAL product); raw (only with ref_pixel) is DISP referenced to that pixel, and coherence its temporal coherence.

Source code in src/disp2vlm_validation/stack.py

def read_epoch(
    epoch: Epoch, ref_pixel: tuple[int, int] | None = None
) -> tuple[np.ndarray | None, np.ndarray | None, np.ndarray | None]:
    """Return ``(calibrated, raw, coherence)`` for one epoch (float32, metres).

    ``calibrated`` is ``DISP - calibration`` (None without a CAL product);
    ``raw`` (only with `ref_pixel`) is DISP referenced to that pixel, and
    ``coherence`` its temporal coherence.
    """
    disp = _read2d(epoch.disp, "displacement")
    calibrated = None
    if epoch.cal is not None:
        calibrated = (disp - _read2d(epoch.cal, "calibration")).astype(np.float32)
    raw = coh = None
    if ref_pixel is not None:
        raw = (disp - disp[tuple(ref_pixel)]).astype(np.float32)
        coh = _read2d(epoch.disp, "temporal_coherence").astype(np.float32)
    return calibrated, raw, coh

stream_epochs(epochs, ref_pixel=None, workers=4)

Yield (epoch, read_epoch(epoch)) in order, reading ahead in processes.

workers=0 reads in this process (tests, small stacks).

Source code in src/disp2vlm_validation/stack.py

def stream_epochs(
    epochs: Sequence[Epoch],
    ref_pixel: tuple[int, int] | None = None,
    workers: int = 4,
) -> Iterator[
    tuple[Epoch, tuple[np.ndarray | None, np.ndarray | None, np.ndarray | None]]
]:
    """Yield ``(epoch, read_epoch(epoch))`` in order, reading ahead in processes.

    ``workers=0`` reads in this process (tests, small stacks).
    """
    if workers <= 0:
        for e in epochs:
            yield e, read_epoch(e, ref_pixel)
        return
    with ProcessPoolExecutor(workers) as pool:
        futs = [pool.submit(read_epoch, e, ref_pixel) for e in epochs[:workers]]
        for i, e in enumerate(epochs):
            res = futs[i].result()
            futs[i] = None  # type: ignore[call-overload]
            if i + workers < len(epochs):
                futs.append(pool.submit(read_epoch, epochs[i + workers], ref_pixel))
            yield e, res

Independent GNSS: UNR daily stations, MIDAS rates, LOS projection, sampling.

Only daily station positions and their MIDAS rates are used, never the UNR gridded field the calibration is fitted to (GATE.independent_gnss_only). Ported from trade studies/e2e_validation/validate_stack.py (plan T21.2, T21.3).

StationSampler

Median of the valid (land, finite) pixels in a window around each station.

Source code in src/disp2vlm_validation/gnss.py

class StationSampler:
    """Median of the valid (land, finite) pixels in a window around each station."""

    def __init__(
        self,
        stations: pd.DataFrame,
        shape: tuple[int, int],
        land: np.ndarray,
        window: int = 7,
        min_pixels: int = 10,
    ) -> None:
        h = window // 2
        ny, nx = shape
        self.stations = list(stations.index)
        self.windows = [
            (
                slice(max(r - h, 0), min(r + h + 1, ny)),
                slice(max(c - h, 0), min(c + h + 1, nx)),
            )
            for r, c in zip(stations["row"], stations["col"], strict=True)
        ]
        self.land = [land[w] for w in self.windows]
        self.min_pixels = min_pixels

    def __call__(self, arr: np.ndarray) -> np.ndarray:
        """One value per station; NaN with fewer than `min_pixels` valid pixels."""
        out = np.full(len(self.windows), np.nan)
        for i, (w, land) in enumerate(zip(self.windows, self.land, strict=True)):
            v = arr[w][land]
            v = v[np.isfinite(v)]
            if v.size >= self.min_pixels:
                out[i] = np.median(v)
        return out

__call__(arr)

One value per station; NaN with fewer than min_pixels valid pixels.

Source code in src/disp2vlm_validation/gnss.py

def __call__(self, arr: np.ndarray) -> np.ndarray:
    """One value per station; NaN with fewer than `min_pixels` valid pixels."""
    out = np.full(len(self.windows), np.nan)
    for i, (w, land) in enumerate(zip(self.windows, self.land, strict=True)):
        v = arr[w][land]
        v = v[np.isfinite(v)]
        if v.size >= self.min_pixels:
            out[i] = np.median(v)
    return out

StationSplit dataclass

Validation stations and the held-out ones (for fitting, e.g. TS-G1's k).

Source code in src/disp2vlm_validation/gnss.py

@dataclass(frozen=True)
class StationSplit:
    """Validation stations and the held-out ones (for fitting, e.g. TS-G1's k)."""

    validation: list[str]
    held_out: list[str]
    excluded: list[str]

fetch_gnss(bbox, min_span_years=2.0, max_workers=8)

UNR stations in bbox (W, S, E, N): MIDAS ENU rates and daily ENU series.

Returns

midas : DataFrame Index station; lon, lat, v{e,n,u}_mm_yr, sigma_v{e,n,u}_mm_yr, n_gps_epochs, span_years, first, last (IGS20). daily : DataFrame station, date, east, north, up in metres, zeroed by the mean.

Source code in src/disp2vlm_validation/gnss.py

def fetch_gnss(
    bbox: Sequence[float], min_span_years: float = 2.0, max_workers: int = 8
) -> tuple[pd.DataFrame, pd.DataFrame]:
    """UNR stations in `bbox` (W, S, E, N): MIDAS ENU rates and daily ENU series.

    Returns
    -------
    midas : DataFrame
        Index ``station``; ``lon, lat, v{e,n,u}_mm_yr, sigma_v{e,n,u}_mm_yr,
        n_gps_epochs, span_years, first, last`` (IGS20).
    daily : DataFrame
        ``station, date, east, north, up`` in metres, zeroed by the mean.

    """
    from geepers.gps_sources.unr import UnrSource
    from geepers.midas import midas

    gps = UnrSource().timeseries_many(
        bbox=tuple(bbox), frame="ENU", zero_by="mean", max_workers=max_workers
    )
    rows = {}
    for sid, group in gps.groupby("id"):
        d = group.sort_values("date")
        span = (d["date"].max() - d["date"].min()).days / DAYS_PER_YEAR
        if len(d) < 2 or span < min_span_years:
            continue
        t = (d["date"] - d["date"].iloc[0]).dt.total_seconds().to_numpy() / 86400.0
        t = t / DAYS_PER_YEAR
        row = {"lon": float(d["lon"].iloc[0]), "lat": float(d["lat"].iloc[0])}
        for comp, k in (("east", "e"), ("north", "n"), ("up", "u")):
            r = midas(t, d[comp].to_numpy())
            row[f"v{k}_mm_yr"] = r.velocity * 1000.0
            row[f"sigma_v{k}_mm_yr"] = r.velocity_uncertainty * 1000.0
        if all(np.isfinite(row[f"v{k}_mm_yr"]) for k in "enu"):
            row.update(
                n_gps_epochs=len(d),
                span_years=span,
                first=d["date"].min().date(),
                last=d["date"].max().date(),
            )
            rows[sid] = row
    midas_df = pd.DataFrame.from_dict(rows, orient="index").rename_axis("station")
    daily = pd.DataFrame(gps[["id", "date", "east", "north", "up"]]).rename(
        columns={"id": "station"}
    )
    return midas_df, daily[daily["station"].isin(midas_df.index)]

gnss_los_series(daily, stations, dates)

GNSS LOS displacement (mm) at the acquisition dates: dates x stations.

Source code in src/disp2vlm_validation/gnss.py

def gnss_los_series(
    daily: pd.DataFrame, stations: pd.DataFrame, dates: pd.DatetimeIndex
) -> pd.DataFrame:
    """GNSS LOS displacement (mm) at the acquisition dates: dates x stations."""
    d = daily[daily["station"].isin(stations.index)].copy()
    d["date"] = pd.to_datetime(d["date"]).dt.normalize()
    e, n, u = (d["station"].map(stations[f"los_{k}"]) for k in "enu")
    d["los"] = (e * d["east"] + n * d["north"] + u * d["up"]) * 1000.0
    wide = d.pivot_table(index="date", columns="station", values="los", aggfunc="mean")
    return wide.reindex(index=dates, columns=stations.index)

load_or_fetch_gnss(cache_dir, bbox, *, min_span_years=2.0, refresh=False)

MIDAS and daily GNSS, frozen in cache_dir on first use.

Every run of a frame is then compared with the same GNSS numbers; the trade-study cache layout (gnss_midas_enu.csv, gnss_enu_daily.csv.gz) is read as is.

Source code in src/disp2vlm_validation/gnss.py

def load_or_fetch_gnss(
    cache_dir: Path | str,
    bbox: Sequence[float],
    *,
    min_span_years: float = 2.0,
    refresh: bool = False,
) -> tuple[pd.DataFrame, pd.DataFrame]:
    """MIDAS and daily GNSS, frozen in `cache_dir` on first use.

    Every run of a frame is then compared with the same GNSS numbers; the
    trade-study cache layout (``gnss_midas_enu.csv``,
    ``gnss_enu_daily.csv.gz``) is read as is.
    """
    cache = Path(cache_dir)
    midas_csv = cache / "gnss_midas_enu.csv"
    daily_csv = cache / "gnss_enu_daily.csv.gz"
    if refresh or not (midas_csv.exists() and daily_csv.exists()):
        cache.mkdir(parents=True, exist_ok=True)
        midas_df, daily = fetch_gnss(bbox, min_span_years)
        midas_df.to_csv(midas_csv)
        daily.to_csv(daily_csv, index=False)
    return (
        pd.read_csv(midas_csv, index_col=0),
        pd.read_csv(daily_csv, parse_dates=["date"]),
    )

split_held_out(stations, *, exclude=(), held_out_fraction=0.0, seed='disp2vlm')

Split stations into validation and held-out sets, deterministically.

exclude (config, not code: e.g. UH01, a bad record) never enters either set. The held-out choice hashes each station id with seed, so it is stable across runs and machines and independent of the station order: stations used to fit parameters (TS-G1) can never enter validation.

Source code in src/disp2vlm_validation/gnss.py

def split_held_out(
    stations: Iterable[str],
    *,
    exclude: Iterable[str] = (),
    held_out_fraction: float = 0.0,
    seed: str = "disp2vlm",
) -> StationSplit:
    """Split stations into validation and held-out sets, deterministically.

    `exclude` (config, not code: e.g. ``UH01``, a bad record) never enters
    either set. The held-out choice hashes each station id with `seed`, so it
    is stable across runs and machines and independent of the station order:
    stations used to fit parameters (TS-G1) can never enter validation.
    """
    if not 0.0 <= held_out_fraction < 1.0:
        msg = "held_out_fraction must be in [0, 1)"
        raise ValueError(msg)
    excluded = sorted(set(exclude))
    validation: list[str] = []
    held: list[str] = []
    for s in stations:
        if s in excluded:
            continue
        h = int(hashlib.sha256(f"{seed}:{s}".encode()).hexdigest()[:8], 16) / 0xFFFFFFFF
        (held if h < held_out_fraction else validation).append(s)
    return StationSplit(sorted(validation), sorted(held), excluded)

station_pixels(stations, crs, transform, los)

Locate stations on the DISP grid and project MIDAS rates onto the LOS.

los is the 3-band (east, north, up) LOS unit-vector stack of the DISP-S1-STATIC product on the same grid. Stations outside the frame or without a look are dropped. Adds row, col, los_{e,n,u}, gps_hz_los_mm_yr, gps_los_mm_yr.

Source code in src/disp2vlm_validation/gnss.py

def station_pixels(
    stations: pd.DataFrame, crs: object, transform: object, los: np.ndarray
) -> pd.DataFrame:
    """Locate stations on the DISP grid and project MIDAS rates onto the LOS.

    `los` is the 3-band ``(east, north, up)`` LOS unit-vector stack of the
    DISP-S1-STATIC product on the same grid. Stations outside the frame or
    without a look are dropped. Adds ``row, col, los_{e,n,u},
    gps_hz_los_mm_yr, gps_los_mm_yr``.
    """
    import rasterio.transform
    from pyproj import Transformer

    to_xy = Transformer.from_crs("EPSG:4326", crs, always_xy=True)
    xs, ys = to_xy.transform(stations["lon"].to_numpy(), stations["lat"].to_numpy())
    rows, cols = rasterio.transform.rowcol(transform, xs, ys)
    st = stations.copy()
    st["row"], st["col"] = np.asarray(rows), np.asarray(cols)
    ny, nx = los.shape[1:]
    inside = (st["row"] >= 0) & (st["row"] < ny) & (st["col"] >= 0) & (st["col"] < nx)
    st = st[inside].copy()
    for k, b in zip("enu", range(3), strict=True):
        st[f"los_{k}"] = los[b, st["row"].to_numpy(), st["col"].to_numpy()]
    st = st[np.isfinite(st["los_u"]) & (st["los_u"] != 0)].copy()
    st["gps_hz_los_mm_yr"] = st["los_e"] * st["ve_mm_yr"] + st["los_n"] * st["vn_mm_yr"]
    st["gps_los_mm_yr"] = st["gps_hz_los_mm_yr"] + st["los_u"] * st["vu_mm_yr"]
    return st

Reports and CLI

Reports and traceability (plan T24).

Per frame: stations.csv, sill_draws.npz, SVG figures and report.html; across frames: release_summary.html. Every figure is followed by its data table (the table view the chart conventions require, and the only form a reviewer can quote). Colours: the validated categorical order blue (calibrated), orange (raw), aqua (baseline); status colours only for the verdict, always with a text label.

write_outputs(metrics, stations, dd, out)

stations.csv, sill_draws.npz, figures and report.html in out.

Source code in src/disp2vlm_validation/report.py

def write_outputs(
    metrics: dict[str, Any],
    stations: pd.DataFrame,
    dd: Any,
    out: Path,
) -> None:
    """``stations.csv``, ``sill_draws.npz``, figures and ``report.html`` in `out`."""
    out = Path(out)
    figs = out / "figures"
    figs.mkdir(parents=True, exist_ok=True)
    stations.to_csv(out / "stations.csv")
    np.savez_compressed(
        out / "sill_draws.npz", **{k: s.draws for k, s in dd.sills.items()}
    )
    sv_df = fig_semivariogram(
        dd.semivariograms, figs / "semivariogram.svg", metrics["sill_min_km"]
    )
    common = stations[np.all(np.isfinite(stations.filter(like="pix::")), axis=1)]
    st_df = fig_scatter(common, figs / "velocity_scatter.svg")
    by_class = pd.DataFrame(metrics["velocity_by_class"]).set_index("class")
    fig_class_rmse(by_class, figs / "class_rmse.svg")
    page = (
        _env()
        .get_template("frame.html.j2")
        .render(
            m=metrics,
            trace=trace_table(metrics),
            sv_table=_table_html(
                sv_df.pivot_table(
                    index="distance_km", columns="variant", values="gamma_mm2"
                )
                if not sv_df.empty
                else sv_df
            ),
            st_table=_table_html(st_df),
        )
    )
    (out / "report.html").write_text(page)

write_summary(run_dirs, out, title='Release summary')

release_summary.html across frame runs (their metrics.json).

Source code in src/disp2vlm_validation/report.py

def write_summary(
    run_dirs: Sequence[Path | str], out: Path | str, title: str = "Release summary"
) -> Path:
    """``release_summary.html`` across frame runs (their ``metrics.json``)."""
    out = Path(out)
    out.parent.mkdir(parents=True, exist_ok=True)
    runs = []
    for d in run_dirs:
        m = json.loads((Path(d) / "metrics.json").read_text())
        m["_link"] = os.path.relpath(Path(d) / "report.html", out.parent)
        runs.append(m)
    entries = load_traceability()
    trace_rows = [
        {
            "id": e["id"],
            "requirement": e["requirement"],
            "statuses": [TRACE_CHECKS[e["check"]](m)[0] for m in runs],
        }
        for e in entries
    ]
    page = (
        _env()
        .get_template("summary.html.j2")
        .render(
            title=title,
            runs=runs,
            n_pass=sum(m["verdict"] == "PASS" for m in runs),
            trace_rows=trace_rows,
        )
    )
    out.write_text(page)
    return out

load_traceability()

Entries of the packaged traceability.yaml.

Source code in src/disp2vlm_validation/report.py

def load_traceability() -> list[dict[str, Any]]:
    """Entries of the packaged ``traceability.yaml``."""
    text = (
        resources.files("disp2vlm_validation").joinpath("traceability.yaml").read_text()
    )
    return yaml.safe_load(text)

disp2vlm-validate command line.

gate prints the release gate and the benchmark frames; run validates one candidate on one frame (run.yaml, see pipeline.RunConfig); vlm validates a VLM run against GNSS (vlm.yaml, see vlm.VlmRunConfig); summary gathers frame runs into release_summary.html; demo writes a small synthetic frame and its configs to try the others on.

demo(out_dir, /, seed=1)

Write a small synthetic frame and ready-to-run configs (no download).

Then run, e.g., disp2vlm-validate run --config OUT_DIR/configs/perfect.yaml.

Source code in src/disp2vlm_validation/cli.py

def demo(out_dir: Path, /, seed: int = 1) -> None:
    """Write a small synthetic frame and ready-to-run configs (no download).

    Then run, e.g., ``disp2vlm-validate run --config OUT_DIR/configs/perfect.yaml``.
    """
    from .demo import write_demo

    for path in write_demo(out_dir, seed):
        print(f"wrote {path}")
    print(f"next: disp2vlm-validate run --config {out_dir}/configs/perfect.yaml")

gate(fmt='text')

Print the CalVal release gate and the benchmark frames (PRD 2.7).

Source code in src/disp2vlm_validation/cli.py

def gate(
    fmt: Annotated[Literal["text", "json"], tyro.conf.arg(name="format")] = "text",
) -> None:
    """Print the CalVal release gate and the benchmark frames (PRD 2.7)."""
    if fmt == "json":
        print(json.dumps({"gate": asdict(GATE), "frames": BENCHMARK_FRAMES}, indent=1))
        return
    print(f"disp2vlm-validation {__version__}: CalVal release gate")
    for key, value in asdict(GATE).items():
        print(f"  {key:30s} {value}")
    print("benchmark frames:")
    for num, f in BENCHMARK_FRAMES.items():
        print(f"  {num}. {f['frame']:7s} {f['name']:32s} {f['condition']}")

main()

Entry point.

Source code in src/disp2vlm_validation/cli.py

def main() -> None:
    """Entry point."""
    tyro.extras.subcommand_cli_from_dict(
        {"gate": gate, "run": run, "vlm": vlm, "summary": summary, "demo": demo},
        description=__doc__,
    )

run(config, no_cache=False, verbose=False)

Validate one candidate on one frame; write metrics, figures, report.html.

Parameters

config run.yaml (disp2vlm_validation.pipeline.RunConfig). no_cache Recompute even when metrics.json has the same config hash. verbose Debug logging.

Source code in src/disp2vlm_validation/cli.py

def run(config: Path, no_cache: bool = False, verbose: bool = False) -> None:
    """Validate one candidate on one frame; write metrics, figures, report.html.

    Parameters
    ----------
    config
        ``run.yaml`` (`disp2vlm_validation.pipeline.RunConfig`).
    no_cache
        Recompute even when ``metrics.json`` has the same config hash.
    verbose
        Debug logging.

    """
    from .pipeline import load_run_config, run_validation

    logging.basicConfig(
        level=logging.DEBUG if verbose else logging.INFO,
        format="%(asctime)s %(levelname)s %(message)s",
    )
    m = run_validation(load_run_config(config), use_cache=not no_cache)
    print(f"{m['frame']} {m['label']}: {m['verdict']}")
    for r in m["verdict_reasons"]:
        print(f"  - {r}")

summary(runs, /, out=Path('release_summary.html'), title='Release summary')

Gather frame runs (folders with metrics.json) into one HTML page.

Source code in src/disp2vlm_validation/cli.py

def summary(
    runs: list[Path],
    /,
    out: Path = Path("release_summary.html"),
    title: str = "Release summary",
) -> None:
    """Gather frame runs (folders with ``metrics.json``) into one HTML page."""
    from .report import write_summary

    path = write_summary(runs, out, title)
    for d in runs:
        m = json.loads((Path(d) / "metrics.json").read_text())
        print(f"{m['frame']:8s} {m['label']:24s} {m['verdict']}")
    print(f"wrote {path}")

vlm(config, verbose=False)

Validate one VLM run against independent GNSS (vlm.yaml).

See disp2vlm_validation.vlm.VlmRunConfig for the fields.

Source code in src/disp2vlm_validation/cli.py

def vlm(config: Path, verbose: bool = False) -> None:
    """Validate one VLM run against independent GNSS (``vlm.yaml``).

    See `disp2vlm_validation.vlm.VlmRunConfig` for the fields.
    """
    from .vlm import VlmRunConfig, run_vlm

    logging.basicConfig(
        level=logging.DEBUG if verbose else logging.INFO,
        format="%(asctime)s %(levelname)s %(message)s",
    )
    from .config_io import read_config

    cfg = VlmRunConfig(**read_config(config))
    m = run_vlm(cfg)
    s = m["vertical_rate"]["gated_both_looks"]
    print(
        f"{m['frame']} {m['label']}: {m['verdict']} (n={s.get('n')}, median"
        f" {s.get('median', float('nan')):+.2f}, NMAD"
        f" {s.get('nmad', float('nan')):.2f}, RMSE {s.get('rmse', float('nan')):.2f}"
        " mm/yr)"
    )
    for r in m["verdict_reasons"]:
        print(f"  - {r}")