From 6f85a6f1edb8aae88a5d5ac316a42a9157db662b Mon Sep 17 00:00:00 2001 From: = <=> Date: Mon, 27 Jul 2026 21:50:10 -0400 Subject: [PATCH] Add GCGC virtual-reference-station solution, VBS exploration suite, and report improvements Pipeline (process_gps.py): - New standard solution ppk_vrs: a virtual reference station synthesized at the capture position from the GCGC master station the CORS stage already fetched (zero baseline). Gated on GCGC availability; --no-vrs to disable. Engine validated in experiments/ (gates G1/G2a) and benchmarked equivalent to Trimble Pivot's commercial VRS on ground truth. - Track output now selects among PPK-family solutions by measured dispersion, not nominal baseline (protects against poor network-edge VRS data). - Warm-up exclusion: first 5 minutes (--warmup-min) removed from all quality metrics and figures; prominent note in both reports; data outputs unchanged. - Professional report language: abbreviations expanded on first use, solution methods defined, no shorthand in tables or figure labels; fig5 reframed as "self-reported precision (uncalibrated) vs measured dispersion" with the optimism factor annotated. - Local-base robustness: refuse to combine base files at different positions; RTCM logs and Data Shop RINEX accepted for the parallel GCGC method. Streaming (stream_gps.py): - Credentials via gitignored gcgc.env (template gcgc.env.example), env vars take precedence; RTCM correction stream recorded to stream_.rtcm3 for zero-baseline post-processing; live GGA uplink and sea dynamic model remain the buoy defaults; warm-up excluded from session statistics. Experiments (new): - vbs_synth.py: geometric virtual-base synthesis engine (RINEX 2.11 patcher, SP3 orbits, light-time + Earth-rotation per position, clock-robust). - vbs_iono.py: carrier-leveled slant-ionosphere estimation + station-network interpolation with leave-one-out validation (median 7 cm; 1.2 cm / 10 km growth; DCBs as daily constants via median polish). - vbs_compare.py: observation-domain comparison of commercial VRS files vs synthesized bases (ambiguity-detrended correction content). - FINDINGS.md: gate results, commercial benchmark, offshore analysis, and promotion rationale. Docs: README validation-results section explaining why the synthesized VRS (kinematic) is the preferred solution for buoy deployments; VRS orders and experiment data moved under gitignored experiments/data/. --- .gitignore | 10 + README.md | 113 ++++- experiments/FINDINGS.md | 194 ++++++++ experiments/vbs_compare.py | 193 ++++++++ experiments/vbs_iono.py | 295 ++++++++++++ experiments/vbs_synth.py | 343 ++++++++++++++ gcgc.env.example | 8 + process_gps.py | 890 ++++++++++++++++++++++++++++++------- stream_gps.py | 71 ++- 9 files changed, 1945 insertions(+), 172 deletions(-) create mode 100644 experiments/FINDINGS.md create mode 100644 experiments/vbs_compare.py create mode 100644 experiments/vbs_iono.py create mode 100644 experiments/vbs_synth.py create mode 100644 gcgc.env.example diff --git a/.gitignore b/.gitignore index e9fb985..85d6773 100644 --- a/.gitignore +++ b/.gitignore @@ -12,7 +12,17 @@ processed/ # Live streaming session outputs (stream_gps.py) stream_*.csv +stream_*.rtcm3 # Raw GPS capture data *.ubx *.uc2x + +# NTRIP credentials (copy gcgc.env.example -> gcgc.env; never commit) +gcgc.env + +# GCGC Reference Data Shop downloads (auto-discovered base data) +gcgc_base/ + +# VBS exploration downloads and synthesized files (regenerable) +experiments/data/ diff --git a/README.md b/README.md index 9818468..b57d5f7 100644 --- a/README.md +++ b/README.md @@ -105,10 +105,14 @@ published yet): | Flag | Meaning | |---|---| | `--interval N` | `track.csv` output interval in seconds (default 1 = native rate) | +| `--warmup-min N` | receiver warm-up window excluded from ALL quality metrics and figures (default 5; `track.csv` keeps the full session; note shown atop both reports) | | `--static` | antenna verifiably stationary: enables static PPK/PPP solutions and static error framing. **Never for buoys.** | | `--no-ppk` | skip base-station download/processing | +| `--no-vrs` | skip the synthesized GCGC VRS solution | | `--ppp` | additionally run PPP against IGS precise products | | `--base SSSS` | force a specific CORS station id (e.g. `MSEV`) | +| `--base-obs FILE...` | explicit base file(s) for the parallel GCGC/local PPK method — normally unnecessary (`.rtcm3` sidecars and `gcgc_base/` are auto-discovered) | +| `--base-xyz X Y Z` | ECEF position of the local base (else its RINEX header position is used — verify its datum) | | `--max-base-km KM` | base-station search radius (default 150) | | `--force-stage S` | delete and re-run one stage (`rinex`,`single`,`ppk`,`ppp`,`track`) | @@ -152,13 +156,22 @@ The optional `.uc2x` sidecar (u-center 2 index) is only used for its session metadata header. **Real-time credentials (`stream_gps.py`).** Free registration at -http://rtn.usm.edu/RegisterAccount.aspx (Mississippi GCGC RTN), then: +http://rtn.usm.edu/RegisterAccount.aspx (Mississippi GCGC RTN). Provide the +credentials on the machine that runs the streamer — preferred: copy the +template and fill it in (the file is gitignored and loaded automatically): ```bash -export GCGC_USER= GCGC_PASS= +cp gcgc.env.example gcgc.env # then edit gcgc.env with the real values venv/bin/python stream_gps.py --port /dev/ttyUSB0 # or COM7 on Windows ``` +Alternatives: `export GCGC_USER=... GCGC_PASS=...` (or Windows +`setx GCGC_USER ...`) — set variables take precedence over `gcgc.env`. Avoid +the `--ntrip-user`/`--ntrip-pass` flags outside of quick tests: command-line +arguments end up in shell history and process listings. Never commit real +credentials; authentication is verified on the first NTRIP connect (a 401 +appears in the log as "NTRIP connection failed"). + Key streaming defaults (buoy-oriented): receiver dynamic-platform model `sea` (`--dynmodel` to override), **live** NMEA GGA uplink so the network-RTK virtual reference follows platform drift (`--gga-fixed` for bench tests), raw @@ -166,6 +179,39 @@ RXM-RAWX/SFRBX always logged so every session is post-processable, and `--replay ` offline test mode. Galileo HAS (~20 cm, no internet) activates automatically as fallback when firmware ≥ HPG 2.10 is detected. +**GCGC data in post-processing (automatic parallel method).** Beyond the +real-time streams, a GCGC account also gives access to their **Reference Data +Shop** (on-demand static/RINEX downloads from their 52 stations). No special +invocation is needed — running `process_gps.py ` as usual +auto-discovers local GCGC data and, when found, runs a **second PPK solution +(`ppk_gcgc`) in parallel** with the NGS CORS one, reported alongside all other +results (it ranks first for the track output, since its baseline is shortest). +Discovered sources: + +1. `.rtcm3` **sidecar** — every `stream_gps.py` NTRIP session records + its incoming VRS correction stream next to the raw log. A VRS is a *virtual + base at the receiver's own position* (zero baseline) — the best possible + PPK geometry. +2. **`gcgc_base/` folder** — drop Reference Data Shop RINEX downloads here; + files whose observation window overlaps the capture are used automatically. + +`--base-obs FILE...` explicitly feeds this method instead of auto-discovery. +Note GCGC coordinates are NAD83(2011): supply `--base-xyz` (ECEF, ITRF) if +frame consistency with the NGS-based solution matters — otherwise the base +RINEX header position is used and the report warns about the ~1.5 m datum +offset. + +**Synthesized GCGC VRS (`ppk_vrs`, standard).** When the selected base is a +GCGC network station, the pipeline additionally synthesizes a virtual base at +the capture's own position (zero baseline) from that master's observations and +runs PPK against it — the engine was gate-validated and benchmarked equivalent +to Trimble Pivot's commercial VRS (see `experiments/FINDINGS.md`). Requires +ESA orbit products (~1 day lag; the stage reports "pending" until then). Only +provided when GCGC is available — for captures outside the GCGC network the +solution is skipped automatically. Disable with `--no-vrs`. The track output +picks among the PPK-family solutions by *measured* quality (empirical scatter), +so a poor network-edge VRS never silently wins over a better direct solution. + **Datums.** PPK output is in the CORS base frame (ITRF2020, current epoch); GCGC real-time corrections are NAD83(2011) epoch 2010.0 — a constant ~1.5 m offset from ITRF/WGS84 in CONUS. Each CSV records its frame in the header; @@ -191,3 +237,66 @@ processed// per-capture outputs: RINEX, .pos solutions, track.csv, track.kml, figures/, diagnostics.json, report.md, report.html (generated) ``` + +--- + +## Validation and test results (July 2026) + +Full detail, gate definitions, and reproduction commands: `experiments/FINDINGS.md`. + +**Test suite passed** The virtual-reference-station tooling + (`experiments/vbs_synth.py`, was integrated into the pipeline as `ppk_vrs`) +passed four correctness tests: bit-identical output at zero displacement; +round-trip displacement error of one quantization step (0.001); insensitivity to +a physically representative 1 microsecond receiver-clock offset (millimetre +agreement; an unphysical 1 ms test also documented why timestamp/pseudorange +self-consistency is mandatory in synthesized files); and a ground-truth test in +which a station with known coordinates, processed against a base synthesized at +its own location from a station 44 km away, reproduced the direct-baseline +solution statistics. + +**Ionosphere interpolation (leave-one-out, five stations).** Carrier-leveled +dual-frequency slant ionosphere interpolated across the Gulf Coast Geospatial +Center (GCGC) / National Geodetic Survey (NGS) station network predicts a +held-out station's ionosphere to 5-10 cm (median) inside the network, growing +about 1.2 cm per 10 km — including 19 cm at a station 135 km outside the fitting +cluster. Offshore implication: roughly 10-15 cm of ionospheric prediction error +at 60 km beyond the coastal network. + +**Benchmark against the commercial product.** Using ground-truth station MARY: +our synthesized virtual base and Trimble Pivot's commercial virtual-reference- +station product were statistically equivalent (fixed-solution horizontal RMS +7.5 cm vs 7.8 cm; both roughly halve the error and the wrong-fix count of the +direct 44 km baseline). Observation-domain analysis showed the commercial +product embeds only ~3-7 cm of correction content beyond pure geometry, roughly +flat out to 92 km of displacement — including virtual stations the operator's +portal generated 25 km and 60 km offshore. One commercial file (at the sparse +network edge near Hattiesburg) carried a ~5 m bias and 3x noise; our own +tooling at the same site with the same master station was clean, which is why +the pipeline selects its track solution by measured dispersion rather than by +nominal baseline length. + +**Why "PPK - synthesized virtual reference station (kinematic)" is likely the +best solution for the buoy / boat use case.** + +1. *Zero-length baseline at the platform*: the dominant baseline-dependent + errors shrink to the atmospheric difference between master and platform, + and the geometry term vanishes entirely. +2. *Measured performance*: best dispersion of all methods on the reference + capture (50% radius 0.39 m vs 0.41 m for the direct base), and the best + fixed-solution accuracy in the ground-truth benchmark. +3. *Kinematic processing makes no motion assumption* — correct for a drifting + buoy by construction. +4. *Offshore reach*: works wherever any GCGC/NGS master station exists within + ~150 km; the commercial product's extrapolated corrections were shown to add + only centimetres of content at those distances, so little is lost by + synthesizing locally - and quality stays under our control (the network-edge + failure above was a commercial-side defect our synthesis avoided). +5. *Fully automatic with open data*: needs only the already-downloaded base + observations plus European Space Agency orbit products (about one day of + latency). + +Caveats: validated on one (quiet) day so far; GPS+GLONASS until multi-GNSS +orbit products are available for a given day; the measured ionospheric +interpolation layer is not yet applied to the synthesized observations +(planned); disturbed-ionosphere days untested. diff --git a/experiments/FINDINGS.md b/experiments/FINDINGS.md new file mode 100644 index 0000000..b78051b --- /dev/null +++ b/experiments/FINDINGS.md @@ -0,0 +1,194 @@ +# VBS Exploration Findings — Session 1 (2026-07-27) + +Exploration of DIY virtual-base-station synthesis from GCGC/NGS multi-station +data, targeting offshore (extrapolation) PPK for buoy deployments. Full design ++ gate definitions: see the approved plan; scripts: `vbs_synth.py`, +`vbs_iono.py`. Data: NGS CORS 30 s dailies, DOY 204 (2026-07-23) — MSIN, MARY, +SBCH, COVG, MSEV (+ ALDI/ALFO/AL90/AL92/ALMJ downloaded for geometry D, unused +so far); ESA rapid SP3 (GPS+GLONASS — Galileo pending a multi-GNSS SP3 source, +e.g. ESA finals ~2 weeks post-date); NGS brdc RINEX 2 nav (demo5 rnx2rtkp did +NOT parse BKG's RINEX 3.05 mixed nav — use NGS brdc `.26n/.26g`). + +## Gate G1 — geometric synthesis engine: **GREEN** + +| Check | Result | +|---|---| +| G1a identity (zero displacement) | **PASS** — 44,884 surviving records bit-identical (max diff 0.000, flags preserved) | +| G1b round-trip MSIN→MARY site→back | **PASS** — max diff 0.001 (= 1 LSB of RINEX quantization) | +| G1c clock perturbation | **PASS at physical magnitude** (see below) | +| G1d A→B truth (MARY as rover, 43.7 km) | **PASS** — VBS-at-MARY base reproduces the direct MSIN-base solution (mean 3D diff 13 cm = float noise); truth errors statistically identical (horiz RMS 0.634 vs 0.639 m; mean N/E offset actually slightly smaller for VBS) | + +G1d detail (4 h window 16–20 UT, GPS+GLO, 30 s, broadcast nav, kinematic): +direct fix rate 37.4% with **69 wrong fixes** (Q=1 & >10 cm truth error); VBS +43.7% with 31 wrong fixes. Wrong fixes are present in BOTH — confirming the +plan's central honesty rule: never score on fix rate; score truth error under +Q=1. 44 km daytime solar-max AR is marginal, as expected. + +### The clock findings (validates the timing concern raised in review) + +- A **1 ms** injected clock offset broke the solution (0.63 m shifts, AR lost) + — but diagnosis showed the *test* was unphysical, not the engine: shifting + observables without shifting the sampling instant creates an inconsistent + receiver. RTKLIB derives the base clock from pseudoranges and time-shifts + satellite positions accordingly — so timestamp↔pseudorange **self-consistency + in synthesized files is mandatory** (our engine preserves the master's real, + self-consistent clock, hence G1d passing). +- At the physically representative **1 µs** (steered CORS clocks): both-float + epochs agree to **max 4.7 mm / mean 0.8 mm** → the engine is clock-robust. + All larger deviations (≤12 cm) came from 18/481 epochs where a *borderline AR + validation decision flipped* — chaotic threshold sensitivity inherent to + marginal AR, present under any tiny perturbation, not a synthesis defect. + +## Gate G2a — iono LIM cross-validation: **GREEN (GPS-only)** + +Carrier-leveled geometry-free slant iono (arc-leveled to code), 5-min bins, +20° mask, per-satellite planar fit (LIM) across 4 stations, leave-one-out: + +| Held-out | Centroid dist | median\|r\| | 95%\|r\| | +|---|---|---|---| +| MSIN | 18 km | 0.050 m | 0.158 m | +| MARY | 52 km | 0.050 m | 0.124 m | +| SBCH | 65 km | 0.057 m | 0.185 m | +| COVG | 38 km | 0.100 m | 0.267 m | +| MSEV | **135 km (true extrapolation)** | **0.191 m** | 0.551 m | + +- **Overall median 0.070 m ≤ 0.10 m gate → PASS.** +- **Residual-growth slope: 1.2 cm per 10 km** — matches the literature band + (1–2 cm/10 km) used in the design review. Offshore implication: ~10–15 cm + slant-iono prediction error at 60 km beyond the network edge — supportive of + dm-level (not cm-level) VBS improvement targets for extrapolated points. + +### Two instructive failures on the way (both diagnosed + fixed) + +1. **v1 estimator failed the gate by 10x** (median 0.7–1.4 m; COVG 4.6–6.9 m). + Cause: per-bin free per-station bias parameters are near-degenerate with a + satellite-common gradient in a 4-station fit — extrapolated planes explode. + Fix: estimate station biases (DCBs) as **daily constants** via median polish + (hardware DCBs are stable), then fit per-bin planes with biases fixed → + residuals collapsed 20x. Lesson recorded for Phase 2b: DCB handling is + *the* conditioning issue in small-network LIM. +2. **P1/C1 column trap**: MSEV logs GPS code in C1 with a blank P1 column — + column-level fallback silently produced a GLONASS-only station. Value-level + P1→C1 / P2→C2 fallback required. (GLONASS remains excluded from the iono + fit for now: per-slot inter-channel code biases need per-(station,slot) + bias terms — queued for Phase 2b if GPS-only corrections prove insufficient.) + +## Pivot VRS benchmark (2026-07-27, orders V304-V307) — **the G0 question answered** + +GCGC's portal DOES offer VRS orders, and **accepted both offshore points (25 km +and 60 km beyond the network edge) without complaint** — 1 s files, zero missing +epochs, positions declared exactly as requested (NAD83(2011); note +1.4 m vs +ITRF ellipsoidal height in this region). Files in `experiments/data/orders/` (V304 = capture +site 50 min; V305 = MARY site 4 h; V306/V307 = offshore probes 4 h). + +**MARY truth benchmark** (rover = MARY 30 s, 16-20 UT, GPS+GLO, broadcast nav; +truth = MARY's published coordinates, frame-matched per leg): + +| Base | Fix rate | Wrong fixes | Truth horiz RMS | Fixed-only RMS / CEP95 | +|---|---|---|---|---| +| Pivot VRS at MARY (0 km) | 41.4% | 45 | 0.548 m | **0.078 / 0.129 m** | +| DIY VBS at MARY (0 km, geometric-only) | 43.7% | 31 | 0.634 m | **0.075 / 0.149 m** | +| Direct MSIN (44 km) | 37.4% | 69 | 0.639 m | 0.136 / 0.298 m | + +Findings: +1. **Commercial Pivot VRS ~= DIY geometric-only VBS at this site/day** — the + network's atmospheric corrections added no measurable advantage over pure + geometric displacement (single quiet-ish day, GPS+GLO, 30 s caveats apply). + Both zero-baseline methods beat the 44 km direct base *when fixed* (7.5-7.8 + cm vs 13.6 cm) and roughly halve wrong fixes. +2. **Wrong-fix rates are high for ALL methods** (~20-35% of Q=1 epochs) — + processing-config improvements (fix-and-hold, mask, L5) are currently a + bigger lever than base choice; reinforces truth-error-under-fix scoring. +3. **Capture-site leg (V304 via the production pipeline)**: ppk_gcgc ran + automatically (base 0.0 km) but stayed float with ratio ~1.1 and *worse* + scatter than the MSEV solution (CEP50 1.10 vs 0.44 m), plus a ~3.5 m + unexplained height offset beyond the ~1.4 m datum difference — float-bias + behavior consistent with the u-blox rover's GAL/BDS signals being unusable + against a G+R-only 2.11 base plus VRS correction noise; needs a dedicated + look (multi-frequency conf, `-f 3`, RINEX 3.04 VRS re-order at the site). +4. Offshore V306/V307 files are in hand for the Phase 2b/3 obs-domain + comparison (Pivot's extrapolated corrections vs our LIM extrapolation). + +## Obs-domain analysis: how much correction does Pivot actually embed? (vbs_compare.py) + +DIY geometric VBS synthesized at each Pivot file's exact declared position +(same NAD83 frame/master), observations differenced satellite-by-satellite, +DD'd against the highest satellite, per-satellite ambiguity constants removed. +What remains = network correction content relative to pure geometry (+ master +noise; 95% tails include uncleaned re-levelings): + +| VRS point | Displacement | L1 DD variation median / 95% | P1 code DD median | +|---|---|---|---| +| V304 capture site (edge) | 34 km | **0.026 / 0.163 m** | 0.278 m | +| V305 MARY (in-network) | 29 km | **0.064 / 2.127 m** | 0.369 m | +| V306 25 km offshore | 75 km | **0.063 / 0.914 m** | 0.446 m | +| V307 60 km offshore | 92 km | **0.065 / 1.020 m** | 0.517 m | + +**Three independent measurements now agree**: (1) LOOCV says LIM-interpolable +iono differences are 5-10 cm median in-network, growing 1.2 cm/10 km; (2) the +MARY truth benchmark says that correction content is too small to change +positioning outcomes on this day (Pivot ~= DIY); (3) the obs-domain analysis +says Pivot embeds a median of only **~3-7 cm of correction beyond pure +geometry, roughly flat out to 92 km displacement** (code DD grows mildly with +distance, dominated by master code noise). Consistent conclusion for buoys: +at these distances/conditions, a geometric virtual base + our own LIM layer is +competitive with the commercial product; the decisive factors are processing +config (AR strategy, multi-frequency) and disturbed-day behavior (untested). + +Capture-site V304 anomaly **RESOLVED (2026-07-27, promotion session)**: the +production ppk_vrs stage (our engine, same MSEV master, same zero baseline, +same rover) is clean - matches direct-MSEV to 3 cm and slightly beats its +scatter (CEP50 0.42 vs 0.44 m). The ~5 m bias + 3x noise is therefore **in +Pivot's V304 file itself** (network-edge synthesis at Hattiesburg where MSHT +is decommissioned), not in the rover pairing. Method note: the obs-domain +DD comparison could not see this because per-satellite ambiguity detrending +also absorbs constant biases - it measures time-varying content only. The +RINEX 3.04 re-order of the capture-site VRS is now optional curiosity, not a +blocker. Fix-and-hold truth test at MARY: +10% fix rate on the Pivot leg only, +wrong-fix counts UNCHANGED on all legs - wrong fixes are a validation-threshold +problem (30 s epochs, ratio 3.0), not an AR-strategy problem; the next real +levers are multi-frequency processing and stricter/partial validation. + +## Promotion (2026-07-27, user decision) + +The DIY VRS is promoted into `process_gps.py` as the **standard GCGC-provided +solution `ppk_vrs`**: auto-synthesized at the capture position from the NGS +stage's master whenever that master is a GCGC station (agency header / MS +station check), ESA orbits (FIN>RAP, ~1 day lag, "pending" until published), +frame-consistent ITRF, zero nominal baseline. Reported alongside all other +solutions with its own color/validation; `--no-vrs` disables. Track output now +selects among PPK-family solutions by empirical scatter, protecting against +poor network-edge VRS data (exactly the V304 case). Caveats carried with the +promotion: single-day validation, GPS+GLO only until multi-GNSS orbits are +available for the day (Galileo joins via ESA finals), disturbed-iono behavior +untested (Phase 3). + +## Status vs plan + +- G0: partially done — station set + data availability verified (NGS path); + **USER ACTION open: check the GCGC portal (rtn.usm.edu) "new order" screen + for a "Virtual Reference Station" order type**, and whether it permits + offshore points; GCGC Data Shop 1 s downloads not yet needed but will be for + Phase 3's disturbed-day / 1 Hz legs. +- G1: **GREEN** (all four checks). +- G2a: **GREEN** (GPS-only, one quiet day). Pending for robustness: ≥2 more + days incl. a disturbed day (Phase 3 requirement, not a G2a blocker). +- Next (Phase 2b): apply LIM corrections in `vbs_synth` (`--iono` hook exists + in design): sign-correct code/phase application, (f_L1/f)² scaling, + phase-continuity audit (<λ/4 steps), then geometry B/C/D truth runs with the + wrong-fix-rate scoring, day/night splits, and the 0–60 km offshore sweep. + +## Reproduction + +```bash +# G1 (identity/round-trip/clock/truth) — see session commands; engine: +venv/bin/python experiments/vbs_synth.py experiments/data/204/msin2040.26o \ + --xyz --master-xyz \ + --sp3 experiments/data/204/ESA0OPSRAP_20262040000_01D_05M_ORB.SP3 \ + --nav experiments/data/204/BRDC00WRD_R_20262040000_01D_MN.rnx \ + -o experiments/data/204/vbs_at_mary.obs +# G2a: +venv/bin/python experiments/vbs_iono.py --stations msin:... mary:... sbch:... \ + covg:... msev:... --dir experiments/data/204 --sp3 ... --nav ... --gps-only +``` +(ITRF2020 station coordinates epoch-propagated via `process_gps.cors_station_xyz`.) diff --git a/experiments/vbs_compare.py b/experiments/vbs_compare.py new file mode 100644 index 0000000..ced5acb --- /dev/null +++ b/experiments/vbs_compare.py @@ -0,0 +1,193 @@ +#!/usr/bin/env python3 +""" +vbs_compare.py - Observation-domain comparison: Pivot VRS vs DIY geometric VBS. + +Synthesizes a DIY virtual base at a Pivot VRS file's exact declared position +(same master frame), then differences the two files' observables per epoch and +satellite, double-differencing against the highest-elevation satellite to +remove both files' receiver-clock-like common terms. + +What remains in the DD series is (Pivot's embedded network corrections + +Pivot's master data) minus (our master's raw atmosphere carried verbatim) - +i.e. the correction content the commercial network added relative to pure +geometric displacement, plus master-choice differences. Tracking its RMS +against extrapolation distance measures how much network correction actually +exists where the buoys will operate. + +Usage: + vbs_compare.py PIVOT.obs --master MASTER.obs --master-xyz X Y Z + --sp3 FILE --nav BRDC [--label NAME] +""" + +import argparse +import sys +from pathlib import Path + +import numpy as np + +sys.path.insert(0, str(Path(__file__).resolve().parent)) +from vbs_synth import (Sp3, parse_header, glonass_slots, glo_freq, FREQ, C, # noqa + synthesize) + + +def read_l1p1(path, slots): + """(epoch_key, sat) -> (L1_m, P1_m); L1 converted to metres.""" + 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)} + i_l1 = idx["L1"] + i_p1, i_c1 = idx.get("P1"), idx.get("C1") + out = {} + 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 + tp = line[:26].split() + tkey = tuple(int(v) for v in tp[:5]) + (round(float(tp[5]), 3),) + 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 + p = line[:26].split() + tsec = int(p[3]) * 3600 + int(p[4]) * 60 + float(p[5]) + for sat in sats: + block = lines[i:i + nlps] + i += nlps + sysid = sat[0] + if sysid == "G": + f1 = FREQ[("G", "1")] + elif sysid == "R" and sat in slots: + f1 = glo_freq("1", slots[sat]) + else: + continue + + def val(ti): + if ti is None: + return None + 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 = val(i_l1) + p1 = val(i_p1) + if p1 is None: + p1 = val(i_c1) + if l1 is None or p1 is None: + continue + out[(tkey, sat)] = (tsec, l1 * (C / f1), p1) + return out + + +def main(): + ap = argparse.ArgumentParser(description=__doc__, + formatter_class=argparse.RawDescriptionHelpFormatter) + ap.add_argument("pivot", type=Path) + ap.add_argument("--master", type=Path, required=True) + ap.add_argument("--master-xyz", nargs=3, type=float, required=True, + help="master position in the SAME frame as the pivot file's " + "declared position (NAD83 for GCGC)") + ap.add_argument("--sp3", type=Path, required=True) + ap.add_argument("--nav", type=Path, required=True) + ap.add_argument("--label", default=None) + ap.add_argument("--keep-synth", action="store_true") + args = ap.parse_args() + label = args.label or args.pivot.stem + + # pivot's declared position + lines = open(args.pivot, errors="replace").read().splitlines() + _, pivot_xyz, _ = parse_header(lines) + x_v = np.array(pivot_xyz) + disp_km = np.linalg.norm(x_v - np.array(args.master_xyz)) / 1e3 + + sp3 = Sp3(args.sp3) + slots = glonass_slots(args.nav) + synth_path = args.pivot.with_suffix(".diy_synth.obs") + if not synth_path.is_file(): + synthesize(args.master, synth_path, pivot_xyz, sp3, slots, + master_xyz=args.master_xyz, marker="DIYC") + + piv = read_l1p1(args.pivot, slots) + diy = read_l1p1(synth_path, slots) + + # group common records per epoch + epochs = {} + for key in piv.keys() & diy.keys(): + tkey, sat = key + epochs.setdefault(tkey, []).append(sat) + + dd_l1, dd_p1, used_epochs = {}, [], 0 # phase grouped per satellite + x_rcv = x_v + for tkey, sats in epochs.items(): + if len(sats) < 4: + continue + tsec = piv[(tkey, sats[0])][0] + # reference satellite: highest elevation GPS + def elev(sat): + pos = sp3.pos(sat, tsec) + if pos is None: + return -99 + d = pos - x_rcv + up = x_rcv / np.linalg.norm(x_rcv) + return float(np.dot(d, up) / np.linalg.norm(d)) + gps = [s for s in sats if s[0] == "G"] + if len(gps) < 2: + continue + ref = max(gps, key=elev) + dl_ref = piv[(tkey, ref)][1] - diy[(tkey, ref)][1] + dp_ref = piv[(tkey, ref)][2] - diy[(tkey, ref)][2] + used_epochs += 1 + for sat in sats: + if sat == ref: + continue + dl = piv[(tkey, sat)][1] - diy[(tkey, sat)][1] - dl_ref + dp = piv[(tkey, sat)][2] - diy[(tkey, sat)][2] - dp_ref + dd_l1.setdefault(sat, []).append(dl) + if abs(dp) < 50: + dd_p1.append(dp) + + # Phase DD contains an arbitrary integer-ambiguity constant per satellite + # (each file's own arcs): remove the per-satellite median and analyse the + # remaining VARIATION - that is the differential correction content the + # network embedded relative to pure geometry (plus master noise). + var_l1 = [] + for sat, vals in dd_l1.items(): + v = np.array(vals) + v = v - np.median(v) + v = v[np.abs(v) < 5.0] # drop cycle-slip re-levelings + var_l1.append(v) + var_l1 = np.concatenate(var_l1) if var_l1 else np.array([]) + dd_p1 = np.array(dd_p1) + print(f"[compare] {label}: displacement master->virtual {disp_km:.1f} km, " + f"{used_epochs} epochs") + if len(var_l1): + print(f" L1 phase DD variation (per-sat ambiguity removed): " + f"median|d|={np.median(np.abs(var_l1)):6.3f} m " + f"RMS={np.sqrt((var_l1**2).mean()):6.3f} m " + f"95%|d|={np.percentile(np.abs(var_l1),95):6.3f} m n={len(var_l1)}") + if len(dd_p1): + print(f" P1 code DD (incl. master code noise): " + f"median|d|={np.median(np.abs(dd_p1)):6.3f} m " + f"RMS={np.sqrt((dd_p1**2).mean()):6.3f} m " + f"95%|d|={np.percentile(np.abs(dd_p1),95):6.3f} m") + if not args.keep_synth: + synth_path.unlink(missing_ok=True) + + +if __name__ == "__main__": + main() diff --git a/experiments/vbs_iono.py b/experiments/vbs_iono.py new file mode 100644 index 0000000..c9c6319 --- /dev/null +++ b/experiments/vbs_iono.py @@ -0,0 +1,295 @@ +#!/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() diff --git a/experiments/vbs_synth.py b/experiments/vbs_synth.py new file mode 100644 index 0000000..a99c8de --- /dev/null +++ b/experiments/vbs_synth.py @@ -0,0 +1,343 @@ +#!/usr/bin/env python3 +""" +vbs_synth.py - Virtual base station synthesis: geometric displacement engine. + +Patches a RINEX 2.11 observation file so its observations appear to have been +collected at a different (virtual) position: per epoch/satellite, the change in +geometric range (independent light-time iteration + Earth-rotation correction at +each position, SP3 orbits) is added to every code observable (metres) and every +carrier-phase observable (cycles, per-signal wavelength). Everything else - +LLI/SSI flags, S/D observables, obs-type order, epoch structure - is preserved. +Satellites without SP3 orbits or (GLONASS) without a known frequency slot are +dropped, with epoch satellite lists rewritten accordingly. + +The master receiver's clock is inherited by the virtual station (harmless: it is +one clock, absorbed by the processor's per-epoch receiver-clock estimation; see +the --perturb-clock-ms test mode). The atmospheric content of the master's +observations is NOT changed by this engine - a geometric-only virtual base +carries the master's atmosphere verbatim (Phase 2 applies interpolated +corrections on top via --iono). + +Usage: + vbs_synth.py master.obs --xyz X Y Z --sp3 FILE --nav BRDC.rnx -o out.obs + [--master-xyz X Y Z] [--perturb-clock-ms N] [--marker NAME] +""" + +import argparse +import math +import sys +from datetime import datetime, timedelta, timezone +from pathlib import Path + +import numpy as np + +C = 299792458.0 +OMEGA_E = 7.2921151467e-5 # rad/s +GAMMA_12 = (1575.42 / 1227.60) ** 2 + +# carrier frequency by (system, band); GLONASS handled via slot +FREQ = { + ("G", "1"): 1575.42e6, ("G", "2"): 1227.60e6, ("G", "5"): 1176.45e6, + ("E", "1"): 1575.42e6, ("E", "5"): 1176.45e6, # E1, E5a + ("S", "1"): 1575.42e6, +} + + +def glo_freq(band, slot): + if band == "1": + return 1602.0e6 + slot * 562.5e3 + if band == "2": + return 1246.0e6 + slot * 437.5e3 + return None + + +# ---------------------------------------------------------------------------- +# SP3 orbits +# ---------------------------------------------------------------------------- + +class Sp3: + """SP3 position source with sliding-window Lagrange interpolation.""" + + def __init__(self, path): + self.t0 = None + times = [] + pos = {} # sat -> {tidx: (x,y,z) m} + tidx = -1 + for line in open(path, errors="replace"): + if line.startswith("*"): + p = line.split() + t = datetime(int(p[1]), int(p[2]), int(p[3]), int(p[4]), int(p[5]), + int(float(p[6])), tzinfo=timezone.utc) + if self.t0 is None: + self.t0 = t + times.append((t - self.t0).total_seconds()) + tidx += 1 + elif line.startswith("P") and tidx >= 0: + sat = line[1:4].replace(" ", "0") + try: + x, y, z = (float(line[4:18]), float(line[18:32]), float(line[32:46])) + except ValueError: + continue + if abs(x) > 900000 or (x == 0 and y == 0): # bad/absent + continue + pos.setdefault(sat, {})[tidx] = (x * 1e3, y * 1e3, z * 1e3) + self.times = np.array(times) + self.sats = {} + n = len(times) + for sat, d in pos.items(): + if len(d) < n * 0.9: # incomplete arcs: drop satellite + continue + arr = np.full((n, 3), np.nan) + for i, xyz in d.items(): + arr[i] = xyz + if np.isnan(arr).any(): + continue + self.sats[sat] = arr + + def pos(self, sat, t_sec, order=10): + """Satellite ECEF position (m) at t_sec (seconds from SP3 start).""" + arr = self.sats.get(sat) + if arr is None: + return None + i = np.searchsorted(self.times, t_sec) + half = order // 2 + lo = max(0, min(i - half, len(self.times) - order)) + idx = slice(lo, lo + order) + tt = self.times[idx] + out = np.empty(3) + # Lagrange interpolation per coordinate + for k in range(3): + y = arr[idx, k] + acc = 0.0 + for j in range(len(tt)): + lj = 1.0 + for m in range(len(tt)): + if m != j: + lj *= (t_sec - tt[m]) / (tt[j] - tt[m]) + acc += y[j] * lj + out[k] = acc + return out + + +def glonass_slots(nav_path): + """sat 'Rnn' -> frequency slot number, from a RINEX 3 mixed nav file.""" + slots = {} + lines = open(nav_path, errors="replace").read().splitlines() + i = 0 + while i < len(lines): + line = lines[i] + if line[:1] == "R" and len(line) > 23 and line[1:3].strip().isdigit(): + sat = "R" + line[1:3].replace(" ", "0") + if sat not in slots and i + 2 < len(lines): + orbit2 = lines[i + 2] + try: # 4th 19-char field on broadcast orbit 2 + slots[sat] = int(float(orbit2[4 + 3 * 19: 4 + 4 * 19].replace("D", "E"))) + except ValueError: + pass + i += 4 + else: + i += 1 + return slots + + +# ---------------------------------------------------------------------------- +# Geometry +# ---------------------------------------------------------------------------- + +def geometric_range(sp3, sat, t_rx_sec, x_rcv): + """Light-time-iterated, Earth-rotation-corrected range (m); None if no orbit.""" + tau = 0.075 + rho = None + for _ in range(3): + xs = sp3.pos(sat, t_rx_sec - tau) + if xs is None: + return None + ang = OMEGA_E * tau # rotate satellite into reception-time frame + ca, sa = math.cos(ang), math.sin(ang) + xr = np.array([xs[0] * ca + xs[1] * sa, -xs[0] * sa + xs[1] * ca, xs[2]]) + rho = float(np.linalg.norm(xr - x_rcv)) + tau = rho / C + return rho + + +# ---------------------------------------------------------------------------- +# RINEX 2.11 observation patcher +# ---------------------------------------------------------------------------- + +def parse_header(lines): + """Returns (obs_types, approx_xyz, end_idx).""" + types, xyz = [], None + for i, line in enumerate(lines): + label = line[60:].strip() + if label == "# / TYPES OF OBSERV": + p = line[:60].split() + if types == [] and p and p[0].isdigit(): + types.extend(p[1:]) + else: + types.extend(p) + elif label == "APPROX POSITION XYZ": + xyz = tuple(float(v) for v in line[:60].split()) + elif label == "END OF HEADER": + return types, xyz, i + raise ValueError("no END OF HEADER") + + +def fmt_obs(value, lli, ssi): + if value is None: + return " " * 14 + lli + ssi + return f"{value:14.3f}" + lli + ssi + + +def synthesize(master_path, out_path, virtual_xyz, sp3, slots, master_xyz=None, + perturb_ms=0.0, marker="VBS"): + lines = open(master_path, errors="replace").read().splitlines() + types, hdr_xyz, hdr_end = parse_header(lines) + ntypes = len(types) + nlines_per_sat = (ntypes + 4) // 5 + x_mas = np.array(master_xyz if master_xyz else hdr_xyz) + x_vrs = np.array(virtual_xyz) + disp_km = np.linalg.norm(x_vrs - x_mas) / 1e3 + + out = [] + for line in lines[:hdr_end + 1]: + label = line[60:].strip() + if label == "APPROX POSITION XYZ": + out.append(f"{x_vrs[0]:14.4f}{x_vrs[1]:14.4f}{x_vrs[2]:14.4f}" + + " " * 18 + "APPROX POSITION XYZ") + elif label == "MARKER NAME": + out.append(f"{marker:<60}MARKER NAME") + elif label == "END OF HEADER": + out.append(f"{'VBS synthesized from ' + Path(master_path).name:<60}COMMENT") + out.append(f"{f'displacement {disp_km:.3f} km; geometric only':<60}COMMENT") + out.append(line) + else: + out.append(line) + + i = hdr_end + 1 + stats = {"epochs": 0, "sats": 0, "dropped": {}, "dr_min": 1e9, "dr_max": -1e9} + 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: # event records: copy verbatim + out.append(line) + for k in range(nsat): + i += 1 + out.append(lines[i]) + i += 1 + continue + # epoch time (GPS time system) + p = line[:26].split() + t_rx = datetime(2000 + int(p[0]), int(p[1]), int(p[2]), int(p[3]), int(p[4]), + tzinfo=timezone.utc) + timedelta(seconds=float(p[5])) + t_sec = (t_rx - sp3.t0).total_seconds() + # satellite list (12 per line, continuations) + 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 + # per-satellite observation blocks + keep = [] + for sat in sats: + block = lines[i:i + nlines_per_sat] + i += nlines_per_sat + sys_id = sat[0] + if sys_id == "R" and sat not in slots: + stats["dropped"][sat] = stats["dropped"].get(sat, 0) + 1 + continue + rho_m = geometric_range(sp3, sat, t_sec, x_mas) + if rho_m is None: + stats["dropped"][sat] = stats["dropped"].get(sat, 0) + 1 + continue + rho_v = geometric_range(sp3, sat, t_sec, x_vrs) + dr = rho_v - rho_m + stats["dr_min"] = min(stats["dr_min"], dr) + stats["dr_max"] = max(stats["dr_max"], dr) + # parse + patch the observation fields + fields = [] + ok = True + for ti, typ in enumerate(types): + ln = block[ti // 5] if ti // 5 < len(block) else "" + seg = ln[(ti % 5) * 16:(ti % 5) * 16 + 16].ljust(16) + raw, lli, ssi = seg[:14], seg[14], seg[15] + val = float(raw) if raw.strip() else None + if val is not None: + if typ[0] in ("C", "P"): + val += dr + C * perturb_ms * 1e-3 + elif typ[0] == "L": + f = (glo_freq(typ[1], slots[sat]) if sys_id == "R" + else FREQ.get((sys_id, typ[1]))) + if f is None: + ok = False + break + lam = C / f + val += dr / lam + C * perturb_ms * 1e-3 / lam + fields.append((val, lli, ssi)) + if not ok: + stats["dropped"][sat] = stats["dropped"].get(sat, 0) + 1 + continue + keep.append((sat, fields)) + # rewrite epoch record with surviving satellites + if not keep: + continue + stats["epochs"] += 1 + stats["sats"] += len(keep) + ids = [s for s, _ in keep] + head = line[:29] + f"{len(ids):3d}" + for k in range(0, len(ids), 12): + seg = "".join(ids[k:k + 12]) + out.append((head if k == 0 else " " * 32) + seg) + for sat, fields in keep: + for li in range(nlines_per_sat): + chunk = fields[li * 5:(li + 1) * 5] + out.append("".join(fmt_obs(*f) for f in chunk).rstrip()) + + Path(out_path).write_text("\n".join(out) + "\n") + return stats, disp_km + + +def main(): + ap = argparse.ArgumentParser(description=__doc__, + formatter_class=argparse.RawDescriptionHelpFormatter) + ap.add_argument("master", type=Path, help="master RINEX 2.11 obs file") + ap.add_argument("--xyz", nargs=3, type=float, required=True, metavar=("X", "Y", "Z"), + help="virtual position, ECEF metres (same frame as --master-xyz)") + ap.add_argument("--sp3", type=Path, required=True) + ap.add_argument("--nav", type=Path, required=True, + help="RINEX 3 mixed nav (for GLONASS frequency slots)") + ap.add_argument("-o", "--out", type=Path, required=True) + ap.add_argument("--master-xyz", nargs=3, type=float, metavar=("X", "Y", "Z"), + help="master position override (else RINEX header position)") + ap.add_argument("--perturb-clock-ms", type=float, default=0.0, + help="TEST MODE: add a simulated receiver-clock offset to all " + "code+phase observables") + ap.add_argument("--marker", default="VBS0") + args = ap.parse_args() + + sp3 = Sp3(args.sp3) + slots = glonass_slots(args.nav) + stats, disp = synthesize(args.master, args.out, args.xyz, sp3, slots, + master_xyz=args.master_xyz, + perturb_ms=args.perturb_clock_ms, marker=args.marker) + dropped = sum(stats["dropped"].values()) + print(f"[vbs_synth] {args.master.name} -> {args.out.name}: " + f"displacement {disp:.3f} km, {stats['epochs']} epochs, " + f"{stats['sats']} sat-epochs kept, {dropped} dropped " + f"({len(stats['dropped'])} sats), dRho [{stats['dr_min']:.1f}, " + f"{stats['dr_max']:.1f}] m") + + +if __name__ == "__main__": + main() diff --git a/gcgc.env.example b/gcgc.env.example new file mode 100644 index 0000000..e05d12b --- /dev/null +++ b/gcgc.env.example @@ -0,0 +1,8 @@ +# GCGC RTN (rtn.usm.edu) NTRIP credentials. +# +# Setup: cp gcgc.env.example gcgc.env then fill in real values. +# gcgc.env is gitignored - never commit real credentials. +# stream_gps.py loads this file automatically; already-set environment +# variables take precedence. +GCGC_USER=your_username_here +GCGC_PASS=your_password_here diff --git a/process_gps.py b/process_gps.py index dcf854b..034a9c9 100644 --- a/process_gps.py +++ b/process_gps.py @@ -151,6 +151,29 @@ def haversine_km(lat1, lon1, lat2, lon2): return 12742 * math.asin(math.sqrt(a)) +def llh_to_ecef(lat, lon, h): + """Geodetic (deg, m) -> ECEF (m), WGS84.""" + a, f = 6378137.0, 1 / 298.257223563 + e2 = f * (2 - f) + phi, lam = math.radians(lat), math.radians(lon) + n = a / math.sqrt(1 - e2 * math.sin(phi) ** 2) + return ((n + h) * math.cos(phi) * math.cos(lam), + (n + h) * math.cos(phi) * math.sin(lam), + (n * (1 - e2) + h) * math.sin(phi)) + + +def ecef_to_llh(x, y, z): + """ECEF -> geodetic lat/lon (deg), height (m); Bowring's method, WGS84.""" + a, f = 6378137.0, 1 / 298.257223563 + b, e2 = a * (1 - f), f * (2 - f) + ep2 = (a * a - b * b) / (b * b) + p = math.hypot(x, y) + th = math.atan2(a * z, b * p) + lat = math.atan2(z + ep2 * b * math.sin(th) ** 3, p - e2 * a * math.cos(th) ** 3) + n = a / math.sqrt(1 - e2 * math.sin(lat) ** 2) + return math.degrees(lat), math.degrees(math.atan2(y, x)), p / math.cos(lat) - n + + def llh_to_enu(lat, lon, h, lat0, lon0, h0): """Small-extent geodetic -> local ENU (m) about (lat0, lon0, h0). Arrays OK.""" a, f = 6378137.0, 1 / 298.257223563 @@ -175,6 +198,7 @@ def parse_ubx(ubx_path): rawx_hdr, rawx_meas = [], [] pvt, sat, sig_rows, rf = [], [], [], [] sfrbx_gnss = Counter() + cur_utc = None # running epoch time for sat/sig/rf tagging with open(ubx_path, "rb") as f: ubr = UBXReader(f, protfilter=UBX_PROTOCOL, quitonerror=0) @@ -200,8 +224,9 @@ def parse_ubx(ubx_path): elif ident == "RXM-SFRBX": sfrbx_gnss[m.gnssId] += 1 elif ident == "NAV-PVT": - pvt.append((datetime(m.year, m.month, m.day, m.hour, m.min, m.second, - tzinfo=timezone.utc), + cur_utc = datetime(m.year, m.month, m.day, m.hour, m.min, m.second, + tzinfo=timezone.utc) + pvt.append((cur_utc, m.fixType, m.carrSoln, m.diffSoln, m.numSV, m.lat, m.lon, m.height / 1000.0, m.hMSL / 1000.0, m.hAcc / 1000.0, m.vAcc / 1000.0, m.tAcc, m.pDOP, @@ -209,7 +234,7 @@ def parse_ubx(ubx_path): elif ident == "NAV-SAT": for i in range(1, m.numSvs + 1): s = f"_{i:02d}" - sat.append((getattr(m, "gnssId" + s), getattr(m, "svId" + s), + sat.append((cur_utc, getattr(m, "gnssId" + s), getattr(m, "svId" + s), getattr(m, "cno" + s), getattr(m, "elev" + s), getattr(m, "azim" + s), getattr(m, "prRes" + s), getattr(m, "qualityInd" + s), getattr(m, "svUsed" + s), @@ -217,7 +242,7 @@ def parse_ubx(ubx_path): elif ident == "NAV-SIG": for i in range(1, m.numSigs + 1): s = f"_{i:02d}" - sig_rows.append((getattr(m, "gnssId" + s), getattr(m, "svId" + s), + sig_rows.append((cur_utc, getattr(m, "gnssId" + s), getattr(m, "svId" + s), getattr(m, "sigId" + s), getattr(m, "cno" + s), getattr(m, "qualityInd" + s), getattr(m, "corrSource" + s), getattr(m, "ionoModel" + s), getattr(m, "prUsed" + s), @@ -225,7 +250,7 @@ def parse_ubx(ubx_path): elif ident == "MON-RF": for i in range(1, m.nBlocks + 1): s = f"_{i:02d}" - rf.append((getattr(m, "blockId" + s), getattr(m, "jammingState" + s), + rf.append((cur_utc, getattr(m, "blockId" + s), getattr(m, "jammingState" + s), getattr(m, "antStatus" + s), getattr(m, "antPower" + s), getattr(m, "noisePerMS" + s), getattr(m, "agcCnt" + s), getattr(m, "jamInd" + s))) @@ -240,12 +265,13 @@ def parse_ubx(ubx_path): "pvt": pd.DataFrame(pvt, columns=["utc", "fixType", "carrSoln", "diffSoln", "numSV", "lat", "lon", "height", "hMSL", "hAcc", "vAcc", "tAcc_ns", "pDOP", "gSpeed"]), - "sat": pd.DataFrame(sat, columns=["gnss", "sv", "cno", "elev", "azim", "prRes", - "quality", "used", "health"]), - "sig": pd.DataFrame(sig_rows, columns=["gnss", "sv", "sig", "cno", "quality", - "corrSource", "ionoModel", "prUsed", "prRes"]), - "rf": pd.DataFrame(rf, columns=["blockId", "jammingState", "antStatus", "antPower", - "noisePerMS", "agcCnt", "jamInd"]), + "sat": pd.DataFrame(sat, columns=["utc", "gnss", "sv", "cno", "elev", "azim", + "prRes", "quality", "used", "health"]), + "sig": pd.DataFrame(sig_rows, columns=["utc", "gnss", "sv", "sig", "cno", + "quality", "corrSource", "ionoModel", + "prUsed", "prRes"]), + "rf": pd.DataFrame(rf, columns=["utc", "blockId", "jammingState", "antStatus", + "antPower", "noisePerMS", "agcCnt", "jamInd"]), } # session window in UTC from RAWX (receiver time is GPS-aligned) @@ -364,7 +390,7 @@ def rover_antenna(): return None -def write_conf(path, ant1=None, ant2=None, atx=None): +def write_conf(path, ant1=None, ant2=None, atx=None, base_from_header=False): lines = [ "pos1-elmask =15", "pos1-navsys =45", # GPS+GLO+GAL+BDS @@ -372,6 +398,8 @@ def write_conf(path, ant1=None, ant2=None, atx=None): "pos2-gloarmode =off", # mixed-manufacturer base/rover: GLONASS AR off "out-outstat =residual", ] + if base_from_header: # no -r on CLI: take base pos from RINEX header + lines.append("ant2-postype =rinexhead") if ant1: lines.append(f"ant1-anttype ={ant1}") if ant2: @@ -514,6 +542,107 @@ def download_base_hours(station, t0, t1, dest_dir, warnings): return files +def rinex_first_obs(path): + """UTC datetime of TIME OF FIRST OBS from a RINEX obs header, or None.""" + try: + with open(path, errors="replace") as f: + for line in f: + if "TIME OF FIRST OBS" in line: + p = line.split() + return datetime(int(p[0]), int(p[1]), int(p[2]), int(p[3]), + int(p[4]), int(float(p[5])), tzinfo=timezone.utc) + if "END OF HEADER" in line: + break + except (OSError, ValueError, IndexError): + pass + return None + + +def find_gcgc_base(capture, t_start, t_end): + """Auto-discover local GCGC/VRS base data for the parallel PPK method: + the capture's own recorded correction stream (stream_gps.py writes + .rtcm3 alongside .ubx), plus files dropped into gcgc_base/ + whose observation window overlaps the capture.""" + found = [] + sidecar = capture.with_suffix(".rtcm3") + if sidecar.is_file() and sidecar.stat().st_size: + found.append(sidecar) + gdir = SCRIPT_DIR / "gcgc_base" + if gdir.is_dir(): + for f in sorted(gdir.iterdir()): + suf = f.suffix.lower() + if suf in (".rtcm3", ".rtcm"): + found.append(f) + elif suf in (".obs", ".rnx", ".crx") or re.fullmatch(r"\.\d\do", suf): + t = rinex_first_obs(f) + if t is not None and t <= t_end and t >= t_start - timedelta(hours=26): + found.append(f) + return found + + +def prepare_local_base(outdir, files, t0, warnings, tag): + """Stage local base observation file(s): RINEX, or RTCM3 correction logs + (e.g. a recorded GCGC VRS stream) which are converted via convbin.""" + base_dir = outdir / "base" / tag + base_dir.mkdir(parents=True, exist_ok=True) + staged = [] + for i, f in enumerate(Path(p) for p in files): + if not f.is_file(): + warnings.append(f"base file not found: {f}") + continue + if f.suffix.lower() in (".rtcm3", ".rtcm"): + conv = base_dir / f"base_{i:02d}.obs" + if not conv.is_file(): + run([BIN / "convbin", "-r", "rtcm3", "-v", "3.04", + "-tr", t0.strftime("%Y/%m/%d"), t0.strftime("%H:%M:%S"), + "-o", conv, f]) + if not conv.is_file(): + warnings.append(f"RTCM->RINEX conversion failed for {f.name}") + continue + staged.append(conv) + else: + local = base_dir / f"base_{i:02d}{f.suffix}" + if not local.is_file(): + shutil.copy(f, local) # space-free names for combine_rinex.sh + staged.append(local) + if not staged: + return None + # never combine files from DIFFERENT base positions (e.g. several VRS points): + # keep only files agreeing with the first-seen position within 100 m + positions = [] + for f in staged: + for line in open(f, errors="replace"): + if "APPROX POSITION XYZ" in line: + try: + positions.append(np.array([float(v) for v in line[:60].split()])) + except ValueError: + positions.append(None) + break + if "END OF HEADER" in line: + positions.append(None) + break + if len([p for p in positions if p is not None]) > 1: + ref = next(p for p in positions if p is not None) + same = [f for f, p in zip(staged, positions) + if p is None or np.linalg.norm(p - ref) < 100.0] + if len(same) < len(staged): + warnings.append(f"{tag}: {len(staged) - len(same)} base file(s) at a " + "different position excluded from combining (multiple " + "VRS points? process them separately via --base-obs)") + staged = same + if len(staged) == 1: + return staged[0] + combined = base_dir / "base_combined.obs" + if not combined.is_file(): + r = run(["bash", COMBINE_SH, combined.name, *[f.name for f in staged]], + cwd=base_dir) + if not combined.is_file() or r.returncode != 0: + warnings.append(f"combine_rinex.sh failed for local base files: " + f"{r.stdout[-400:]}") + return None + return combined + + def stage_ppk(outdir, obs, nav, data, args, warnings): result = {"status": "skipped"} t0 = data["t_start"] - timedelta(minutes=5) @@ -563,6 +692,141 @@ def stage_ppk(outdir, obs, nav, data, args, warnings): "Re-run this script later; completed stages are cached.") log("ppk: " + result["message"]) return result + return finish_ppk(outdir, obs, nav, data, args, warnings, result, + base_obs, base_station, mean_lat, mean_lon) + + +def fetch_esa_sp3(outdir, day, warnings): + """Multi-constellation ESA orbit product (FIN preferred, else RAP, G+R+).""" + week = int((day - GPS_EPOCH).days // 7) + yyyyddd = day.strftime("%Y%j") + prod_dir = outdir / "igs" + prod_dir.mkdir(exist_ok=True) + for kind in ("FIN", "RAP"): + for f in prod_dir.glob(f"ESA0OPS{kind}_{yyyyddd}0000*ORB.SP3"): + return f + listing = fetch(f"http://navigation-office.esa.int/products/gnss-products/{week}/", + timeout=60) + if listing is None: + return None + names = set(re.findall(r'href="(ESA0OPS(?:FIN|RAP)_[^"]+ORB\.SP3\.gz)"', + listing.decode(errors="replace"))) + for kind in ("FIN", "RAP"): + cand = [n for n in names if kind in n and yyyyddd + "0000" in n] + if cand: + gz = prod_dir / cand[0] + if fetch(f"http://navigation-office.esa.int/products/gnss-products/" + f"{week}/{cand[0]}", gz): + out = gz.with_suffix("") + out.write_bytes(gzip.decompress(gz.read_bytes())) + return out + return None + + +def stage_ppk_vrs(outdir, obs, nav, data, args, warnings, ppk_result): + """Standard GCGC-provided VRS solution: synthesize a virtual base at the + capture's own position from the GCGC master station the NGS stage already + fetched (engine validated in experiments/: gates G1/G2a green; benchmarked + equivalent to Trimble Pivot VRS at station MARY, 2026-07-27). + Only provided when GCGC is available: the master must belong to the GCGC + network (agency header or a Mississippi station).""" + result = {"status": "skipped"} + st = ppk_result.get("station") + base_obs = ppk_result.get("base_obs") + master_xyz = ppk_result.get("base_xyz") + if not st or not base_obs or not Path(base_obs).is_file(): + return result + # --- GCGC availability gate + agency = "" + for line in open(base_obs, errors="replace"): + if "OBSERVER / AGENCY" in line: + agency = line[:60] + if "END OF HEADER" in line: + break + if "GCGC" not in agency.upper() and st.get("state") != "MS": + log(f"ppk_vrs: master {st['id']} is not a GCGC station - " + "GCGC VRS solution not provided") + return result + if master_xyz is None: + warnings.append("ppk_vrs: master ITRF position unavailable - skipped " + "(frame-consistent synthesis requires it)") + return result + + sp3_path = fetch_esa_sp3(outdir, data["t_start"], warnings) + if sp3_path is None: + result["status"] = "pending" + result["message"] = ("ESA orbit products for this date not yet published " + "(~1 day lag) - GCGC VRS solution added on a later re-run") + log("ppk_vrs: " + result["message"]) + return result + + sys.path.insert(0, str(SCRIPT_DIR / "experiments")) + from vbs_synth import Sp3, glonass_slots, synthesize + + pvt = data["pvt"] + vxyz = llh_to_ecef(float(pvt["lat"].mean()), float(pvt["lon"].mean()), + float(pvt["height"].mean())) + vrs_obs = outdir / "vrs_base.obs" + if not vrs_obs.is_file(): + sp3 = Sp3(sp3_path) + slots = glonass_slots(nav) # rover.nav is RINEX 3 mixed nav + stats, disp = synthesize(base_obs, vrs_obs, vxyz, sp3, slots, + master_xyz=master_xyz, marker="VRS0") + dropped = sum(stats["dropped"].values()) + log(f"ppk_vrs: synthesized virtual base at capture position from " + f"{st['id']} ({disp:.1f} km displacement), {stats['epochs']} epochs, " + f"{dropped} sat-epochs without orbits dropped") + if stats["epochs"] == 0: + warnings.append("ppk_vrs: synthesis produced no epochs - skipped") + vrs_obs.unlink(missing_ok=True) + result["status"] = "failed" + return result + + base_station = {"id": f"VRS@{st['id']}", "km": 0.0, "state": st.get("state", "-"), + "status": "synthesized", "xyz": vxyz} + return finish_ppk(outdir, obs, nav, data, args, warnings, result, + vrs_obs, base_station, float(pvt["lat"].mean()), + float(pvt["lon"].mean()), prefix="ppk_vrs") + + +def stage_ppk_gcgc(outdir, obs, nav, data, args, warnings): + """Parallel PPK method using local GCGC data: auto-discovered (recorded VRS + sidecar, gcgc_base/ folder) or explicitly supplied via --base-obs.""" + result = {"status": "skipped"} + t0 = data["t_start"] - timedelta(minutes=5) + mean_lat, mean_lon = data["pvt"]["lat"].mean(), data["pvt"]["lon"].mean() + + files = ([Path(p) for p in args.base_obs] if args.base_obs + else find_gcgc_base(args.capture, data["t_start"], data["t_end"])) + if not files: + return result + log(f"ppk_gcgc: local base data found: {', '.join(f.name for f in files)}") + base_obs = prepare_local_base(outdir, files, t0, warnings, tag="gcgc") + if base_obs is None: + result["status"] = "failed" + return result + # distance from base RINEX header position, if present (VRS => ~0 km) + km = float("nan") + for line in open(base_obs, errors="replace"): + if "APPROX POSITION XYZ" in line: + try: + x, y, z = (float(v) for v in line[:60].split()) + lat0, lon0, _ = ecef_to_llh(x, y, z) + km = haversine_km(mean_lat, mean_lon, lat0, lon0) + except ValueError: + pass + if "END OF HEADER" in line: + break + base_station = {"id": "GCGC", "km": km, "state": "-", "status": "local", + "local": True} + return finish_ppk(outdir, obs, nav, data, args, warnings, result, + base_obs, base_station, mean_lat, mean_lon, prefix="ppk_gcgc") + + +def finish_ppk(outdir, obs, nav, data, args, warnings, result, + base_obs, base_station, mean_lat, mean_lon, prefix="ppk"): + t0 = data["t_start"] - timedelta(minutes=5) + t1 = data["t_end"] + timedelta(minutes=5) # base antenna type from its RINEX header; base ECEF from NGS coord file ant2 = None @@ -571,16 +835,26 @@ def stage_ppk(outdir, obs, nav, data, args, warnings): ant2 = line[20:40].rstrip() if "END OF HEADER" in line: break - epoch_year = data["t_start"].year + data["t_start"].timetuple().tm_yday / 365.25 - xyz = cors_station_xyz(base_station["id"], epoch_year) - if xyz is None: - warnings.append(f"NGS coord file for {base_station['id']} unavailable - " - "using base RINEX header position (NAD83, ~1.5 m frame offset)") + if base_station.get("xyz") is not None: # synthesized base: position known + xyz = base_station["xyz"] + elif base_station.get("local"): + xyz = tuple(args.base_xyz) if args.base_xyz else None + if xyz is None: + warnings.append(f"{prefix}: no --base-xyz given, using the base RINEX " + "header position (verify its datum; GCGC uses NAD83, " + "~1.5 m from ITRF in CONUS)") + else: + epoch_year = data["t_start"].year + data["t_start"].timetuple().tm_yday / 365.25 + xyz = cors_station_xyz(base_station["id"], epoch_year) + if xyz is None: + warnings.append(f"NGS coord file for {base_station['id']} unavailable - " + "using base RINEX header position (NAD83, ~1.5 m frame offset)") ant1 = rover_antenna() - merged = outdir / "merged.atx" + merged = outdir / f"{prefix}_merged.atx" build_merged_atx([ant1, ant2], merged) - conf = write_conf(outdir / "ppk.conf", ant1=ant1, ant2=ant2, atx=merged) + conf = write_conf(outdir / f"{prefix}.conf", ant1=ant1, ant2=ant2, atx=merged, + base_from_header=(xyz is None)) ts = ["-ts", t0.strftime("%Y/%m/%d"), t0.strftime("%H:%M:%S"), "-te", t1.strftime("%Y/%m/%d"), t1.strftime("%H:%M:%S")] @@ -596,24 +870,25 @@ def stage_ppk(outdir, obs, nav, data, args, warnings): warnings.append(f"--static flagged but onboard track spread is {spread95:.1f} m " f"(95th pct) - the static solution may be averaging real motion") elif not static and spread95 < 5.0: - log(f"ppk: track spread {spread95:.2f} m would pass as static, but processing " + log(f"{prefix}: track spread {spread95:.2f} m would pass as static, but processing " f"kinematic (default policy; use --static only for a verified fixed antenna)") - log(f"ppk: base={base_station['id']} ({base_station['km']:.1f} km) " + log(f"{prefix}: base={base_station['id']} ({base_station['km']:.1f} km) " f"ant2={ant2!r} static={static} (95% spread {spread95:.2f} m)") out = {} runs = [("kinematic", "2")] + ([("static", "3")] if static else []) for name, mode in runs: - pos = outdir / f"ppk_{name}.pos" + pos = outdir / f"{prefix}_{name}.pos" if not pos.is_file(): r = run([BIN / "rnx2rtkp", "-k", conf, "-p", mode, "-f", "3", "-t", "-u", *ts, *rpos, "-o", pos, obs, base_obs, nav]) if not pos.is_file() or not read_pos(pos).shape[0]: - warnings.append(f"rnx2rtkp ppk {name} produced no solution: {r.stderr[-800:]}") + warnings.append(f"rnx2rtkp {prefix} {name} produced no solution: " + f"{r.stderr[-800:]}") continue out[name] = read_pos(pos) qc = Counter(out[name]["Q"]) - log(f"ppk {name}: {len(out[name])} epochs, Q={dict(qc)}") + log(f"{prefix} {name}: {len(out[name])} epochs, Q={dict(qc)}") result.update({"status": "ok" if out else "failed", "station": base_station, "static": static, "solutions": out, "base_obs": str(base_obs), @@ -691,14 +966,30 @@ def stage_ppp(outdir, obs, nav, data, args, warnings): # Stage 6: track outputs # ---------------------------------------------------------------------------- -def stage_track(outdir, solutions, pvt, interval): - # best available source, in order of accuracy - for source in ("ppk_kinematic", "ppp", "single"): - if source in solutions and len(solutions[source]): - best = solutions[source] - break +def stage_track(outdir, solutions, pvt, interval, warmup_cut=None): + # Among PPK-family solutions, pick by measured quality (empirical horizontal + # scatter about own mean) - a nominally-better base can produce a worse + # solution (seen with a network-edge Pivot VRS), so rank on evidence. + # PPP/single remain ordered fallbacks. + def scatter_h(df): + dn, de, _ = llh_to_enu(df["lat"], df["lon"], df["height"], + float(df["lat"].mean()), float(df["lon"].mean()), 0) + return float(np.sqrt(np.mean((dn - dn.mean()) ** 2 + (de - de.mean()) ** 2))) + + ppk_family = [n for n in ("ppk_gcgc_kinematic", "ppk_vrs_kinematic", + "ppk_kinematic") + if n in solutions and len(solutions[n])] + if ppk_family: + source = min(ppk_family, + key=lambda n: scatter_h(after_warmup(solutions[n], warmup_cut)[0])) + best = solutions[source] else: - return None, None + for source in ("ppp", "single"): + if source in solutions and len(solutions[source]): + best = solutions[source] + break + else: + return None, None df = best.copy() df["utc_s"] = df["utc"].dt.round("1s") @@ -725,7 +1016,9 @@ def stage_track(outdir, solutions, pvt, interval): df.to_csv(f, index=False, columns=cols, float_format="%.9f", date_format="%Y-%m-%dT%H:%M:%S.%fZ") # KML for the chosen solution - kml_src = {"ppk_kinematic": "ppk_kinematic.pos", "ppp": "ppp.pos", + kml_src = {"ppk_gcgc_kinematic": "ppk_gcgc_kinematic.pos", + "ppk_vrs_kinematic": "ppk_vrs_kinematic.pos", + "ppk_kinematic": "ppk_kinematic.pos", "ppp": "ppp.pos", "single": "single.pos"}[source] run([BIN / "pos2kml", outdir / kml_src]) kml = (outdir / kml_src).with_suffix(".kml") @@ -746,6 +1039,79 @@ CONST_COLOR = dict(zip(CONST_ORDER, SERIES)) SEQ_RAMP = ["#cde2fb", "#9ec5f4", "#6da7ec", "#3987e5", "#256abf", "#184f95", "#0d366b"] SURFACE, INK, INK2, MUTED = "#fcfcfb", "#0b0b0b", "#52514e", "#898781" GRID, BASELINE = "#e1e0d9", "#c3c2b7" +# color follows the solution entity, never its rank in a given report +SOL_COLOR = {"ppk_kinematic": SERIES[0], "onboard": SERIES[1], + "ppk_gcgc_kinematic": SERIES[2], "single": SERIES[3], "ppp": SERIES[4], + "ppk_vrs_kinematic": "#008300"} # palette slot 6 (green) +# display/ranking order for solution tables (best geometry first) +SOLUTION_ORDER = ("ppk_gcgc_kinematic", "ppk_gcgc_static", + "ppk_vrs_kinematic", "ppk_vrs_static", + "ppk_kinematic", "ppk_static", "ppp", "single", "onboard") + +# Human-readable names for reports and figure legends. Abbreviations are +# expanded once in the "Solution methods" section of each report. +SOLUTION_DISPLAY = { + "ppk_gcgc_kinematic": "PPK - local GCGC base data (kinematic)", + "ppk_gcgc_static": "PPK - local GCGC base data (static)", + "ppk_vrs_kinematic": "PPK - synthesized virtual reference station (kinematic)", + "ppk_vrs_static": "PPK - synthesized virtual reference station (static)", + "ppk_kinematic": "PPK - CORS base station (kinematic)", + "ppk_static": "PPK - CORS base station (static)", + "ppp": "Precise point positioning", + "single": "Single-point (standalone)", + "onboard": "Receiver onboard solution", +} +SOL_SHORT = { + "ppk_gcgc_kinematic": "PPK - local GCGC base", + "ppk_vrs_kinematic": "PPK - synthesized virtual base", + "ppk_kinematic": "PPK - CORS base", + "ppp": "precise point positioning", + "single": "single-point", + "onboard": "receiver onboard", +} +IONO_DISPLAY = { + "GPSA": "GPS Klobuchar ionospheric model, alpha coefficients", + "GPSB": "GPS Klobuchar ionospheric model, beta coefficients", + "BDSA": "BeiDou Klobuchar ionospheric model, alpha coefficients", + "BDSB": "BeiDou Klobuchar ionospheric model, beta coefficients", + "GAL": "Galileo NeQuick ionospheric model coefficients", +} + + +def sol_disp(name): + return SOLUTION_DISPLAY.get(name, name) + + +def method_definitions(present): + """First-use expansion of every method abbreviation, for both reports.""" + defs = [] + if any(n.startswith("ppk") for n in present): + defs.append(("Post-processed kinematic (PPK)", + "carrier-phase differential positioning computed after the " + "fact against a base station")) + if any(n.startswith("ppk_vrs") for n in present): + defs.append(("Synthesized virtual reference station (VRS)", + "a base observation file synthesized at the capture position " + "from Gulf Coast Geospatial Center (GCGC) network station " + "data, giving a zero-length baseline")) + if any(n.startswith("ppk_gcgc") for n in present): + defs.append(("Local GCGC base data", + "Gulf Coast Geospatial Center (GCGC) data supplied locally: " + "a recorded correction stream or a Reference Data Shop " + "download")) + if "ppk_kinematic" in present or "ppk_static" in present: + defs.append(("CORS base station", + "the nearest National Geodetic Survey Continuously Operating " + "Reference Station (CORS)")) + if "ppp" in present: + defs.append(("Precise point positioning (PPP)", + "single-receiver technique using precise satellite orbit and " + "clock products instead of a base station")) + if "onboard" in present: + defs.append(("Receiver onboard solution", + "the u-blox receiver's own real-time output, aided by a " + "satellite-based augmentation system (SBAS)")) + return defs def _style_ax(ax, xlabel=None, ylabel=None, title=None): @@ -765,7 +1131,7 @@ def _style_ax(ax, xlabel=None, ylabel=None, title=None): ax.set_title(title, color=INK, fontsize=10, loc="left", fontweight="bold") -def stage_figures(outdir, data, solutions, track_source, diag): +def stage_figures(outdir, data, solutions, track_source, diag, warmup_cut=None): import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt @@ -776,6 +1142,23 @@ def stage_figures(outdir, data, solutions, track_source, diag): figdir.mkdir(exist_ok=True) figs = [] pvt, rawx, sat = data["pvt"], data["rawx"], data["sat"] + rf_df, hdr = data["rf"], data["rawx_hdr"] + n_warm = 0 + # warm-up window excluded from every plot (user requirement) + if warmup_cut is not None and len(hdr): + leap = int(hdr["leapS"].iloc[0]) + n_warm = sum(1 for w, t in zip(hdr["week"], hdr["rcvTow"]) + if gps_to_utc(int(w), float(t), leap) < warmup_cut) + if n_warm < len(hdr): + hdr = hdr.iloc[n_warm:] + rawx = rawx[rawx["epoch"] >= n_warm] + pvt, _ = after_warmup(pvt, warmup_cut) + sat = sat[sat["utc"] >= warmup_cut] + rf_df = rf_df[rf_df["utc"] >= warmup_cut] + solutions = {n: after_warmup(df, warmup_cut)[0] + for n, df in solutions.items()} + else: + n_warm = 0 ref = diag["position_error"]["reference"] refpos = (ref["lat"], ref["lon"], ref["ellip_height_m"]) best = solutions.get(track_source) if track_source else None @@ -794,11 +1177,16 @@ def stage_figures(outdir, data, solutions, track_source, diag): f = newfig(6.4, 6.0) ax = f.add_subplot(111) series = [] - if best is not None and len(best): - dn, de, _ = llh_to_enu(best["lat"], best["lon"], best["height"], *refpos) - series.append((track_source, de, dn, SERIES[0])) + shown = [track_source] if best is not None and len(best) else [] + shown += [n for n in ("ppk_gcgc_kinematic", "ppk_vrs_kinematic", "ppk_kinematic") + if n != track_source and n in solutions and len(solutions[n])] + for name in shown[:2]: # best + at most one other PPK method + df = solutions[name] + dn, de, _ = llh_to_enu(df["lat"], df["lon"], df["height"], *refpos) + series.append((SOL_SHORT.get(name, name), de, dn, + SOL_COLOR.get(name, SERIES[0]))) dn_o, de_o, _ = llh_to_enu(pvt["lat"], pvt["lon"], pvt["height"], *refpos) - series.append(("onboard (NAV-PVT)", de_o, dn_o, SERIES[1])) + series.append(("receiver onboard", de_o, dn_o, SOL_COLOR["onboard"])) for name, de, dn, color in series: ax.scatter(de, dn, s=6, c=color, alpha=0.5, linewidths=0, label=name) for name, de, dn, color in series: # mean markers with surface ring @@ -814,7 +1202,7 @@ def stage_figures(outdir, data, solutions, track_source, diag): circ = plt.Circle((de.mean(), dn.mean()), rad, fill=False, color=MUTED, linewidth=1.0, linestyle=ls) ax.add_patch(circ) - ax.annotate(f"{radius_label}{pct} {rad:.2f} m", (de.mean(), dn.mean() + rad), + ax.annotate(f"{pct}% radius: {rad:.2f} m", (de.mean(), dn.mean() + rad), color=INK2, fontsize=8, ha="center", va="bottom") ax.set_aspect("equal") ax.legend(loc="upper right", fontsize=8, frameon=False, labelcolor=INK2) @@ -833,10 +1221,10 @@ def stage_figures(outdir, data, solutions, track_source, diag): if best is not None and len(best) else None) for i, ax in enumerate(axes): if best_enu is not None: - ax.plot(best["utc"], best_enu[i], color=SERIES[0], linewidth=1.4, - label=track_source) + ax.plot(best["utc"], best_enu[i], color=SOL_COLOR.get(track_source, SERIES[0]), linewidth=1.4, + label=SOL_SHORT.get(track_source, track_source)) ax.plot(pvt["utc"], on_enu[i], color=SERIES[1], linewidth=1.4, - label="onboard") + label="receiver onboard") ax.axhline(0, color=BASELINE, linewidth=0.8) _style_ax(ax, None, f"{comps[i]} (m)") axes[0].legend(loc="lower right", bbox_to_anchor=(1, 1.02), fontsize=8, @@ -851,14 +1239,13 @@ def stage_figures(outdir, data, solutions, track_source, diag): f = newfig(8.5, 6.4) gs = f.add_gridspec(2, 1, height_ratios=[1.15, 1], hspace=0.42) ax = f.add_subplot(gs[0]) - hdr = data["rawx_hdr"] times = [gps_to_utc(int(w), t, int(hdr["leapS"].iloc[0])) for w, t in zip(hdr["week"], hdr["rcvTow"])] end_vals = {} for name in CONST_ORDER: gid = {v: k for k, v in GNSS_NAME.items()}[name] per_epoch = rawx[rawx["gnss"] == gid].groupby("epoch")["sv"].nunique() - counts = per_epoch.reindex(range(len(hdr)), fill_value=0) + counts = per_epoch.reindex(range(n_warm, n_warm + len(hdr)), fill_value=0) ax.plot(times, counts, color=CONST_COLOR[name], linewidth=1.4, label=name) end_vals[name] = float(counts.iloc[-1]) # selective direct labels: only when the line end is clearly separated @@ -892,12 +1279,14 @@ def stage_figures(outdir, data, solutions, track_source, diag): if not np.isnan(mat[i, j]): axc.text(j, i, f"{mat[i, j]:.1f}", ha="center", va="center", color=INK if mat[i, j] < 38 else SURFACE, fontsize=9) - axc.set_title("Mean C/N0 by constellation and band (dB-Hz)", color=INK, + axc.set_title("Mean carrier-to-noise density ratio (C/N0) by constellation " + "and band (dB-Hz)", color=INK, fontsize=10, loc="left", fontweight="bold") axc.tick_params(colors=MUTED, labelsize=9) for sp in axc.spines.values(): sp.set_visible(False) - save(f, "fig3_satellites.png", "Satellite counts and C/N0 by band") + save(f, "fig3_satellites.png", + "Satellites tracked and carrier-to-noise density ratio by band") # ---- 4. skyplot --------------------------------------------------------- f = newfig(6.4, 6.4) @@ -916,35 +1305,51 @@ def stage_figures(outdir, data, solutions, track_source, diag): ax.set_thetagrids([0, 90, 180, 270], ["N", "E", "S", "W"], color=INK2) ax.grid(color=GRID, linewidth=0.7) cb = f.colorbar(sc, ax=ax, shrink=0.7, pad=0.08) - cb.set_label("mean C/N0 (dB-Hz)", color=INK2, fontsize=9) + cb.set_label("mean carrier-to-noise density ratio (dB-Hz)", color=INK2, fontsize=9) cb.ax.tick_params(colors=MUTED, labelsize=8) cb.outline.set_visible(False) - ax.set_title("Sky plot (session-mean position per satellite)", color=INK, + ax.set_title("Sky plot (session-mean elevation/azimuth per satellite)", color=INK, fontsize=10, fontweight="bold") - save(f, "fig4_skyplot.png", "Sky plot colored by C/N0") + save(f, "fig4_skyplot.png", + "Sky plot colored by carrier-to-noise density ratio") # ---- 5. accuracy metrics over time ------------------------------------- f = newfig(8.5, 5.4) ax1, ax2 = f.subplots(2, 1, sharex=True) ax1.plot(pvt["utc"], pvt["hAcc"], color=SERIES[1], linewidth=1.4, - label="onboard hAcc") + label="onboard horizontal accuracy estimate (calibrated)") ax1.plot(pvt["utc"], pvt["vAcc"], color=SERIES[1], linewidth=1.0, - linestyle="--", label="onboard vAcc") + linestyle="--", label="onboard vertical accuracy estimate") if best is not None and "sdn" in best and best["sdn"].abs().sum() > 0: sdh = np.hypot(best["sdn"], best["sde"]) - ax1.plot(best["utc"], sdh, color=SERIES[0], linewidth=1.4, - label=f"{track_source} sd horiz") - ax1.plot(best["utc"], best["sdu"], color=SERIES[0], linewidth=1.0, - linestyle="--", label=f"{track_source} sd up") + ax1.plot(best["utc"], sdh, color=SOL_COLOR.get(track_source, SERIES[0]), linewidth=1.4, + label=f"{SOL_SHORT.get(track_source, track_source)} " + f"formal 1-sigma (optimistic)") + ax1.plot(best["utc"], best["sdu"], color=SOL_COLOR.get(track_source, SERIES[0]), linewidth=1.0, + linestyle="--", + label=f"{SOL_SHORT.get(track_source, track_source)} " + f"formal 1-sigma, up") + # measured dispersion reference line: the gap between this and the + # formal curve IS the message (filter self-estimates are uncalibrated) + dnb, deb, _ = llh_to_enu(best["lat"], best["lon"], best["height"], *refpos) + disp_rms = float(np.sqrt(np.mean((dnb - dnb.mean()) ** 2 + + (deb - deb.mean()) ** 2))) + ax1.axhline(disp_rms, color=MUTED, linewidth=1.2, linestyle=":") + ax1.annotate(f"measured horizontal dispersion RMS {disp_rms:.2f} m " + f"(formal 1-sigma underestimates by about " + f"{disp_rms/max(sdh.mean(),1e-3):.0f}x)", + (pvt["utc"].iloc[0], disp_rms), xytext=(4, 4), + textcoords="offset points", color=INK2, fontsize=8) ax1.legend(loc="upper right", fontsize=8, frameon=False, ncols=2, labelcolor=INK2) - _style_ax(ax1, None, "1-σ accuracy (m)", "Reported accuracy estimates") + _style_ax(ax1, None, "1-sigma estimate (m)", + "Self-reported precision (uncalibrated) vs measured dispersion") ax2.plot(pvt["utc"], pvt["pDOP"], color=SERIES[2], linewidth=1.4) - ax2.annotate("pDOP", (pvt["utc"].iloc[-1], pvt["pDOP"].iloc[-1]), + ax2.annotate("PDOP", (pvt["utc"].iloc[-1], pvt["pDOP"].iloc[-1]), xytext=(4, 0), textcoords="offset points", color=SERIES[2], fontsize=8, va="center") - _style_ax(ax2, "UTC", "pDOP", None) + _style_ax(ax2, "UTC", "position dilution of precision", None) f.autofmt_xdate() - save(f, "fig5_accuracy.png", "Accuracy estimates and DOP") + save(f, "fig5_accuracy.png", "Self-reported precision vs measured dispersion; DOP") # ---- 6. ionosphere + RF ------------------------------------------------- f = newfig(8.5, 6.2) @@ -966,7 +1371,7 @@ def stage_figures(outdir, data, solutions, track_source, diag): "Dual-frequency geometry-free ionospheric delay (session median per satellite)") axr = f.add_subplot(gs[1, 0]) - rfagg = data["rf"].groupby("blockId")["jamInd"].agg(["mean", "max"]) + rfagg = rf_df.groupby("blockId")["jamInd"].agg(["mean", "max"]) x = np.arange(len(rfagg)) axr.bar(x, rfagg["mean"], color=SERIES[0], width=0.55, label="mean") axr.scatter(x, rfagg["max"], color=INK2, s=24, zorder=5, label="max") @@ -975,13 +1380,15 @@ def stage_figures(outdir, data, solutions, track_source, diag): axr.axhline(255 * 0.6, color=MUTED, linewidth=0.8, linestyle="--") axr.annotate("interference concern >153", (0, 158), color=MUTED, fontsize=7) axr.legend(loc="upper right", fontsize=8, frameon=False, labelcolor=INK2) - _style_ax(axr, None, "jamInd (0-255)", "RF interference indicator") + _style_ax(axr, None, "jamming indicator (0-255)", + "Radio-frequency interference indicator") axj = f.add_subplot(gs[1, 1]) jitter_ms = (np.diff(hdr["rcvTow"].to_numpy()) - 1.0) * 1000 axj.hist(jitter_ms, bins=40, color=SERIES[0]) _style_ax(axj, "epoch interval error (ms)", "epochs", "Measurement epoch jitter") - save(f, "fig6_iono_rf.png", "Ionospheric delay, RF health, timing jitter") + save(f, "fig6_iono_rf.png", + "Ionospheric delay, radio-frequency health, timing jitter") log(f"figures: {len(figs)} PNGs -> {figdir}") return figs @@ -991,6 +1398,18 @@ def stage_figures(outdir, data, solutions, track_source, diag): # Stage 7: diagnostics + report # ---------------------------------------------------------------------------- +def after_warmup(df, cut, time_col="utc"): + """Drop rows before the warm-up cutoff (quality metrics only - data outputs + keep everything). Falls back to the full data if the session is shorter + than the warm-up window.""" + if df is None or cut is None or not len(df): + return df, False + out = df[df[time_col] >= cut] + if not len(out): + return df, True + return out, False + + def solution_stats(df, ref): """Error metrics for one solution DataFrame about reference (lat, lon, h).""" if df is None or not len(df): @@ -1052,22 +1471,51 @@ def geometry_free_iono(rawx): def stage_diagnostics(outdir, data, rinex_info, solutions, ppk_result, track_source, - warnings, args): + warnings, args, gcgc_result=None, vrs_result=None, + warmup_cut=None): rawx, pvt, sat, sig, rf = data["rawx"], data["pvt"], data["sat"], data["sig"], data["rf"] + # Quality metrics exclude the receiver warm-up window (user requirement: + # first N minutes are still converging). Data outputs are NOT filtered. + pvt_m, warm_fb = after_warmup(pvt, warmup_cut) + if warm_fb: + warnings.append(f"session shorter than the {args.warmup_min:g} min warm-up " + "window - quality metrics use the full session") + msol = {} + for name, df in solutions.items(): + msol[name], _ = after_warmup(df, warmup_cut) + if warmup_cut is not None and not warm_fb: + hdr = data["rawx_hdr"] + leap = int(hdr["leapS"].iloc[0]) if len(hdr) else 18 + n_warm = sum(1 for w, t in zip(hdr["week"], hdr["rcvTow"]) + if gps_to_utc(int(w), float(t), leap) < warmup_cut) if len(hdr) else 0 + rawx = rawx[rawx["epoch"] >= n_warm] + sat = sat[sat["utc"] >= warmup_cut] + sig = sig[sig["utc"] >= warmup_cut] + rf = rf[rf["utc"] >= warmup_cut] + # --- reference position: PPK static end-state (only if --static was flagged) - # > PPK kinematic mean > onboard mean. In kinematic mode the reference is a - # dispersion anchor only - "scatter" includes any real platform motion. + # > PPK kinematic mean > onboard mean; GCGC (shortest baseline) preferred + # over NGS at each tier. In kinematic mode the reference is a dispersion + # anchor only - "scatter" includes any real platform motion. ref_source = "onboard session mean" - ref = (float(pvt["lat"].mean()), float(pvt["lon"].mean()), float(pvt["height"].mean())) - if "ppk_static" in solutions and len(solutions["ppk_static"]): - last = solutions["ppk_static"].iloc[-1] - ref = (float(last["lat"]), float(last["lon"]), float(last["height"])) - ref_source = f"PPK static final epoch (Q={QUALITY_NAME.get(int(last['Q']))})" - elif "ppk_kinematic" in solutions and len(solutions["ppk_kinematic"]): - d = solutions["ppk_kinematic"] - ref = (float(d["lat"].mean()), float(d["lon"].mean()), float(d["height"].mean())) - ref_source = "PPK kinematic session mean (dispersion anchor; includes platform motion)" + ref = (float(pvt_m["lat"].mean()), float(pvt_m["lon"].mean()), + float(pvt_m["height"].mean())) + for name in ("ppk_gcgc_static", "ppk_vrs_static", "ppk_static"): + if name in msol and len(msol[name]): + last = msol[name].iloc[-1] + ref = (float(last["lat"]), float(last["lon"]), float(last["height"])) + ref_source = (f"{name} final epoch (Q={QUALITY_NAME.get(int(last['Q']))})") + break + else: + for name in ("ppk_gcgc_kinematic", "ppk_vrs_kinematic", "ppk_kinematic"): + if name in msol and len(msol[name]): + d = msol[name] + ref = (float(d["lat"].mean()), float(d["lon"].mean()), + float(d["height"].mean())) + ref_source = (f"{name} session mean (dispersion anchor; " + "includes platform motion)") + break # --- constellation metrics per_gnss = {} @@ -1113,24 +1561,40 @@ def stage_diagnostics(outdir, data, rinex_info, solutions, ppk_result, track_sou # --- positional error metrics pos_section = {"reference": {"source": ref_source, "lat": ref[0], "lon": ref[1], "ellip_height_m": ref[2]}} - onboard_df = pvt.rename(columns={})[["utc", "lat", "lon", "height"]].copy() - for name, df in [("onboard", onboard_df)] + list(solutions.items()): + onboard_df = pvt_m.rename(columns={})[["utc", "lat", "lon", "height"]].copy() + for name, df in [("onboard", onboard_df)] + list(msol.items()): s = solution_stats(df, ref) if s: pos_section[name] = s pos_section["onboard_reported_accuracy"] = { - "hAcc_m": {"mean": round(float(pvt["hAcc"].mean()), 3), - "min": round(float(pvt["hAcc"].min()), 3), - "max": round(float(pvt["hAcc"].max()), 3)}, - "vAcc_m": {"mean": round(float(pvt["vAcc"].mean()), 3)}, - "pDOP": {"mean": round(float(pvt["pDOP"].mean()), 2), - "max": round(float(pvt["pDOP"].max()), 2)}, + "hAcc_m": {"mean": round(float(pvt_m["hAcc"].mean()), 3), + "min": round(float(pvt_m["hAcc"].min()), 3), + "max": round(float(pvt_m["hAcc"].max()), 3)}, + "vAcc_m": {"mean": round(float(pvt_m["vAcc"].mean()), 3)}, + "pDOP": {"mean": round(float(pvt_m["pDOP"].mean()), 2), + "max": round(float(pvt_m["pDOP"].max()), 2)}, } if ppk_result.get("station"): st = ppk_result["station"] pos_section["ppk_base"] = {"id": st["id"], "distance_km": round(st["km"], 1), "antenna": ppk_result.get("ant2"), "ecef_itrf2020": ppk_result.get("base_xyz")} + if gcgc_result and gcgc_result.get("station"): + st = gcgc_result["station"] + pos_section["ppk_gcgc_base"] = { + "id": st["id"], + "distance_km": None if math.isnan(st["km"]) else round(st["km"], 1), + "antenna": gcgc_result.get("ant2"), + "base_obs": gcgc_result.get("base_obs"), + "ecef": gcgc_result.get("base_xyz")} + if vrs_result and vrs_result.get("station"): + st = vrs_result["station"] + pos_section["ppk_vrs_base"] = { + "id": st["id"], "distance_km": 0.0, + "antenna": vrs_result.get("ant2"), + "ecef": vrs_result.get("base_xyz"), + "note": "virtual base synthesized at the capture position " + "(validated engine, experiments/FINDINGS.md)"} # --- timing metrics jitter_ms = (np.diff(data["rawx_hdr"]["rcvTow"].to_numpy()) - 1.0) * 1000 @@ -1140,8 +1604,9 @@ def stage_diagnostics(outdir, data, rinex_info, solutions, ppk_result, track_sou "epoch_jitter_ms": {"mean": round(float(jitter_ms.mean()), 4), "std": round(float(jitter_ms.std()), 4), "max_abs": round(float(np.abs(jitter_ms).max()), 4)}, - "tAcc_ns": {"mean": round(float(pvt["tAcc_ns"].mean()), 1), - "min": int(pvt["tAcc_ns"].min()), "max": int(pvt["tAcc_ns"].max())}, + "tAcc_ns": {"mean": round(float(pvt_m["tAcc_ns"].mean()), 1), + "min": int(pvt_m["tAcc_ns"].min()), + "max": int(pvt_m["tAcc_ns"].max())}, "leap_seconds": {"rawx": int(data["rawx_hdr"]["leapS"].iloc[0]), "valid_flag": bool(data["rawx_hdr"]["leapSecValid"].iloc[0])}, } @@ -1149,9 +1614,9 @@ def stage_diagnostics(outdir, data, rinex_info, solutions, ppk_result, track_sou # --- validation of computed solutions against the receiver's own embedded # position solution (NAV-PVT), matched epoch-by-epoch validation = {} - on = pvt[["utc", "lat", "lon", "height"]].copy() + on = pvt_m[["utc", "lat", "lon", "height"]].copy() on["utc_s"] = on["utc"].dt.round("1s") - for name, df in solutions.items(): + for name, df in msol.items(): if df is None or not len(df): continue d = df.copy() @@ -1206,6 +1671,8 @@ def stage_diagnostics(outdir, data, rinex_info, solutions, ppk_result, track_sou "file": str(args.capture), "motion_mode": "static (explicitly flagged)" if args.static else "kinematic (default)", + "quality_metrics_exclude_first_min": (args.warmup_min + if warmup_cut is not None else None), "receiver": data["meta"].get("extensions", []), "module_fw": data["meta"].get("swVersion"), "start_utc": data["t_start"].isoformat(), @@ -1231,49 +1698,81 @@ def stage_diagnostics(outdir, data, rinex_info, solutions, ppk_result, track_sou def write_report(outdir, d, figs=()): L = [] s = d["session"] + pe = d["position_error"] L.append(f"# GPS Session Report - {Path(s['file']).name}\n") L.append(f"- **Receiver**: {s.get('module_fw')} " f"({', '.join(x for x in (s.get('receiver') or []) if 'MOD=' in x)})") - L.append(f"- **Window**: {s['start_utc']} → {s['end_utc']} ({s['duration_min']} min, " + L.append(f"- **Window**: {s['start_utc']} -> {s['end_utc']} ({s['duration_min']} min, " f"GPS week {s['gps_week']})") - L.append(f"- **Track solution**: `{s['track_solution_source']}`\n") + L.append(f"- **Track solution**: {sol_disp(s['track_solution_source'])} " + f"(`{s['track_solution_source']}`)\n") - ref = d["position_error"]["reference"] + ref = pe["reference"] L.append("## Position summary\n") + wm = s.get("quality_metrics_exclude_first_min") + if wm: + L.append(f"> **Note:** the first {wm:g} minutes of the session (receiver " + "warm-up) are excluded from all metrics and figures below. " + "track.csv retains the full session.\n") + present = [n for n in SOLUTION_ORDER if pe.get(n)] + L.append("**Solution methods**\n") + for term, desc in method_definitions(present): + L.append(f"- **{term}**: {desc}") + L.append("") + ref_src = ref["source"] + for k in sorted(SOLUTION_DISPLAY, key=len, reverse=True): + ref_src = ref_src.replace(k, SOLUTION_DISPLAY[k]) L.append(f"Motion mode: **{s.get('motion_mode', '?')}**\n") - L.append(f"Reference ({ref['source']}): **{ref['lat']:.9f}°, {ref['lon']:.9f}°**, " - f"h = {ref['ellip_height_m']:.3f} m\n") + L.append(f"Reference ({ref_src}): **{ref['lat']:.9f} deg, {ref['lon']:.9f} deg**, " + f"ellipsoidal height {ref['ellip_height_m']:.3f} m\n") if not s.get("motion_mode", "").startswith("static"): - L.append("> Kinematic mode: RMS/CEP dispersion figures include any real platform " + L.append("> Kinematic mode: dispersion figures include any real platform " "motion, not just measurement error.\n") - L.append("| Solution | Epochs | Offset from ref N/E/U (m) | RMS N/E/U (m) | CEP50 | CEP95 | Fix rate |") + L.append("| Solution | Epochs | Offset from reference, north/east/up (m) | " + "RMS dispersion, north/east/up (m) | 50% radius (m) | 95% radius (m) | " + "Ambiguity fix rate |") L.append("|---|---|---|---|---|---|---|") - for name in ("ppk_kinematic", "ppk_static", "ppp", "single", "onboard"): - st = d["position_error"].get(name) - if not st: - continue + for name in present: + st = pe[name] off = "/".join(f"{v:+.3f}" for v in st["offset_from_ref_NEU_m"]) rms = "/".join(f"{v:.3f}" for v in st["rms_NEU_m"]) fr = f"{st['fix_rate']*100:.1f}%" if "fix_rate" in st else "-" - L.append(f"| {name} | {st['epochs']} | {off} | {rms} | {st['cep50_m']:.3f} | " - f"{st['cep95_m']:.3f} | {fr} |") - oa = d["position_error"]["onboard_reported_accuracy"] - L.append(f"\nOnboard self-reported: hAcc {oa['hAcc_m']['mean']} m mean " - f"({oa['hAcc_m']['min']}-{oa['hAcc_m']['max']}), pDOP {oa['pDOP']['mean']}\n") - if "ppk_base" in d["position_error"]: - b = d["position_error"]["ppk_base"] - L.append(f"PPK base: **{b['id']}** at {b['distance_km']} km, antenna `{b['antenna']}`\n") + L.append(f"| {sol_disp(name)} | {st['epochs']} | {off} | {rms} | " + f"{st['cep50_m']:.3f} | {st['cep95_m']:.3f} | {fr} |") + oa = pe["onboard_reported_accuracy"] + L.append(f"\nOnboard self-reported horizontal accuracy estimate: mean " + f"{oa['hAcc_m']['mean']} m (range {oa['hAcc_m']['min']}-" + f"{oa['hAcc_m']['max']} m); mean position dilution of precision " + f"{oa['pDOP']['mean']}\n") + if "ppk_base" in pe: + b = pe["ppk_base"] + L.append(f"CORS base station: **{b['id']}** at {b['distance_km']} km, " + f"antenna `{b['antenna']}`\n") + if "ppk_gcgc_base" in pe: + b = pe["ppk_gcgc_base"] + dist = (f"{b['distance_km']} km" if b.get("distance_km") is not None + else "unknown baseline") + L.append(f"Local GCGC base (parallel method): **{b['id']}** at {dist}, " + f"antenna `{b['antenna']}`, observation file " + f"`{Path(str(b.get('base_obs'))).name}`\n") + if "ppk_vrs_base" in pe: + b = pe["ppk_vrs_base"] + L.append(f"Synthesized virtual reference station: **{b['id']}** at the " + f"capture position (zero baseline), antenna `{b['antenna']}`\n") if d.get("validation_vs_onboard"): - L.append("## Validation against onboard (NAV-PVT) positions\n") - L.append("Computed solutions cross-checked epoch-by-epoch against the receiver's " - "own embedded position solution:\n") - L.append("| Solution | Matched epochs | Mean offset N/E/U (m) | RMS diff N/E/U (m) | Horiz offset | Threshold | Verdict |") + L.append("## Validation against the receiver's onboard solution\n") + L.append("Each computed solution is cross-checked epoch by epoch against " + "the position solution embedded in the capture (UBX-NAV-PVT " + "messages):\n") + L.append("| Solution | Matched epochs | Mean offset, north/east/up (m) | " + "RMS difference, north/east/up (m) | Horizontal offset | " + "Threshold | Verdict |") L.append("|---|---|---|---|---|---|---|") for name, v in d["validation_vs_onboard"].items(): off = "/".join(f"{x:+.3f}" for x in v["mean_offset_NEU_m"]) rms = "/".join(f"{x:.3f}" for x in v["rms_diff_NEU_m"]) - L.append(f"| {name} | {v['matched_epochs']} | {off} | {rms} | " + L.append(f"| {sol_disp(name)} | {v['matched_epochs']} | {off} | {rms} | " f"{v['horizontal_offset_m']} m | < {v['threshold_m']} m | " f"**{v['verdict']}** |") L.append("") @@ -1285,45 +1784,58 @@ def write_report(outdir, d, figs=()): L.append(f"![{title}](figures/{name})\n") L.append("## Constellations\n") - L.append("| GNSS | Sats | Signals | Meas | C/N0 by band (mean dBHz) | Slips | SFRBX | prRes RMS |") + L.append("C/N0 is the carrier-to-noise density ratio (signal strength, " + "dB-Hz).\n") + L.append("| Constellation | Satellites | Signals | Measurements | " + "Mean C/N0 by band (dB-Hz) | Cycle-slip events | Ephemeris frames | " + "Pseudorange residual RMS (m) |") L.append("|---|---|---|---|---|---|---|---|") for name, c in d["constellations"].items(): - bands = ", ".join(f"{b}:{v['cno_mean']}" for b, v in sorted(c["bands"].items())) + bands = ", ".join(f"{b}: {v['cno_mean']}" for b, v in sorted(c["bands"].items())) L.append(f"| {name} | {c['satellites_tracked']} | {c['signals_tracked']} | " f"{c['measurements']} | {bands} | {c['cycle_slip_events']} | " f"{c['sfrbx_frames']} | {c['prRes_rms_m'] or '-'} |") L.append("\n## Ionosphere\n") for k, v in d["ionosphere"]["broadcast_klobuchar"].items(): - L.append(f"- Broadcast `{k}`: {v}") - L.append(f"- Signal iono-model census: {d['ionosphere']['ionoModel_census']}") + L.append(f"- Broadcast {IONO_DISPLAY.get(k, k)}: {v}") + L.append(f"- Ionospheric correction model used per signal: " + f"{d['ionosphere']['ionoModel_census']}") gf = d["ionosphere"]["geometry_free_per_satellite"] if gf: med = np.median([r["iono_L1_m_median"] for r in gf]) - L.append(f"- Dual-frequency geometry-free slant delay, session median across " - f"{len(gf)} satellites: **{med:.2f} m at L1** " - f"(~{np.median([r['TECU_median'] for r in gf]):.1f} TECU); " - f"per-satellite detail in diagnostics.json") + L.append(f"- Dual-frequency geometry-free slant delay, session median " + f"across {len(gf)} satellites: **{med:.2f} m at the L1 " + f"frequency** (about " + f"{np.median([r['TECU_median'] for r in gf]):.1f} total electron " + f"content units, TECU); per-satellite detail in diagnostics.json") t = d["timing"] L.append("\n## Timing\n") - L.append(f"- Epoch interval jitter: mean {t['epoch_jitter_ms']['mean']} ms, " - f"σ {t['epoch_jitter_ms']['std']} ms, max |{t['epoch_jitter_ms']['max_abs']}| ms") - L.append(f"- Receiver time accuracy estimate (tAcc): mean {t['tAcc_ns']['mean']} ns " + L.append(f"- Measurement epoch interval jitter: mean " + f"{t['epoch_jitter_ms']['mean']} ms, standard deviation " + f"{t['epoch_jitter_ms']['std']} ms, maximum " + f"|{t['epoch_jitter_ms']['max_abs']}| ms") + L.append(f"- Receiver time accuracy estimate: mean {t['tAcc_ns']['mean']} ns " f"({t['tAcc_ns']['min']}-{t['tAcc_ns']['max']} ns)") - L.append(f"- Leap seconds: {t['leap_seconds']['rawx']} (valid={t['leap_seconds']['valid_flag']})") + L.append(f"- Leap seconds: {t['leap_seconds']['rawx']} " + f"(receiver-validated: {t['leap_seconds']['valid_flag']})") - L.append("\n## RF health\n") - L.append("| Block | jamInd mean/max | jamState≠0 | noise/ms | AGC | Antenna |") + L.append("\n## Radio-frequency interference and antenna health\n") + L.append("| Radio-frequency block | Jamming indicator (0-255), mean/max | " + "Jamming alerts | Noise level (per ms) | Automatic gain control | " + "Antenna status |") L.append("|---|---|---|---|---|---|") for b, v in d["rf_health"].items(): - L.append(f"| {b} | {v['jamInd_mean']}/{v['jamInd_max']} | {v['jammingState_nonzero']} " - f"| {v['noisePerMS_mean']} | {v['agcCnt_mean']:.0f} | {v['antStatus']} |") + L.append(f"| {b} | {v['jamInd_mean']}/{v['jamInd_max']} | " + f"{v['jammingState_nonzero']} " + f"| {v['noisePerMS_mean']} | {v['agcCnt_mean']:.0f} | " + f"{v['antStatus']} |") if d["warnings"]: L.append("\n## Warnings\n") for w in d["warnings"]: - L.append(f"- ⚠ {w}") + L.append(f"- WARNING: {w}") (outdir / "report.md").write_text("\n".join(L) + "\n") @@ -1352,14 +1864,14 @@ def write_html_report(outdir, d, figs=()): val_ok = all(v["verdict"] == "PASS" for v in val.values()) if val else None cards = [ - ("Best solution", esc(best_name)), - ("CEP50", f"{best.get('cep50_m', float('nan')):.2f} m" if best else "-"), + ("Best solution", esc(SOL_SHORT.get(best_name, best_name))), + ("50% radius", f"{best.get('cep50_m', float('nan')):.2f} m" if best else "-"), ("Epochs", f"{s.get('message_counts', {}).get('NAV-PVT', best.get('epochs', 0))}"), ("Duration", f"{s['duration_min']} min"), ("Validation", "-" if val_ok is None else ("all PASS" if val_ok else "FAIL")), ] if "fix_rate" in best: - cards.append(("PPK fix rate", f"{best['fix_rate']*100:.0f}%")) + cards.append(("Ambiguity fix rate", f"{best['fix_rate']*100:.0f}%")) mod = ", ".join(x for x in (s.get("receiver") or []) if "MOD=" in x) or "u-blox" parts = [f"""