ublox-gps-data-processing/experiments/vbs_iono.py

#!/usr/bin/env python3
"""
vbs_iono.py - Carrier-leveled slant ionosphere estimation + LIM cross-validation.

Phase 2a of the VBS exploration (gate G2a): measures whether multi-station
ionosphere interpolation can be accurate enough to help a virtual base station,
BEFORE building the application layer.

Per station/satellite arc: geometry-free carrier LG = lambda1*L1 - lambda2*L2
(clock-immune by construction) leveled to the geometry-free code PG = P2 - P1
(arc median), giving slant iono delay at L1 (metres; contains per-station DCB,
handled as a per-station bias in the fit).

LIM (Wanninger linear interpolation model): per satellite, per time bin, fit
I = a0 + a1*dN + a2*dE across stations, plus one shared per-station bias b_r
(gauge: b=0 at the first station). Leave-one-out: predict each held-out station
from the others; report residual scatter (after removing the held-out station's
own median offset, which is its unobservable DCB) vs its distance from the
network centroid.

Usage:
  vbs_iono.py --stations msin:X:Y:Z mary:X:Y:Z ... --dir data/204 --sp3 FILE
              [--bin-s 300] [--elev-mask 20]
"""

import argparse
import math
import sys
from collections import defaultdict
from datetime import datetime, timezone
from pathlib import Path

import numpy as np

sys.path.insert(0, str(Path(__file__).resolve().parent))
sys.path.insert(0, str(Path(__file__).resolve().parent.parent))
from vbs_synth import Sp3, parse_header, FREQ, glo_freq, glonass_slots  # noqa: E402

C = 299792458.0
F1, F2 = 1575.42e6, 1227.60e6
GAMMA = (F1 / F2) ** 2
L1M = C / F1
L2M = C / F2


def read_gf(path, slots):
    """Per (sat) -> time-sorted arrays (t_sec_of_day, LG_m, PG_m) for GPS+GLO."""
    lines = open(path, errors="replace").read().splitlines()
    types, xyz, hdr_end = parse_header(lines)
    nlps = (len(types) + 4) // 5
    idx = {t: k for k, t in enumerate(types)}
    need = all(t in idx for t in ("L1", "L2", "P2")) and ("P1" in idx or "C1" in idx)
    if not need:
        raise ValueError(f"{path}: missing L1/L2/P1/P2")
    i_l1, i_l2 = idx["L1"], idx["L2"]
    i_p1, i_c1 = idx.get("P1"), idx.get("C1")
    i_p2, i_c2 = idx.get("P2"), idx.get("C2")
    out = defaultdict(list)
    i = hdr_end + 1
    while i < len(lines):
        line = lines[i]
        if len(line) < 32 or not line[:26].strip():
            i += 1
            continue
        try:
            flag = int(line[26:29]); nsat = int(line[29:32])
        except ValueError:
            i += 1
            continue
        if flag > 1:
            i += nsat + 1
            continue
        p = line[:26].split()
        tsec = int(p[3]) * 3600 + int(p[4]) * 60 + float(p[5])
        sats = []
        nlin = (nsat + 11) // 12
        for k in range(nlin):
            seg = lines[i + k][32:68]
            for j in range(12):
                s = seg[j * 3:(j + 1) * 3]
                if s.strip():
                    sats.append(s.replace(" ", "0"))
        i += nlin
        for sat in sats:
            block = lines[i:i + nlps]
            i += nlps
            sysid = sat[0]
            if sysid == "G":
                l1f, l2f = F1, F2
            elif sysid == "R" and sat in slots:
                l1f, l2f = glo_freq("1", slots[sat]), glo_freq("2", slots[sat])
            else:
                continue

            def val(ti):
                ln = block[ti // 5] if ti // 5 < len(block) else ""
                seg = ln[(ti % 5) * 16:(ti % 5) * 16 + 14]
                return float(seg) if seg.strip() else None

            l1, l2 = val(i_l1), val(i_l2)
            # value-level P/C fallback: receivers fill P1 or C1 (and P2 or C2)
            # depending on satellite/tracking mode
            p1 = val(i_p1) if i_p1 is not None else None
            if p1 is None and i_c1 is not None:
                p1 = val(i_c1)
            p2 = val(i_p2) if i_p2 is not None else None
            if p2 is None and i_c2 is not None:
                p2 = val(i_c2)
            if None in (l1, l2, p1, p2):
                continue
            gamma = (l1f / l2f) ** 2
            lg = (C / l1f) * l1 - (C / l2f) * l2       # metres
            pg = p2 - p1
            # normalize to L1(GPS)-equivalent iono metres
            out[sat].append((tsec, lg / (gamma - 1) * (l1f / F1) ** 2,
                             pg / (gamma - 1) * (l1f / F1) ** 2))
    return {s: np.array(v) for s, v in out.items()}, np.array(xyz)


def leveled_stec(gf, gap_s=120.0, jump_m=0.5, min_arc=10):
    """Level carrier GF to code GF per arc -> (t, stec_L1_m) arrays per sat."""
    out = {}
    for sat, arr in gf.items():
        t, lg, pg = arr[:, 0], arr[:, 1], arr[:, 2]
        order = np.argsort(t)
        t, lg, pg = t[order], lg[order], pg[order]
        arcs = []
        start = 0
        for k in range(1, len(t)):
            if t[k] - t[k - 1] > gap_s or abs(lg[k] - lg[k - 1]) > jump_m:
                arcs.append((start, k))
                start = k
        arcs.append((start, len(t)))
        ts, ss = [], []
        for a, b in arcs:
            if b - a < min_arc:
                continue
            bias = np.median(lg[a:b] - pg[a:b])
            ts.append(t[a:b])
            ss.append(lg[a:b] - bias)
        if ts:
            out[sat] = (np.concatenate(ts), np.concatenate(ss))
    return out


def elevation(sp3, sat, tsec_gps_day0, xyz):
    pos = sp3.pos(sat, tsec_gps_day0)
    if pos is None:
        return -90.0
    d = pos - xyz
    up = xyz / np.linalg.norm(xyz)
    return math.degrees(math.asin(np.dot(d, up) / np.linalg.norm(d)))


def main():
    ap = argparse.ArgumentParser(description=__doc__,
                                 formatter_class=argparse.RawDescriptionHelpFormatter)
    ap.add_argument("--stations", nargs="+", required=True,
                    metavar="name:X:Y:Z", help="station obs basenames + ECEF")
    ap.add_argument("--dir", type=Path, required=True)
    ap.add_argument("--suffix", default="2040.26o")
    ap.add_argument("--sp3", type=Path, required=True)
    ap.add_argument("--nav", type=Path, required=True)
    ap.add_argument("--bin-s", type=float, default=300.0)
    ap.add_argument("--elev-mask", type=float, default=20.0)
    ap.add_argument("--gps-only", action="store_true",
                    help="exclude GLONASS (per-slot inter-channel code biases "
                         "break the single per-station bias model)")
    args = ap.parse_args()

    sp3 = Sp3(args.sp3)
    slots = glonass_slots(args.nav)
    stations = {}
    for spec in args.stations:
        name, x, y, z = spec.split(":")
        stations[name] = np.array([float(x), float(y), float(z)])

    # ENU offsets (km) about network centroid, for the LIM fit
    from process_gps import ecef_to_llh, llh_to_enu  # noqa: E402
    cen = np.mean(list(stations.values()), axis=0)
    cen_llh = ecef_to_llh(*cen)
    offs = {}
    for name, xyz in stations.items():
        llh = ecef_to_llh(*xyz)
        dn, de, _ = llh_to_enu(llh[0], llh[1], llh[2], *cen_llh)
        offs[name] = (float(dn) / 1e3, float(de) / 1e3)   # km N, km E

    print("[vbs_iono] loading + leveling", flush=True)
    stec = {}
    for name in stations:
        gf, _ = read_gf(args.dir / f"{name}{args.suffix}", slots)
        if args.gps_only:
            gf = {s: v for s, v in gf.items() if s.startswith("G")}
        stec[name] = leveled_stec(gf)
        n = sum(len(v[0]) for v in stec[name].values())
        print(f"  {name}: {len(stec[name])} sats, {n} leveled samples")

    # time-bin medians: station -> sat -> {bin: stec}
    binned = {}
    for name, sats in stec.items():
        b = {}
        for sat, (t, s) in sats.items():
            k = (t // args.bin_s).astype(int)
            d = {}
            for kk in np.unique(k):
                d[int(kk)] = float(np.median(s[k == kk]))
            b[sat] = d
        binned[name] = b

    names = list(stations)

    # Stage 1: per-station bias (DCB) as a DAILY constant via median polish.
    # Leaving these as free per-bin parameters makes the 4-station LIM fit
    # ill-conditioned (bias ~ degenerate with a satellite-common gradient) and
    # extrapolation explodes; hardware DCBs are stable over a day, so estimate
    # them once with full-day redundancy.
    bias = {n: 0.0 for n in names}
    for _ in range(4):
        net_med = {}
        for n in names:
            for sat, d in binned[n].items():
                for k, v in d.items():
                    net_med.setdefault((sat, k), []).append(v - bias[n])
        net_med = {sk: float(np.median(v)) for sk, v in net_med.items() if len(v) >= 3}
        for n in names:
            diffs = [binned[n][sat][k] - m for (sat, k), m in net_med.items()
                     if sat in binned[n] and k in binned[n][sat]]
            if diffs:
                bias[n] = float(np.median(diffs))
    print("  station daily biases (m):",
          {n: round(b, 2) for n, b in bias.items()})
    for n in names:
        for sat in binned[n]:
            for k in binned[n][sat]:
                binned[n][sat][k] -= bias[n]

    print(f"\n[vbs_iono] LOOCV over {len(names)} stations, bin={args.bin_s:.0f}s, "
          f"elev mask {args.elev_mask} deg")
    results = {}
    for held in names:
        fit_names = [n for n in names if n != held]
        residuals = []
        allbins = set()
        for sat in binned[held]:
            allbins.update(binned[held][sat].keys())
        for bin_k in sorted(allbins):
            t_mid = (bin_k + 0.5) * args.bin_s
            # collect sats present at held-out + >=3 fit stations, above mask
            rows, sats_ok = [], []
            for sat in binned[held]:
                if bin_k not in binned[held][sat]:
                    continue
                have = [n for n in fit_names
                        if sat in binned[n] and bin_k in binned[n][sat]]
                if len(have) < 3:
                    continue
                if elevation(sp3, sat, t_mid, stations[held]) < args.elev_mask:
                    continue
                sats_ok.append((sat, have))
            if len(sats_ok) < 3:
                continue
            # per-satellite LIM planes; station biases already removed (stage 1)
            dn, de = offs[held]
            for sat, have in sats_ok:
                A = np.array([(1.0, offs[n][0], offs[n][1]) for n in have])
                y = np.array([binned[n][sat][bin_k] for n in have])
                sol, *_ = np.linalg.lstsq(A, y, rcond=None)
                pred = sol[0] + sol[1] * dn + sol[2] * de
                residuals.append(pred - binned[held][sat][bin_k])
        residuals = np.array(residuals)
        if len(residuals) == 0:
            continue
        # remove held-out station's unobservable DCB offset
        residuals = residuals - np.median(residuals)
        dist_cen = math.hypot(*offs[held])
        results[held] = (dist_cen, residuals)
        print(f"  hold {held}: n={len(residuals):5d}  centroid dist {dist_cen:5.1f} km  "
              f"median|r|={np.median(np.abs(residuals)):.3f} m  "
              f"95%|r|={np.percentile(np.abs(residuals), 95):.3f} m")

    med = {n: float(np.median(np.abs(r))) for n, (d, r) in results.items()}
    dists = {n: d for n, (d, r) in results.items()}
    if len(results) >= 3:
        x = np.array([dists[n] for n in results])
        yv = np.array([med[n] for n in results])
        slope = np.polyfit(x, yv, 1)[0] * 10  # m per 10 km
        print(f"\n  residual-vs-distance slope: {slope*100:.1f} cm per 10 km")
    overall = np.concatenate([r for _, r in results.values()])
    print(f"  OVERALL: median|r|={np.median(np.abs(overall)):.3f} m, "
          f"95%|r|={np.percentile(np.abs(overall), 95):.3f} m  "
          f"(gate G2a: median <= ~0.10 m)")


if __name__ == "__main__":
    main()