From 87ff645e84227843052277c46a17b3afec8eea4c Mon Sep 17 00:00:00 2001 From: alfieprojectsdev Date: Wed, 5 Aug 2026 22:47:33 +0800 Subject: [PATCH 01/17] docs(bernese): why 5.2 tolerated what 5.4 rejects -- the I14 model set is retired Answers the question the previous commit left open by reading her own retained processing summary, R2S251210.PRC, 56,781 lines. Her run completed with 26,172 warnings and 3 errors, and produced F1_251210.SNX regardless. In Bernese ### is a warning and *** an error; she had all three *** and the run still finished. So 5.2 pressed on through conditions 5.4 treats as fatal. The version difference is one of tolerance, not capability. The RCVOBS warnings also confirm her data was multi-GNSS -- Galileo and GLONASS PRNs on a TRIMBLE ALLOY -- so V_SATSYS=GPS meant those observations were skipped with a warning rather than being absent. Two independent blockers follow, neither a misconfiguration. The satellite tables end in 2023: hers 2023-01-31, 5.4's I14 2023-08-10, 5.4's I20 2024-09-17. AIUB no longer publishes SATELLIT_I14.SAT at all, 404 for I14 against 200 for I20. Against 2025 data the PRN-to-SVN resolution lands on stale entries and CODSPP stops on BLOCK IIR-A 044, a satellite that is present in the antenna file. And the I14 ANTEX fails 5.4's consistency validation outright: given SVN and PRN inconsistent, PRN 22 against SVN G041, file not converted. All three variants in her tree fail it. I14.ATX carries four G041 entries with different PRN mappings across epochs; the other two carry one each and fail identically. A file 5.2 consumed without complaint is invalid to 5.4. Setting V_MYATX therefore fails at PID 002 where leaving it blank reaches PID 232, so it is left blank with the reasoning recorded at the override rather than as a bare default. The recommendation is to run I20 first, explicitly as a pipeline test and not as a comparison. That establishes the chain executes end to end and is BRN-001 acceptance evidence in its own right. The finding that I14 cannot be run on 5.4 at this epoch is itself a result worth putting to Abegail, since it bears on how the LUZON series can be continued. --- docs/bernese54_luzon_reprocessing_runbook.md | 71 ++++++++++++++++++++ scripts/derive_luzon_pcf.py | 18 +++++ 2 files changed, 89 insertions(+) diff --git a/docs/bernese54_luzon_reprocessing_runbook.md b/docs/bernese54_luzon_reprocessing_runbook.md index 6f8beba..f7d34ba 100644 --- a/docs/bernese54_luzon_reprocessing_runbook.md +++ b/docs/bernese54_luzon_reprocessing_runbook.md @@ -580,6 +580,77 @@ exclusions and documenting them. Option 1 first. It is the only one that does not change what is being measured. +### 4b.6 Answered: how 5.2 tolerated it, and why 5.4 will not + +Her retained processing summary — `SAVEDISK/2025/OUT/R2S251210.PRC`, 56,781 +lines — settles the question §4b.5 posed. + +**Her run completed with 26,172 warnings and 3 errors.** + +``` +### SR RCVOBS: Satellite/system not found (x thousands) + Receiver name : TRIMBLE ALLOY + PRN : E02 / E03 / R01 + +*** SR R2RDOH : NUMBER OF SAT. (NUMLST) > MAXSAT 136 > 135 +*** PG RXOBV3: TOO MANY OUTPUT FILES DEFINED 31 defined, 30 found +``` + +`###` is a warning in Bernese, `***` an error. **She had all three `***` errors +and still produced `F1_251210.SNX`.** So 5.2 pressed on through conditions 5.4 +treats as fatal — the version difference is real, and it is one of tolerance +rather than capability. + +The `RCVOBS` warnings also confirm her data *was* multi-GNSS (Galileo and +GLONASS PRNs on a TRIMBLE ALLOY): `V_SATSYS=GPS` meant those observations were +warned about and skipped, not absent. + +#### The I14 model set is retired, and 5.4 enforces what 5.2 did not + +Two independent blockers, neither a configuration error: + +**1. Satellite tables end in 2023.** Hers 2023-01-31; 5.4's `SATELLIT_I14.SAT` +2023-08-10; `SATELLIT_I20.SAT` 2024-09-17. **AIUB no longer publishes +`SATELLIT_I14.SAT`** — `BSWUSER54/CONFIG/` returns 404 for I14, 200 for I20. +Against 2025 data the PRN→SVN resolution lands on stale entries, and CODSPP +stops on `BLOCK IIR-A 044` — a satellite that *is* in the antenna file. + +**2. The I14 ANTEX fails 5.4's consistency check.** Her run logs *"Antenna phase +center model updated with: I14.ATX"*, so `V_MYATX=I14.ATX` and ATX2PCV merged it. +On 5.4 that is rejected outright: + +``` +*** PG ATX2PCV: Given SVN and PRN inconsistent in ANTEX file. + File not converted! PRN: 22 SVN: G041 +``` + +All three variants in her tree fail it — `I14.ATX` carries four `G041` entries +with different PRN mappings across epochs; `I14-orig.ATX` and `I14_1.ATX` carry +one each and fail the same way. A file 5.2 consumed without complaint is invalid +to 5.4. + +**Setting `V_MYATX` therefore fails EARLIER (PID 002) than leaving it blank +(PID 232).** It is left blank, with that reasoning recorded at the override. + +#### What this means for the exercise + +**Reproducing her I14 numbers on 5.4 is not a configuration problem to be +solved — it runs into a retired model that 5.4's stricter validation rejects.** +The honest options: + +| Option | Cost | +|---|---| +| **Run under I20** | Tables current to 2024-09 and ANTEX consistent. Answers *"does the pipeline work"* — **not** *"does it reproduce her numbers"*, since it introduces the §1.4 frame difference the exercise exists to isolate. | +| Source updated I14 tables | They are not published. Would need another archive or a hand-repaired ANTEX. | +| Relax 5.4's validation | Not obviously possible, and it would mean processing data the software considers inconsistent. | + +**Recommendation: run I20 first, explicitly as a pipeline test rather than a +comparison.** It establishes that the 31-day chain executes end to end on this +machine, which is BRN-001 acceptance evidence in its own right. The I14 +comparison then becomes a separate question — and the finding that I14 cannot be +run on 5.4 at this epoch is itself a result worth reporting to Abegail, since it +bears on how the LUZON series can be continued at all. + ## 5. Open questions — resolvable only by running it Most of the original list closed during the 2026-08-05 configuration survey diff --git a/scripts/derive_luzon_pcf.py b/scripts/derive_luzon_pcf.py index 070bd15..694f380 100755 --- a/scripts/derive_luzon_pcf.py +++ b/scripts/derive_luzon_pcf.py @@ -35,6 +35,7 @@ V_GNSSAR GRE -> ALL all *selected* systems get AR V_SATSYS GRE -> GPS she processed GPS only -- see the note in OVERRIDES; this is why I14 works + V_MYATX (blank) -> I14.ATX ATX2PCV updates the PCV model, as hers did V_CLU 8 -> 10 her cluster size V_DEL Y -> N keep results in the campaign to inspect V_SAVOBS Y -> N ditto @@ -89,6 +90,23 @@ # So this is not a workaround for a missing file. Processing GPS only is # what she did, and it is *why* I14 remains valid against 2025 data. "V_SATSYS": "GPS", + # V_MYATX IS DELIBERATELY LEFT BLANK, though her run logs + # "Antenna phase center model updated with: I14.ATX". + # + # Setting it makes ATX2PCV (PID 002) merge the ANTEX into the phase-centre + # file. On 5.4 that FAILS for every I14 ANTEX in her tree: + # + # *** PG ATX2PCV: Given SVN and PRN inconsistent in ANTEX file. + # PRN: 22 SVN: G041 + # + # I14.ATX carries four G041 entries (different PRN mappings across epochs); + # I14-orig.ATX and I14_1.ATX carry one each and fail the same check. 5.4 + # validates SVN-to-PRN consistency that 5.2 did not, so a file her run + # consumed without complaint is rejected outright. + # + # Leaving it blank gets to PID 232 (CODSPP); setting it stops at PID 002. + # Neither reaches a solution -- see runbook §4b.6 -- but the blank form + # fails later and more informatively. "V_CLU": "10", "V_DEL": "N", "V_SAVOBS": "N", From 3de664c3bac6eb8d9a53184dc653e8cad62b00fe Mon Sep 17 00:00:00 2001 From: alfieprojectsdev Date: Wed, 5 Aug 2026 22:52:35 +0800 Subject: [PATCH 02/17] =?UTF-8?q?docs(gps3):=20session=20log=20=C2=A718=20?= =?UTF-8?q?=E2=80=94=20first=20Bernese=205.4=20run=20attempts?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Nine attempts at DOY 121 of 2025. No solution yet; the run reaches PID 232 CODSPP after importing RINEX, preparing orbits and converting 92 observation files in about 100 seconds. The plan was wrong in its premise rather than its details. Adapting the 5.2 PCF by renaming scripts and repairing WAIT lists cannot work, because the PCF file format changed between versions, and 5.4 answers a 5.2 PCF with a segmentation fault rather than a parse error. Four attempts passed before that surfaced. Section 14's conclusion that FTP_DWLD could be dropped was also wrong: the products are local in 5.2-era legacy naming and 5.4 reads long-name, so they were present and unusable at once. Two further model files had to be fetched from AIUB, including SAT_2025.CRX where the installed set stopped at 2019 -- readiness gap #6 arriving as predicted. These were reported here as needing credentials; they do not. ftp.aiub.unibe.ch is firewalled from this machine and the timeout was taken as proof of inaccessibility rather than as a reason to look for another route. The most instructive error was V_SATSYS. The override set copied V_GNSSAR=ALL on the reasoning that she resolved ambiguities across all constellations, which reads the variable backwards: V_GNSSAR selects which of the already-selected systems get resolution, while V_SATSYS selects the systems, and hers reads GPS. The run attempted GLONASS and died on a satellite absent from every I14 table. It presented as a missing file and was a selection error, and switching to I20 to "fix" it would have introduced the exact frame confound the exercise exists to isolate, with plausible-looking results. Her own retained log settles why 5.2 managed what 5.4 will not: 26,172 warnings and 3 errors, and it produced a solution regardless. The version difference is one of tolerance. Behind it lie two real blockers -- satellite tables that end in 2023 and are no longer published for I14, and an I14 ANTEX that fails 5.4's SVN-to-PRN consistency check in all three of her variants. Also records three tooling defects found by use, three Bernese environment variables that collided with script locals, and one commit pushed directly to main in breach of Rule 1. --- .../SESSION_LOG_20260729_storage.md | 101 ++++++++++++++++++ 1 file changed, 101 insertions(+) diff --git a/docs/gps3-sessions/SESSION_LOG_20260729_storage.md b/docs/gps3-sessions/SESSION_LOG_20260729_storage.md index 8b32f80..1355f04 100644 --- a/docs/gps3-sessions/SESSION_LOG_20260729_storage.md +++ b/docs/gps3-sessions/SESSION_LOG_20260729_storage.md @@ -1381,3 +1381,104 @@ file report as missing. | `main` | `19c68cf`; PRs #61–#63 still open (T420's) | | LUZON campaign | Staged, registered, PCF derived, driver fixed | | First 5.4 run | **Blocked on one OSB bias product** | + +--- + +## 18. First Bernese 5.4 run attempts — 2026-08-05 evening + +Nine attempts at DOY 121 of 2025. No solution yet, but the failures were +informative and three of them were defects in this project's own tooling. Full +detail in `docs/bernese54_luzon_reprocessing_runbook.md` §4b. + +### 18.1 How far it gets + +RINEX import, orbit preparation and observation conversion all succeed — **92 +observation files in about 100 seconds** — before stopping at PID 232 `CODSPP`. + +### 18.2 The premise was wrong, not the details + +The plan was to adapt her 5.2 `PHIVOL_REL.PCF` by renaming scripts and repairing +WAIT lists. **The PCF file format changed between 5.2 and 5.4** — fixed-column +with a ruler line versus free-form `KEY=VALUE;` — and 5.4 answers a 5.2 PCF with +a **segmentation fault**, not a parse error. Four attempts passed before that +surfaced. `scripts/derive_luzon_pcf.py` now derives from 5.4's own +`RNX2SNX.PCF` and refuses to write a file whose process rows lack `CPU=`. + +### 18.3 Products: present and unusable + +§14 concluded `FTP_DWLD` could be dropped because every product was local. They +are local **in 5.2-era legacy naming** (`igs22364.sp3.Z`); 5.4 reads long-name. +Presence had been verified, usability had not. + +Two further gaps, both fetched from AIUB: the **CODE satellite-bias product**, +without which `R2S_COP` cannot generate the `IAR_*.OSB` it treats as mandatory, +and **`SAT_2025.CRX`** — the installed set stopped at 2019, so any 2025 +processing would have hit it. This is readiness **gap #6** arriving as predicted. + +**A correction worth recording:** these were reported here as needing +credentials. They do not. `ftp.aiub.unibe.ch` is firewalled from gps3 and times +out; `www.aiub.unibe.ch` redirects to a SWITCH S3 bucket that serves everything +anonymously. A timeout was taken as proof of inaccessibility rather than as a +reason to look for another route. The user's own research corrected it. + +### 18.4 `V_SATSYS` — the most instructive error + +The override set copied `V_GNSSAR = ALL` from her PCF, reasoning that she +resolved ambiguities across all constellations. **That reads the variable +backwards.** `V_GNSSAR` selects which of the *already-selected* systems get +ambiguity resolution; **`V_SATSYS` selects the systems**, and hers reads `GPS` +where 5.4 ships `GRE`. + +The run therefore attempted GLONASS and died on **GLONASS-M 861** — launched +after I14's epoch and absent from every I14 table. **It presented as a missing +file and was a constellation-selection error.** Had it been "fixed" by switching +to I20, the comparison would have silently acquired the I14/I20 confound the +whole exercise exists to isolate, and the numbers would have looked plausible. + +Processing GPS-only is not a workaround. It is what she did, and it is why I14 +is usable against 2025 data at all. + +### 18.5 The blocker is model retirement, not configuration + +Her retained log (`R2S251210.PRC`, 56,781 lines) settles it: **26,172 warnings +and 3 errors — and it produced `F1_251210.SNX` regardless.** `###` is a warning +in Bernese and `***` an error; she had all three errors and the run finished. The +version difference is one of **tolerance**, not capability. + +Two independent blockers follow: + +- **Satellite tables end in 2023.** Hers 2023-01-31, 5.4's I14 2023-08-10, I20 + 2024-09-17. **AIUB no longer publishes `SATELLIT_I14.SAT`** (404 for I14, 200 + for I20). CODSPP stops on `BLOCK IIR-A 044` — a satellite that *is* in the + antenna file; what fails is PRN→SVN resolution against a stale table. +- **The I14 ANTEX fails 5.4's consistency check.** `*** ATX2PCV: Given SVN and + PRN inconsistent … PRN 22, SVN G041`. All three variants in her tree fail. A + file 5.2 consumed without complaint is invalid to 5.4. + +**Reproducing her I14 numbers on 5.4 is therefore not a configuration problem.** +The recommendation is to run I20 first, explicitly as a *pipeline test* rather +than a comparison — that is BRN-001 acceptance evidence in its own right — and +to put the I14 finding to Abegail, since it bears on how the LUZON series can be +continued at all. + +### 18.6 Three tooling defects, found by use + +| Defect | Consequence | +|---|---| +| `find_dangling_waits()` knew only the `WAIT=` dialect | On a 5.2 PCF it parsed **zero** PIDs and zero WAITs, reported "0 dangling", and signed off a file with four broken WAIT lists | +| `REWAIT` regex captured 4 fields where the row has 3 | Old dependency survived, replacement appended after it | +| `BPE_CAMPAIGN` as a bare name | `startBPE` tests it with `-d` **relative to CWD** — the stock Bernese drivers are silently directory-dependent | + +The first is the eighth instance of §15.5's pattern, in the function written to +prevent it. All fixed; **198 tests pass**. + +Three Bernese environment variables also collided with script variables this +session — `$SRC`, `$S` and `$P` are all exported by `LOADGPS.setvar`, and two of +them silently clobbered locals, making every source file report as missing. +**Do not use bare short names in scripts that source LOADGPS.** + +### 18.7 A process note + +One commit (`c4bc867`) was pushed **directly to `main`**, breaking Rule 1. After +the PR #66 merge the branch had been deleted and work continued on `main` +without creating a new one. The commit is sound; the route was not. From 706755927950d8e0047085864924a97091525da9 Mon Sep 17 00:00:00 2001 From: alfieprojectsdev Date: Wed, 5 Aug 2026 23:22:16 +0800 Subject: [PATCH 03/17] docs(bernese): I20 pipeline test -- the chain runs; ocean loading is the last gap I20 cleared the I14 blocker, confirming the previous commit's diagnosis. Same PCF, same data, same GPS-only selection; only the frame and antenna model changed. The 24-station configuration is the deepest yet: it reaches PID 513 in 4m12s and produces FIN_20251210.NQ0, a final ambiguity-fixed solution. Its only failure is the closing QC gate, HELMR1 reporting no redundancy, and the cause is understood. 5.4's stock RNX2SNX has no V_RX3DIR -- that variable belongs to her extended 5.2 PCF -- so the RINEX 3 fiducials were never staged and none of the 24 locals is an IGS20 reference station. There was nothing to transform against. 5.4's RNX_COP does handle RINEX 3 natively, globbing long names from the same directory, so the fix was to stage both conventions into one place rather than two. All 32 stations then staged and RXOBV3 passed with every station matching a LUZON.STA entry. That surfaced the next stop: GTOCNL cannot find ocean-loading values for ALIC. No BLQ file anywhere in her tree contains it -- all four copies are identical and cover only the local network. So either her run used the fiducials purely as datum constraints rather than processing their observations, or 5.2 warned where 5.4 errors. Given that her run completed carrying three *** errors, the second is more likely, and it is the fourth instance of 5.4 enforcing a check 5.2 did not. Consequently no configuration completes a single day cleanly, and a multi-day batch is not yet worth running: 24 stations reaches the end without datum control, 32 stations stops at baseline editing. The next task is ocean-loading coefficients for the nine fiducials, freely available from the Onsala BLQ service, merged into LUZON.BLQ. That closes the last known gap and makes datum verification and a multi-day run meaningful together. --- docs/bernese54_luzon_reprocessing_runbook.md | 62 ++++++++++++++++++++ scripts/derive_luzon_pcf.py | 17 ++++++ 2 files changed, 79 insertions(+) diff --git a/docs/bernese54_luzon_reprocessing_runbook.md b/docs/bernese54_luzon_reprocessing_runbook.md index f7d34ba..6b43b36 100644 --- a/docs/bernese54_luzon_reprocessing_runbook.md +++ b/docs/bernese54_luzon_reprocessing_runbook.md @@ -651,6 +651,68 @@ comparison then becomes a separate question — and the finding that I14 cannot run on 5.4 at this epoch is itself a result worth reporting to Abegail, since it bears on how the LUZON series can be continued at all. +### 4b.7 I20 pipeline test: the chain runs, and where each configuration stops + +**I20 cleared the I14 blocker, confirming §4b.6.** Same PCF, same data, same +GPS-only selection; only the frame/antenna triple changed. + +| Configuration | Furthest PID | Time | Result | +|---|---|---|---| +| **24 local stations** (RINEX 2 only) | **513 HELMCHK** | 4m12s | **`FIN_20251210.NQ0` produced** | +| **32 stations** (+ 9 RINEX 3 fiducials) | 322 GPSEDT | 3m22s | fails earlier | + +**The 24-station run is the deepest yet and produced a final ambiguity-fixed +solution.** Its only failure is the closing QC gate: + +``` +*** PGM HELMR1: NO REDUNDANCY. NO VERIFICATION OF SITES POSSIBLE +``` + +Cause understood: **5.4's stock `RNX2SNX.PCF` has no `V_RX3DIR`** — that variable +is specific to her extended 5.2 PCF — so the RINEX 3 fiducials were never staged, +and none of the 24 locals is an IGS20 reference station. Nothing to transform +against. + +`RNX_COP` in 5.4 *does* handle RINEX 3 natively (it globs long names from the +same `${rnxDir}`), so the fix is to stage RINEX 2 and RINEX 3 into **one** +directory rather than two. Done, with BASC/CLAV de-duplicated in favour of +RINEX 3 per §1.1a. All 32 stations then staged and **`RXOBV3` passed** — every +station matched a `LUZON.STA` entry, no hard abort. + +#### But the fiducials introduce a new stop: ocean loading + +``` +*** SR GTOCNL: OCEAN LOADING CORRECTION VALUES NOT FOUND + STATION NAME : ALIC 50137M001 + FILE NAME : ${P}/LUZON/STA/LUZON.BLQ +``` + +**No `.BLQ` anywhere in her tree contains ALIC.** All four copies of +`LUZON.BLQ` are identical (1,544 lines) and cover only the local network. So +either her run did not process the fiducials' observations — using them purely +as datum constraints — or 5.2 warned where 5.4 errors. Given §4b.6's finding +that 5.2 completed with 3 `***` errors, the second is more likely. + +**This is the same shape as everything else in §4b: a check 5.4 enforces and 5.2 +did not.** It is now the fourth instance. + +#### Consequence for scheduling + +There is **no configuration that completes a single day cleanly**: + +- 24 stations → reaches the end, produces `FIN_*.NQ0`, fails datum verification +- 32 stations → fails at baseline editing on missing ocean-loading coefficients + +A multi-day batch is therefore **not yet worth running.** What a 31-day run of +the 24-station configuration would establish is throughput and stability — real, +but modest against what is already known from one day — and every solution it +produced would lack datum control. + +**The next task is to obtain ocean-loading coefficients for the nine fiducials** +(free from the Onsala/Chalmers BLQ service) and merge them into `LUZON.BLQ`. +That closes the last known gap and makes both the datum verification and a +multi-day run meaningful at the same time. + ## 5. Open questions — resolvable only by running it Most of the original list closed during the 2026-08-05 configuration survey diff --git a/scripts/derive_luzon_pcf.py b/scripts/derive_luzon_pcf.py index 694f380..8b2bf5d 100755 --- a/scripts/derive_luzon_pcf.py +++ b/scripts/derive_luzon_pcf.py @@ -134,12 +134,29 @@ def main() -> int: ap.add_argument("--src", type=Path, default=U / "PCF" / "RNX2SNX.PCF") ap.add_argument("--dest", type=Path, default=U / "PCF" / "LUZON_DLY.PCF") ap.add_argument("--apply", action="store_true") + ap.add_argument( + "--frame", choices=("i14", "i20"), default="i14", + help="i14 reproduces her configuration but cannot run on 5.4 at this " + "epoch (runbook §4b.6): the satellite tables end in 2023 and the " + "ANTEX fails 5.4's SVN/PRN check. i20 is a PIPELINE TEST ONLY — it " + "changes the reference frame and antenna model, so its coordinates " + "are NOT comparable with her I14 solutions.", + ) args = ap.parse_args() if not args.src.is_file(): print(f"FATAL: {args.src} not found", file=sys.stderr) return 1 + if args.frame == "i20": + # Pipeline-test mode. Only the frame/model triple moves; V_SATSYS stays + # GPS so exactly one thing differs from the I14 configuration and a + # completed run tells us the I14 model set was the blocker. + OVERRIDES.update({"V_PCV": "I20", "V_REFINF": "IGS20", "V_REFPSD": "IGS20"}) + print("*** FRAME: I20 — PIPELINE TEST ONLY ***") + print(" Coordinates from this run are NOT comparable with her I14") + print(" solutions. See runbook §1.4 and §4b.6.\n") + text, changed = apply_overrides(args.src.read_text(errors="replace")) print(f"src : {args.src}\ndest: {args.dest}\n") From 984b8a0cf39ee6a5ae378a1981217f9067c13993 Mon Sep 17 00:00:00 2001 From: alfieprojectsdev Date: Thu, 6 Aug 2026 08:10:35 +0800 Subject: [PATCH 04/17] feat(bernese): prepare the OTL request and a BLQ merge tool LUZON.BLQ covers the 135 local stations and none of the seven IGS fiducials, so with those staged GPSEDT stops on missing ocean-loading values for ALIC. No BLQ anywhere in the 5.2 capture contains it, so the coefficients were never computed for those sites rather than being mislaid. They are free from the Chalmers/Onsala service, but that service delivers by email, so one manual step cannot be automated away. Prepared either side of it instead. REQUEST_fiducials.txt holds the seven stations as longitude, latitude and height, the format the request form expects, derived from the same LUZON.CRD coordinates the BPE processes so the coefficients apply to the exact points in use. Six of the seven were cross-checked against Abegail's own F1_251210.SNX SITE/ID block and agree to the arcsecond; the seventh, PIMO, does not appear in her solution. merge_blq.py merges the reply back. It refuses to add a station already present unless asked, because a duplicate BLQ block is worse than a missing one -- which of two conflicting entries Bernese honours is not evident from the file -- and it validates that each incoming block is six rows of eleven values before writing anything, then backs up. The README records that the older Holt portal has been down with a mail fault since 2024-05 and that requests now go to Barre, and asks for the chosen ocean tide model to be written down, since it is part of the provenance of every coordinate that follows. --- config/bernese/otl/README.md | 57 ++++++++ config/bernese/otl/REQUEST_fiducials.txt | 7 + scripts/merge_blq.py | 157 +++++++++++++++++++++++ 3 files changed, 221 insertions(+) create mode 100644 config/bernese/otl/README.md create mode 100644 config/bernese/otl/REQUEST_fiducials.txt create mode 100755 scripts/merge_blq.py diff --git a/config/bernese/otl/README.md b/config/bernese/otl/README.md new file mode 100644 index 0000000..8b30305 --- /dev/null +++ b/config/bernese/otl/README.md @@ -0,0 +1,57 @@ +# Ocean-loading coefficients for the IGS fiducials + +`LUZON.BLQ` covers the 135 local PHIVOLCS stations and **none of the seven IGS +fiducials** the LUZON network uses as datum control. With those staged, +`GPSEDT` stops: + +``` +*** SR GTOCNL: OCEAN LOADING CORRECTION VALUES NOT FOUND + STATION NAME : ALIC 50137M001 +``` + +No `.BLQ` anywhere in the 5.2 capture contains ALIC, so the coefficients were +never computed for those sites. They are free from the Chalmers/Onsala OTL +service — **but that service delivers by email**, so one manual step is +unavoidable. + +## What is prepared here + +`REQUEST_fiducials.txt` — the seven stations as `longitude latitude height name`, +the format the request form expects. Coordinates are taken from +`$D/REF54/LUZON.CRD`, so the coefficients are computed for exactly the points +the BPE processes. Six of the seven were cross-checked against Abegail's own +`F1_251210.SNX` `SITE/ID` block and agree to the arcsecond. + +## The manual step + +1. Open + (the older Holt portal is down — a mail-system fault since 2024-05.) +2. Paste the contents of `REQUEST_fiducials.txt` into the coordinate box. +3. Ocean tide model: **FES2014b** unless there is a reason to match an older + run. Note what you chose — the model is part of the provenance. +4. Output format: **BLQ**, with radial and horizontal components. +5. Give your email address; the reply arrives as plain text. + +## Then + +Save the reply and merge: + +```bash +scripts/merge_blq.py --blq $D/REF54/LUZON.BLQ --new onsala-reply.txt # dry run +scripts/merge_blq.py --blq $D/REF54/LUZON.BLQ --new onsala-reply.txt --apply +``` + +It refuses to add a station already present (a duplicate BLQ block is worse than +a missing one — which of two conflicting entries wins is not evident from the +file), checks every incoming block is 6 rows x 11 values, and backs up first. + +Afterwards copy the updated file into the campaign: + +```bash +cp $D/REF54/LUZON.BLQ $P/LUZON/STA/LUZON.BLQ +``` + +## Why the fiducials matter + +Without them the run has no IGS reference stations, so `HELMR1` reports +`NO REDUNDANCY` and the solution has no datum verification — see runbook §4b.7. diff --git a/config/bernese/otl/REQUEST_fiducials.txt b/config/bernese/otl/REQUEST_fiducials.txt new file mode 100644 index 0000000..818f416 --- /dev/null +++ b/config/bernese/otl/REQUEST_fiducials.txt @@ -0,0 +1,7 @@ + 130.59959711 31.82405872 314.630 AIRA + 133.88552215 -23.67010932 603.238 ALIC + 127.37448500 36.39942710 116.845 DAEJ + 131.13274477 -12.84369592 125.104 DARW + 153.97864903 24.29009788 35.645 MCIL + 121.07772737 14.63572068 95.547 PIMO + 147.36601221 -2.04323087 116.351 PNGM diff --git a/scripts/merge_blq.py b/scripts/merge_blq.py new file mode 100755 index 0000000..1b2d992 --- /dev/null +++ b/scripts/merge_blq.py @@ -0,0 +1,157 @@ +#!/usr/bin/env python3 +"""Merge ocean-loading (BLQ) station blocks into an existing Bernese BLQ file. + +WHY THIS IS NEEDED +`LUZON.BLQ` covers the 135 local PHIVOLCS stations and none of the IGS +fiducials. With the fiducials staged, `GPSEDT` stops: + + *** SR GTOCNL: OCEAN LOADING CORRECTION VALUES NOT FOUND + STATION NAME : ALIC 50137M001 + FILE NAME : ${P}/LUZON/STA/LUZON.BLQ + +No `.BLQ` anywhere in the 5.2 capture contains ALIC, so this is not a staging +mistake — the coefficients were never computed for those sites. They come from +the Chalmers/Onsala OTL service, which delivers by **email**, so acquiring them +is a manual step. This handles everything on either side of that step. + +WHY IT DOES NOT JUST APPEND +Bernese matches BLQ stations by name, and a duplicate block is worse than a +missing one: which of two conflicting entries wins is not obvious from the file. +This refuses to add a station that is already present unless `--replace` is +given, and reports rather than guessing. + +Usage: + scripts/merge_blq.py --blq --new + scripts/merge_blq.py --blq ... --new ... --apply +""" +from __future__ import annotations + +import argparse +import re +import shutil +import sys +from datetime import datetime +from pathlib import Path + +# A station block header in Scherneck/Onsala BLQ output: two leading spaces then +# the site name, e.g. " ALIC". Comment lines start with $$ and are not blocks. +_STATION_RE = re.compile(r"^\s{2}(\S+)\s*$") + + +def parse_blocks(text: str) -> dict[str, list[str]]: + """Split BLQ text into {station: [lines]} — the block plus its $$ comments. + + A block runs from its name line until the next name line, so leading $$ + comment lines that belong to a station travel with it. + """ + lines = text.splitlines() + blocks: dict[str, list[str]] = {} + current: str | None = None + buf: list[str] = [] + pending: list[str] = [] + + for ln in lines: + m = _STATION_RE.match(ln) + if m and not ln.strip().startswith("$$"): + if current: + blocks[current] = buf + current = m.group(1).upper() + buf = pending + [ln] + pending = [] + elif current is None and ln.strip().startswith("$$"): + pending.append(ln) # file header, kept out of any block + elif current is not None: + if ln.strip().startswith("$$") and _looks_like_header(ln): + pending.append(ln) # header for the NEXT station + else: + buf.append(ln) + if current: + blocks[current] = buf + return blocks + + +def _looks_like_header(line: str) -> bool: + """True for the $$ banner lines Onsala emits before each station block.""" + low = line.lower() + return any(k in low for k in ("column order", "ocean loading", "scherneck", "olfg")) + + +def main() -> int: + ap = argparse.ArgumentParser(description=__doc__.split("\n")[0]) + ap.add_argument("--blq", type=Path, required=True, help="target BLQ to merge into") + ap.add_argument("--new", type=Path, required=True, help="Onsala reply containing new blocks") + ap.add_argument("--apply", action="store_true", help="write (default: dry run)") + ap.add_argument("--replace", action="store_true", help="overwrite stations already present") + args = ap.parse_args() + + for f in (args.blq, args.new): + if not f.is_file(): + print(f"FATAL: not found: {f}", file=sys.stderr) + return 1 + + target_text = args.blq.read_text(encoding="ascii", errors="replace") + existing = parse_blocks(target_text) + incoming = parse_blocks(args.new.read_text(encoding="ascii", errors="replace")) + + print(f"target : {args.blq} ({len(existing)} stations)") + print(f"incoming : {args.new} ({len(incoming)} stations)\n") + + if not incoming: + print("No station blocks parsed from the reply. Check that it is the raw") + print("BLQ text and not an HTML page or a quoted-printable email body.") + return 1 + + to_add, clash = [], [] + for name in incoming: + (clash if name in existing else to_add).append(name) + + if to_add: + print(f"will add ({len(to_add)}): {' '.join(sorted(to_add))}") + if clash: + verb = "will REPLACE" if args.replace else "already present, SKIPPING" + print(f"{verb} ({len(clash)}): {' '.join(sorted(clash))}") + if not args.replace: + print(" (pass --replace to overwrite them)") + + # Sanity-check each incoming block: 6 numeric rows, 11 columns each. + bad = [] + for name in to_add + (clash if args.replace else []): + rows = [ln for ln in incoming[name] + if ln.strip() and not ln.strip().startswith("$$") + and not _STATION_RE.match(ln)] + numeric = [r for r in rows if len(r.split()) == 11] + if len(numeric) != 6: + bad.append(f"{name}: {len(numeric)} rows of 11 values (expected 6)") + if bad: + print("\nMALFORMED BLOCKS — refusing:") + for b in bad: + print(f" {b}") + print("\nBernese needs 3 amplitude rows then 3 phase rows, 11 values each.") + return 1 + print("\nall incoming blocks have 6 rows x 11 values") + + if not args.apply: + print("\nDry run — nothing written. Re-run with --apply.") + return 0 + + stamp = datetime.now().strftime("%Y%m%d-%H%M%S") + backup = args.blq.with_suffix(f".BLQ.bak-{stamp}") + shutil.copy2(args.blq, backup) + + out = target_text.rstrip("\n") + "\n" + for name in to_add: + out += "\n".join(incoming[name]).rstrip("\n") + "\n" + if args.replace: + for name in clash: + old = "\n".join(existing[name]) + out = out.replace(old, "\n".join(incoming[name])) + + args.blq.write_text(out, encoding="ascii") + merged = parse_blocks(args.blq.read_text(encoding="ascii", errors="replace")) + print(f"\nwrote {args.blq} ({len(existing)} -> {len(merged)} stations)") + print(f"backup: {backup}") + return 0 + + +if __name__ == "__main__": + sys.exit(main()) From b8f2649e54f07fbfd50d1f731b0b5253297a564b Mon Sep 17 00:00:00 2001 From: alfieprojectsdev Date: Thu, 6 Aug 2026 08:11:45 +0800 Subject: [PATCH 05/17] fix(bernese): BLQ parser matched the documented format, not the real one The first version matched a station line as two spaces followed by a single token, which is what the format documentation shows. Run against the actual LUZON.BLQ it found ZERO stations in a file holding 135. Real Scherneck output repeats the name on the header line, sometimes separated by a tab: " ABUY ABUY" " ALBU ALBU" Data rows also begin with two spaces, so the discriminator is the token count -- one or two, with a non-numeric first token -- rather than the indentation. Amplitude and phase rows carry eleven numeric values. Now parses all 135 stations and validates each block as six rows of eleven values before writing anything. Worth noting how it was caught: the tool was run against a file whose contents were already known, so a wrong answer was recognisable as wrong. Had it first been run against the Onsala reply, zero parsed stations would have looked like a problem with the reply. --- scripts/merge_blq.py | 35 ++++++++++++++++++++++++++++------- 1 file changed, 28 insertions(+), 7 deletions(-) diff --git a/scripts/merge_blq.py b/scripts/merge_blq.py index 1b2d992..3bbd0e1 100755 --- a/scripts/merge_blq.py +++ b/scripts/merge_blq.py @@ -33,9 +33,30 @@ from datetime import datetime from pathlib import Path -# A station block header in Scherneck/Onsala BLQ output: two leading spaces then -# the site name, e.g. " ALIC". Comment lines start with $$ and are not blocks. -_STATION_RE = re.compile(r"^\s{2}(\S+)\s*$") +# A station name line in real Scherneck/Onsala BLQ output, verified against +# LUZON.BLQ rather than against the documented example: +# +# " ABUY ABUY" name REPEATED +# " ALBU \t ALBU" and sometimes tab-separated +# +# A first attempt matched a single token and found ZERO stations in a file +# holding 135 -- the documented " NAME" form is not what the service emits. +# Data rows also begin with two spaces, so the discriminator is the token count +# (1-2) plus a non-numeric first token; amplitude and phase rows carry 11 +# numeric values. Files are CRLF, which splitlines() handles. +_NUMERIC_RE = re.compile(r"^[-.\d]") + + +def _station_name(line: str) -> str | None: + """Return the site name if this is a station header line, else None.""" + if line.strip().startswith("$$") or not line.startswith(" "): + return None + tok = line.split() + if not tok or len(tok) > 2 or _NUMERIC_RE.match(tok[0]): + return None + if len(tok) == 2 and tok[0].upper() != tok[1].upper(): + return None # two different tokens: not a name line + return tok[0].upper() def parse_blocks(text: str) -> dict[str, list[str]]: @@ -51,11 +72,11 @@ def parse_blocks(text: str) -> dict[str, list[str]]: pending: list[str] = [] for ln in lines: - m = _STATION_RE.match(ln) - if m and not ln.strip().startswith("$$"): + name = _station_name(ln) + if name: if current: blocks[current] = buf - current = m.group(1).upper() + current = name buf = pending + [ln] pending = [] elif current is None and ln.strip().startswith("$$"): @@ -118,7 +139,7 @@ def main() -> int: for name in to_add + (clash if args.replace else []): rows = [ln for ln in incoming[name] if ln.strip() and not ln.strip().startswith("$$") - and not _STATION_RE.match(ln)] + and not _station_name(ln)] numeric = [r for r in rows if len(r.split()) == 11] if len(numeric) != 6: bad.append(f"{name}: {len(numeric)} rows of 11 values (expected 6)") From acba014fa42e98dbe53ae2cfd3191bb5d1994fd0 Mon Sep 17 00:00:00 2001 From: alfieprojectsdev Date: Thu, 6 Aug 2026 09:42:03 +0800 Subject: [PATCH 06/17] docs(bernese): request OTL on FES2004 to match the existing stations Correcting the earlier recommendation of FES2014b. LUZON.BLQ records its own provenance in the header, and the 135 stations already in it were computed with FES2004, CMC NO, and the Gutenberg-Bullen Greens function. Requesting a newer model would put the seven datum fiducials on a different ocean model from the 135 stations they constrain. That inconsistency does not announce itself in the file -- the header would describe the new blocks correctly and the old ones would keep their own banner -- and it would bias the datum relative to the network in a way that is hard to trace back. A newer model is preferable only if the whole file is recomputed, which is a separate decision and a larger one. Also switches the request from geodetic to ECEF XYZ. The form accepts name X Y Z as an alternative to longitude, latitude and height, so the coordinates now go in exactly as they appear in LUZON.CRD. That removes the geodetic conversion entirely, and with it the lon/lat swap the form warns about and asks users to promise against. Station names are bare four-character codes, matching how the existing blocks are written rather than the name-plus-DOMES form the error message displays. Records the full form settings with the reason for each, the queue and why-not links for a reply that does not arrive, and the request to cite Bos and Scherneck if the results reach publication. --- config/bernese/otl/README.md | 111 ++++++++++++++++------- config/bernese/otl/REQUEST_fiducials.txt | 14 +-- 2 files changed, 86 insertions(+), 39 deletions(-) diff --git a/config/bernese/otl/README.md b/config/bernese/otl/README.md index 8b30305..58f760e 100644 --- a/config/bernese/otl/README.md +++ b/config/bernese/otl/README.md @@ -1,8 +1,8 @@ -# Ocean-loading coefficients for the IGS fiducials +# Ocean-loading coefficients for the seven IGS fiducials `LUZON.BLQ` covers the 135 local PHIVOLCS stations and **none of the seven IGS -fiducials** the LUZON network uses as datum control. With those staged, -`GPSEDT` stops: +fiducials** the network uses for datum control. With those staged, `GPSEDT` +stops: ``` *** SR GTOCNL: OCEAN LOADING CORRECTION VALUES NOT FOUND @@ -10,48 +10,95 @@ fiducials** the LUZON network uses as datum control. With those staged, ``` No `.BLQ` anywhere in the 5.2 capture contains ALIC, so the coefficients were -never computed for those sites. They are free from the Chalmers/Onsala OTL -service — **but that service delivers by email**, so one manual step is -unavoidable. +never computed for those sites. They are free from the Chalmers/Onsala service, +which delivers **by email** — so one manual step is unavoidable. Everything +either side of it is prepared here. -## What is prepared here +## ⚠ Match the existing settings — do not use the service defaults -`REQUEST_fiducials.txt` — the seven stations as `longitude latitude height name`, -the format the request form expects. Coordinates are taken from -`$D/REF54/LUZON.CRD`, so the coefficients are computed for exactly the points -the BPE processes. Six of the seven were cross-checked against Abegail's own -`F1_251210.SNX` `SITE/ID` block and agree to the arcsecond. +The 135 stations already in `LUZON.BLQ` were computed with a specific +configuration, recorded in the file's own header: -## The manual step +``` +$$ Ocean tide model: FES2004 +$$ CMC: NO (corr.tide centre of mass) +$$ Gutenberg-Bullen Greens function is used +``` -1. Open - (the older Holt portal is down — a mail-system fault since 2024-05.) -2. Paste the contents of `REQUEST_fiducials.txt` into the coordinate box. -3. Ocean tide model: **FES2014b** unless there is a reason to match an older - run. Note what you chose — the model is part of the provenance. -4. Output format: **BLQ**, with radial and horizontal components. -5. Give your email address; the reply arrives as plain text. +**Request FES2004, not a newer model.** Coefficients on a different ocean model +would put the seven datum stations on a different basis from the 135 they +constrain — a systematic inconsistency that does not announce itself in the +file and is painful to trace afterwards. A newer model is better *only* if the +whole file is recomputed, which is a separate decision. -## Then +## The request -Save the reply and merge: +Paste `REQUEST_fiducials.txt`: -```bash -scripts/merge_blq.py --blq $D/REF54/LUZON.BLQ --new onsala-reply.txt # dry run -scripts/merge_blq.py --blq $D/REF54/LUZON.BLQ --new onsala-reply.txt --apply ``` +AIRA -3530185.90159 4118797.17811 3344036.67302 +ALIC -4052052.79683 4212835.97064 -2545104.50221 +DAEJ -3120042.45801 4084614.65640 3764026.77258 +DARW -4091359.66844 4684606.39892 -1408579.03120 +MCIL -5227187.43755 2551881.35064 2607618.30816 +PIMO -3186293.50508 5286624.44981 1601158.39629 +PNGM -5367943.20000 3437431.20000 -225886.00000 +``` + +These are **ECEF XYZ taken directly from `$D/REF54/LUZON.CRD`** — the same +coordinates the BPE processes. The form accepts `name X Y Z` as an alternative +to lon/lat/height, so no geodetic conversion is performed and the +"I-never-swap-lon/lat" ambiguity does not arise. + +Names are bare 4-character codes, matching how the existing blocks are written +(` ABUY ABUY`). + +## Form settings -It refuses to add a station already present (a duplicate BLQ block is worse than -a missing one — which of two conflicting entries wins is not evident from the -file), checks every incoming block is 6 rows x 11 values, and backs up first. +| Field | Value | Why | +|---|---|---| +| Ocean tide model | **FES2004** | matches the existing 135 — see above | +| Loading phenomenon | **vertical and horizontal displacements** | not gravity/tilt | +| Greens function | **elastic (Farrell, 1972)** | Farrell 1972 uses the Gutenberg-Bullen earth model named in the header | +| CMC correction | **NO** | header records `CMC: NO` | +| Output format | **BLQ** | not HARPOS | +| Plot | NO | not needed | +| Email | your address | results arrive as plain text | -Afterwards copy the updated file into the campaign: +Portal: +(The older Holt portal has been down with a mail fault since 2024-05.) + +If nothing arrives within a few hours, check the +[queue](https://barre.oso.chalmers.se/loading/queue.php) and the +[why-not criteria](https://barre.oso.chalmers.se/loading/wncace.html). + +## Merging the reply ```bash +scripts/merge_blq.py --blq $D/REF54/LUZON.BLQ --new onsala-reply.txt # dry run +scripts/merge_blq.py --blq $D/REF54/LUZON.BLQ --new onsala-reply.txt --apply cp $D/REF54/LUZON.BLQ $P/LUZON/STA/LUZON.BLQ ``` -## Why the fiducials matter +It refuses to add a station already present (a duplicate block is worse than a +missing one — which of two conflicting entries Bernese honours is not evident +from the file), validates each incoming block as six rows of eleven values, and +backs up first. + +**Record the model actually used** in the commit message when the merged file +lands. It is part of the provenance of every coordinate computed afterwards, and +the file header will say so for the new blocks but not retrospectively for the +old ones. + +## Why this matters + +Without the fiducials there are no IGS reference stations in the solution, so +`HELMR1` reports `NO REDUNDANCY` and the result has no datum verification — +runbook §4b.7. With them, and with consistent ocean loading, the 31-day run +becomes both possible and meaningful. + +## Citation -Without them the run has no IGS reference stations, so `HELMR1` reports -`NO REDUNDANCY` and the solution has no datum verification — see runbook §4b.7. +If the results are used in published work, the service asks to be cited as +`https://barre.oso.chalmers.se/loading/l.php`, naming M.S. Bos and +H.-G. Scherneck. diff --git a/config/bernese/otl/REQUEST_fiducials.txt b/config/bernese/otl/REQUEST_fiducials.txt index 818f416..98c54f4 100644 --- a/config/bernese/otl/REQUEST_fiducials.txt +++ b/config/bernese/otl/REQUEST_fiducials.txt @@ -1,7 +1,7 @@ - 130.59959711 31.82405872 314.630 AIRA - 133.88552215 -23.67010932 603.238 ALIC - 127.37448500 36.39942710 116.845 DAEJ - 131.13274477 -12.84369592 125.104 DARW - 153.97864903 24.29009788 35.645 MCIL - 121.07772737 14.63572068 95.547 PIMO - 147.36601221 -2.04323087 116.351 PNGM +AIRA -3530185.90159 4118797.17811 3344036.67302 +ALIC -4052052.79683 4212835.97064 -2545104.50221 +DAEJ -3120042.45801 4084614.65640 3764026.77258 +DARW -4091359.66844 4684606.39892 -1408579.03120 +MCIL -5227187.43755 2551881.35064 2607618.30816 +PIMO -3186293.50508 5286624.44981 1601158.39629 +PNGM -5367943.20000 3437431.20000 -225886.00000 From f42a69d2ee46451e9dc413fec5272bd36c4c5fd5 Mon Sep 17 00:00:00 2001 From: alfieprojectsdev Date: Thu, 6 Aug 2026 14:09:48 +0800 Subject: [PATCH 07/17] feat(bernese): ocean loading merged -- DOY 121 completes end to end Sessions finished: OK: 1, Error: 0, in 5m36s. Thirty stations in FIN_20251210.SNX, the same count as Abegail's solution, with HELMCHK and COMPARF both passing. The last blocker was a BLQ formatting subtlety that took four wrong guesses before reading the source. GTOCNL reads the station name with FORMAT(2X,A4), then re-reads with FORMAT(//,2X,A10) -- the // skips TWO records, so the ten-character key comes from the THIRD line after the name and must carry the station name at columns 4-7, the match being OLNAME == OLNUMB(2:5). Onsala emits three comment lines for interpolated sites and only two for sites needing no interpolation. ALIC came back as "Complete FES2004", one line short, so the read landed on an amplitude row and the station silently failed to match while the six interpolated fiducials were fine. merge_blq.py now positions the key line at exactly the third position, padding before it rather than after. Three further defects in that tool, each found by running it rather than by reading it: It parsed a station line as two spaces and a single token, which is the documented form, and found zero stations in a file holding 135. Real Scherneck output repeats the name, sometimes tab-separated. It then required those two tokens to be equal, which held for all 135 local stations and rejected every fiducial the moment one appeared -- locals repeat the name where IGS sites carry a DOMES. And it appended to end-of-file. Bernese stops reading at the first $$ END TABLE, so the new blocks sat past the terminator: present on inspection, invisible to the software. It also rewrote CRLF as LF, because Path.read_text uses universal newlines and had already normalised them before the check ran. Blocks are now inserted before the terminator and line endings detected from raw bytes. The coefficients themselves are FES2004 with CMC NO, matching the header of the existing 135 stations rather than the newer model first recommended. Mixing ocean models would have put the seven datum fiducials on a different basis from the network they constrain. --- config/bernese/otl/onsala-reply-20260806.txt | 111 +++++++++++++++++++ scripts/merge_blq.py | 88 +++++++++++++-- 2 files changed, 192 insertions(+), 7 deletions(-) create mode 100644 config/bernese/otl/onsala-reply-20260806.txt diff --git a/config/bernese/otl/onsala-reply-20260806.txt b/config/bernese/otl/onsala-reply-20260806.txt new file mode 100644 index 0000000..1055401 --- /dev/null +++ b/config/bernese/otl/onsala-reply-20260806.txt @@ -0,0 +1,111 @@ +$$ Ocean loading displacement +$$ +$$ Calculated on barre.oso.chalmers.se using olfg/olmpp of H.-G. Scherneck +$$ +$$ Greens function: mc00egbc +$$ +$$ COLUMN ORDER: M2 S2 N2 K2 K1 O1 P1 Q1 MF MM SSA +$$ +$$ ROW ORDER: +$$ AMPLITUDES (m) +$$ RADIAL +$$ TANGENTL EW +$$ TANGENTL NS +$$ PHASES (degrees) +$$ RADIAL +$$ TANGENTL EW +$$ TANGENTL NS +$$ +$$ Displacement is defined positive in Upwards, South and West direction. +$$ The phase lag is relative to Greenwich and lags positive. The +$$ Gutenberg-Bullen Greens function is used. In the ocean tide model the +$$ deficit of tidal water mass has been corrected by subtracting a uniform +$$ layer of water with a certain phase lag globally. +$$ +$$ Complete : No interpolation of ocean model was necessary +$$ _PP : Ocean model has been interpolated near the station +$$ (PP = Post-Processing) +$$ +$$ CMC: NO (corr.tide centre of mass) +$$ +$$ Ocean tide model: FES2004 +$$ +$$ END HEADER +$$ + ALIC 50137M001 +$$ Complete FES2004 +$$ Computed by OLFG, H.-G. Scherneck, Onsala Space Observatory 2026-Aug- 6 +$$ ALIC, RADI TANG lon/lat: 133.8855 -23.6701 603.238 + .00249 .00087 .00078 .00022 .00221 .00245 .00074 .00064 .00032 .00016 .00013 + .00166 .00074 .00032 .00018 .00175 .00107 .00058 .00024 .00004 .00001 .00001 + .00034 .00012 .00012 .00003 .00089 .00071 .00029 .00018 .00013 .00008 .00006 + 144.5 -60.1 101.7 -54.8 4.7 9.2 5.4 0.0 -177.1 179.9 -179.4 + 80.4 104.3 71.0 106.0 -169.1 174.7 -169.2 157.7 -54.2 -107.3 -165.2 + -128.0 92.0 -162.1 115.0 -90.9 -115.1 -90.3 -120.6 -170.4 -175.5 -178.8 +$$ + AIRA 21742S001 +$$ FES2004_PP ID: 2026-08-06 03:48:18 +$$ Computed by OLMPP by H G Scherneck, Onsala Space Observatory, 2026 +$$ AIRA, RADI TANG lon/lat: 130.5996 31.8241 314.630 + .01512 .00633 .00309 .00154 .01127 .00876 .00370 .00175 .00012 .00008 .00008 + .00382 .00186 .00073 .00061 .00163 .00128 .00055 .00023 .00008 .00002 .00000 + .00396 .00143 .00082 .00039 .00196 .00159 .00064 .00032 .00004 .00003 .00004 + 122.3 139.6 119.3 143.0 -114.2 -130.8 -113.7 -139.2 -55.1 -129.7 -171.3 + 27.8 43.7 18.8 44.9 -166.4 169.1 -166.8 159.1 -42.4 -69.6 -117.8 + -63.3 -42.2 -76.8 -51.5 93.2 75.5 94.1 64.9 41.7 32.6 4.6 +$$ + DAEJ 23902M002 +$$ FES2004_PP ID: 2026-08-06 03:48:24 +$$ Computed by OLMPP by H G Scherneck, Onsala Space Observatory, 2026 +$$ DAEJ, RADI TANG lon/lat: 127.3745 36.3994 116.845 + .00514 .00268 .00070 .00047 .00604 .00497 .00193 .00096 .00022 .00006 .00002 + .00205 .00038 .00039 .00010 .00191 .00140 .00064 .00024 .00009 .00003 .00000 + .00190 .00077 .00044 .00024 .00181 .00143 .00059 .00028 .00002 .00002 .00002 + 66.7 99.1 86.2 120.5 -102.8 -118.5 -103.0 -131.6 -2.4 -28.0 -16.9 + 153.9 154.7 135.5 82.1 -161.2 173.8 -159.0 170.6 -36.6 -57.6 -64.4 + 0.7 10.0 -16.3 20.4 98.5 80.7 101.3 69.3 99.4 58.1 6.9 +$$ + DARW 50134M001 +$$ FES2004_PP ID: 2026-08-06 03:48:29 +$$ Computed by OLMPP by H G Scherneck, Onsala Space Observatory, 2026 +$$ DARW, RADI TANG lon/lat: 131.1328 -12.8437 125.104 + .00482 .00203 .00117 .00037 .00659 .00519 .00217 .00139 .00076 .00042 .00036 + .00153 .00069 .00037 .00008 .00362 .00209 .00117 .00045 .00006 .00002 .00000 + .00209 .00107 .00044 .00018 .00102 .00073 .00034 .00023 .00016 .00009 .00008 + 56.7 46.5 46.1 61.4 3.5 -1.5 6.1 -12.3 -170.9 -173.9 -178.4 + 121.5 139.1 96.1 114.7 -171.6 169.4 -171.1 157.7 -34.1 -56.6 -68.6 + 12.3 79.1 -20.5 83.2 -0.8 -11.4 2.6 -25.5 -168.5 -173.4 -178.3 +$$ + MCIL 21789S001 +$$ FES2004_PP ID: 2026-08-06 03:48:31 +$$ Computed by OLMPP by H G Scherneck, Onsala Space Observatory, 2026 +$$ MCIL, RADI TANG lon/lat: 153.9787 24.2901 35.645 + .00748 .00507 .00100 .00140 .01293 .00962 .00427 .00187 .00034 .00029 .00027 + .00373 .00162 .00068 .00040 .00116 .00090 .00038 .00017 .00005 .00001 .00000 + .00061 .00099 .00025 .00032 .00164 .00114 .00054 .00021 .00010 .00005 .00005 + -24.4 3.9 -29.7 1.8 -136.8 -159.3 -137.1 -166.4 -132.3 -153.9 -175.1 + -57.4 -27.3 -65.6 -31.5 93.3 86.3 93.9 78.2 -54.0 -88.9 -138.2 + 115.5 129.6 172.2 117.8 112.1 100.4 112.2 90.7 1.4 8.2 1.4 +$$ + PIMO 22003M001 +$$ FES2004_PP ID: 2026-08-06 03:48:39 +$$ Computed by OLMPP by H G Scherneck, Onsala Space Observatory, 2026 +$$ PIMO, RADI TANG lon/lat: 121.0777 14.6357 95.547 + .00976 .00349 .00211 .00078 .00746 .00833 .00256 .00171 .00074 .00046 .00042 + .00466 .00185 .00091 .00048 .00355 .00238 .00116 .00042 .00006 .00002 .00001 + .00144 .00060 .00024 .00017 .00174 .00146 .00061 .00031 .00003 .00001 .00002 + 133.6 145.4 114.7 141.7 -30.6 -54.6 -29.1 -66.2 -164.9 -169.3 -177.7 + 102.1 119.1 92.7 111.6 -136.5 -168.8 -136.7 -176.2 -56.4 -101.3 -164.2 + 92.1 128.7 90.4 129.2 -164.4 168.5 -166.3 155.9 1.7 17.4 2.7 +$$ + PNGM 51006M001 +$$ FES2004_PP ID: 2026-08-06 03:48:42 +$$ Computed by OLMPP by H G Scherneck, Onsala Space Observatory, 2026 +$$ PNGM, RADI TANG lon/lat: 147.3660 -2.0432 116.351 + .00437 .00403 .00240 .00120 .01583 .00932 .00521 .00170 .00114 .00067 .00057 + .00409 .00212 .00088 .00055 .00153 .00105 .00051 .00023 .00004 .00003 .00002 + .00054 .00032 .00023 .00014 .00068 .00063 .00023 .00014 .00007 .00005 .00003 + 63.2 -6.8 54.3 -34.4 -120.3 -135.9 -120.7 -139.7 -167.8 -172.0 -178.1 + -28.7 -15.5 -30.7 -22.3 -152.0 -178.1 -151.8 166.7 -117.0 -151.0 -174.3 + 4.8 -45.6 -20.0 -64.3 -79.0 -114.2 -79.4 -118.2 -160.8 -171.3 -178.2 +$$ END TABLE diff --git a/scripts/merge_blq.py b/scripts/merge_blq.py index 3bbd0e1..68efa2e 100755 --- a/scripts/merge_blq.py +++ b/scripts/merge_blq.py @@ -15,11 +15,23 @@ is a manual step. This handles everything on either side of that step. WHY IT DOES NOT JUST APPEND +Two reasons, and the second cost a run before it was understood. + Bernese matches BLQ stations by name, and a duplicate block is worse than a missing one: which of two conflicting entries wins is not obvious from the file. This refuses to add a station that is already present unless `--replace` is given, and reports rather than guessing. +**Bernese also stops reading at the first `$$ END TABLE`.** Appending to +end-of-file puts new blocks *after* the terminator, where they are silently +invisible — the file looks correct, the station is plainly there on inspection, +and GTOCNL still reports the coefficients missing. New blocks are therefore +inserted BEFORE that marker, and the incoming file's own header (everything +above its first station) is dropped rather than carried in mid-file. + +Line endings are preserved: these files are CRLF, and rewriting them as LF is +an unnecessary change to a file other tools read. + Usage: scripts/merge_blq.py --blq --new scripts/merge_blq.py --blq ... --new ... --apply @@ -54,8 +66,13 @@ def _station_name(line: str) -> str | None: tok = line.split() if not tok or len(tok) > 2 or _NUMERIC_RE.match(tok[0]): return None - if len(tok) == 2 and tok[0].upper() != tok[1].upper(): - return None # two different tokens: not a name line + # The second token is whatever the CRD carries in that column: local sites + # repeat the name ("ABUY ABUY"), IGS sites carry a DOMES ("ALIC 50137M001"). + # An earlier version required the two to be EQUAL, which held for all 135 + # local stations and rejected every fiducial the moment one was added -- + # a rule generalised from the only examples available at the time. + # Any 1-2 token line with a non-numeric first token is a name line; data + # rows carry 11 numeric values and $$ lines are excluded above. return tok[0].upper() @@ -134,6 +151,38 @@ def main() -> int: if not args.replace: print(" (pass --replace to overwrite them)") + # Position the key line. GTOCNL reads the name with FORMAT(2X,A4), then + # re-reads with FORMAT(//,2X,A10): the "//" SKIPS TWO RECORDS, so the 10-char + # key comes from the THIRD line after the name and must carry the station + # name at columns 4-7 (the match is OLNAME == OLNUMB(2:5)). + # + # Onsala emits a "$$ NAME, RADI TANG lon/lat: ..." line that carries it, + # but its position varies: three comment lines for interpolated sites + # ("FES2004_PP ID:", "Computed by OLMPP", "$$ NAME,") and two for sites + # needing none ("Complete FES2004", "$$ NAME,"). Padding must go BEFORE that + # line, not after, or the read lands on a blank comment and the station is + # silently unmatched. + for name in list(incoming): + blk = incoming[name] + name_idx = next((i for i, ln in enumerate(blk) if _station_name(ln)), 0) + key_idx = next( + (i for i, ln in enumerate(blk) + if ln.strip().startswith(f"$$ {name}") or ln.strip().startswith(f"$$ {name.lower()}")), + None, + ) + if key_idx is None: + print(f" WARNING {name}: no '$$ {name}, ...' key line found; " + f"GTOCNL will not match this station") + continue + want = name_idx + 3 + if key_idx < want: + blk = blk[:key_idx] + ["$$"] * (want - key_idx) + blk[key_idx:] + incoming[name] = blk + print(f" {name}: key line moved from position {key_idx - name_idx} to 3") + elif key_idx > want: + print(f" WARNING {name}: key line at position {key_idx - name_idx}, " + f"expected 3") + # Sanity-check each incoming block: 6 numeric rows, 11 columns each. bad = [] for name in to_add + (clash if args.replace else []): @@ -159,15 +208,40 @@ def main() -> int: backup = args.blq.with_suffix(f".BLQ.bak-{stamp}") shutil.copy2(args.blq, backup) - out = target_text.rstrip("\n") + "\n" + # Detect from RAW BYTES. Path.read_text() uses universal newlines and has + # already turned \r\n into \n by this point, so testing target_text for + # "\r\n" is always False and silently rewrites a CRLF file as LF. + crlf = b"\r\n" in args.blq.read_bytes() + eol = "\r\n" if crlf else "\n" + lines = target_text.splitlines() + + # Insert BEFORE the terminator: Bernese stops reading at the first + # "$$ END TABLE", so anything after it is invisible however well-formed. + end_idx = next( + (i for i, ln in enumerate(lines) if ln.strip().startswith("$$ END TABLE")), + len(lines), + ) + new_lines: list[str] = [] for name in to_add: - out += "\n".join(incoming[name]).rstrip("\n") + "\n" + # Drop the reply's file-level header AND its terminators. A stray + # "$$ END TABLE" carried in mid-file would make Bernese stop reading + # there — harmless if it lands after every station, quietly destructive + # if a later merge inserts past it. + blk = [ln.rstrip("\r") for ln in incoming[name] + if not ln.strip().startswith(("$$ END TABLE", "$$ END HEADER"))] + new_lines.extend(blk) + new_lines.append("$$") + merged_lines = lines[:end_idx] + new_lines + lines[end_idx:] + if args.replace: for name in clash: - old = "\n".join(existing[name]) - out = out.replace(old, "\n".join(incoming[name])) + old, new = existing[name], incoming[name] + i = merged_lines.index(old[0]) + merged_lines[i:i + len(old)] = [ln.rstrip("\r") for ln in new] - args.blq.write_text(out, encoding="ascii") + args.blq.write_text(eol.join(merged_lines) + eol, encoding="ascii") + print(f" inserted before line {end_idx + 1} ($$ END TABLE); line endings: " + f"{'CRLF' if crlf else 'LF'}") merged = parse_blocks(args.blq.read_text(encoding="ascii", errors="replace")) print(f"\nwrote {args.blq} ({len(existing)} -> {len(merged)} stations)") print(f"backup: {backup}") From 69a16b455753d49b2234b0ae6eb2e77448eac964 Mon Sep 17 00:00:00 2001 From: alfieprojectsdev Date: Thu, 6 Aug 2026 14:37:21 +0800 Subject: [PATCH 08/17] feat(bernese): month-long LUZON driver, own results tree MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Adds scripts/run_luzon_month.sh — 2025 DOY 121-151 through LUZON_DLY, idempotent, resumable, continue-on-error, one log per day. V_RESULT moves from the stock ${S}/RNX2SNX/$Y+0 to ${S}/LUZON/$Y+0. The stock value is shared by every RNX2SNX-derived campaign, so EXAMPLE's output would land beside LUZON's with nothing in the filename to tell them apart. DOY 121 already ran under the old value and will re-run here; five minutes buys one tree holding one configuration. DOY 139 is excluded, not skipped for convenience. Our copy of her DATAPOOL holds one RINEX2 station for that day where every neighbour holds twenty-five, yet F1_251390.SNX exists — so the observations were there when she processed and our copy of that day is short. Running it would produce a solution from the nine fiducials alone, sitting in SOL/ beside thirty proper ones and distinguishable only by opening it. The run is under I20 and is a PIPELINE TEST. I14's satellite tables end in 2023, AIUB no longer publishes them, and the I14 ANTEX fails 5.4's SVN/PRN check. Coordinates from this month are NOT comparable with the I14 series; both the script header and its final summary say so, because the output will outlive the conversation that produced it. --- scripts/derive_luzon_pcf.py | 4 + scripts/run_luzon_month.sh | 144 ++++++++++++++++++++++++++++++++++++ 2 files changed, 148 insertions(+) create mode 100755 scripts/run_luzon_month.sh diff --git a/scripts/derive_luzon_pcf.py b/scripts/derive_luzon_pcf.py index 8b2bf5d..d7ca38a 100755 --- a/scripts/derive_luzon_pcf.py +++ b/scripts/derive_luzon_pcf.py @@ -90,6 +90,10 @@ # So this is not a workaround for a missing file. Processing GPS only is # what she did, and it is *why* I14 remains valid against 2025 data. "V_SATSYS": "GPS", + # Own results tree. The stock value is ${S}/RNX2SNX/$Y+0, which every + # RNX2SNX-derived campaign shares — EXAMPLE's output would land beside + # LUZON's and neither would be identifiable later. + "V_RESULT": "${S}/LUZON/$Y+0", # V_MYATX IS DELIBERATELY LEFT BLANK, though her run logs # "Antenna phase center model updated with: I14.ATX". # diff --git a/scripts/run_luzon_month.sh b/scripts/run_luzon_month.sh new file mode 100755 index 0000000..8ab08c0 --- /dev/null +++ b/scripts/run_luzon_month.sh @@ -0,0 +1,144 @@ +#!/usr/bin/env bash +# run_luzon_month.sh — process 2025 DOY 121-151 through the LUZON_DLY PCF. +# +# Follows the pattern proven by run_pagenet_week.sh: idempotent, resumable, +# one log per day, safe to re-run after an interruption. +# +# ⚠ WHAT THIS RUN IS AND IS NOT +# It runs under **I20**, not the I14 configuration Abegail used. That is a +# deliberate choice, not an oversight: I14's satellite tables end in 2023 and +# AIUB no longer publishes them, and the I14 ANTEX fails 5.4's SVN/PRN +# consistency check (runbook §4b.6). So: +# +# * This IS a pipeline test — does the chain execute end to end, unattended, +# over a month of real data. That is BRN-001 acceptance evidence. +# * This is NOT a comparison with her results. The coordinates carry an +# I14→I20 frame and antenna-model difference of centimetre scale, +# concentrated in the Up component (runbook §1.4). Do not difference these +# against F1_25*.SNX and report the residual as a Bernese-version effect. +# +# IDEMPOTENCE: a day whose FIN_20250.NQ0 already exists in SAVEDISK is +# skipped. Re-running after a crash resumes rather than redoing. +# +# CONTINUE ON ERROR: a failed day is recorded and the run moves on. For a first +# month-long pass, collecting 30 good days plus one diagnosable failure beats +# halting at hour two and learning nothing about the rest. +set -uo pipefail + +YEAR=2025 +DOY_FROM=121 +DOY_TO=151 +PCF=LUZON_DLY + +# Days excluded for want of data, NOT for convenience. Each needs a reason. +# +# 139 — our copy of her DATAPOOL holds ONE RINEX2 station for this day +# (TGDN) where every neighbouring day holds 25. She did solve it +# (F1_251390.SNX exists), so the observations existed at processing +# time and our copy of that day is incomplete. Running it anyway +# would yield a solution from the nine fiducials alone — TGDN is one +# of the two stations the DOY 121 run dropped — and that file would +# sit in SOL/ beside thirty proper ones, indistinguishable without +# opening it. A missing day is honest; a degenerate one is a trap. +SKIP_DOYS="139" +LOG_DIR="$HOME/luzon-month-logs" +SUMMARY="$LOG_DIR/summary.txt" + +die() { printf '\nERROR: %s\n' "$*" >&2; exit 1; } + +# shellcheck disable=SC1090 +source "$HOME/BERN54/LOADGPS.setvar" >/dev/null 2>&1 || die "cannot source LOADGPS.setvar" +if [ -z "${P:-}" ] || [ -z "${U:-}" ] || [ -z "${S:-}" ]; then + die "P/U/S unset after LOADGPS" +fi +[ -f "$U/PCF/$PCF.PCF" ] || die "$U/PCF/$PCF.PCF not found" +[ -x "$U/SCRIPT/luzon_pcs.pl" ] || die "driver missing: $U/SCRIPT/luzon_pcs.pl" + +mkdir -p "$LOG_DIR" + +# One instance only. A second concurrent BPE on the same campaign would corrupt +# the working directories both are using. +LOCK="$HOME/.run_luzon_month.lock" +if ! mkdir "$LOCK" 2>/dev/null; then + die "another run holds $LOCK — remove it if no BPE is running" +fi +trap 'rmdir "$LOCK" 2>/dev/null' EXIT + +{ + echo "=== LUZON month run: $YEAR DOY $DOY_FROM-$DOY_TO ===" + echo "started : $(date '+%F %T')" + echo "PCF : $PCF (V_PCV=$(grep -oE '^V_PCV *= *\S+' "$U/PCF/$PCF.PCF" | awk '{print $3}' | tr -d ';'))" + echo "campaign: $P/LUZON" + echo +} | tee "$SUMMARY" + +ok=0; failed=0; skipped=0; excluded=0; failed_days="" +start_all=$(date +%s) + +for doy in $(seq "$DOY_FROM" "$DOY_TO"); do + d3=$(printf '%03d' "$doy") + sess="${d3}0" + log="$LOG_DIR/luzon-$YEAR$sess.log" + + case " $SKIP_DOYS " in + *" $doy "*) + printf ' DOY %s EXCLUDED (insufficient observations — see SKIP_DOYS)\n' \ + "$d3" | tee -a "$SUMMARY" + excluded=$((excluded + 1)) + continue + ;; + esac + + # Idempotence: SAVEDISK is where R2S_SAV archives the final solution, so a + # completed day is evident there even after R2S_DEL cleans the campaign. + # Archived name is FIN_0.SNX.gz — four-digit year, and gzipped + # by R2S_SAV. An earlier version looked for a two-digit year in a directory + # that does not exist, so every day would have re-run. + if ls "$S/LUZON/$YEAR/SOL/FIN_${YEAR}${sess}."* >/dev/null 2>&1; then + printf ' DOY %s SKIP (already done)\n' "$d3" | tee -a "$SUMMARY" + skipped=$((skipped + 1)) + continue + fi + + t0=$(date +%s) + printf ' DOY %s start %s ... ' "$d3" "$(date '+%H:%M:%S')" | tee -a "$SUMMARY" + + echo "start: $(date '+%F %T')" > "$log" + perl "$U/SCRIPT/luzon_pcs.pl" "$YEAR" "$sess" "$PCF" >> "$log" 2>&1 + rc=$? + echo "BPE_EXIT=$rc end: $(date '+%F %T')" >> "$log" + dt=$(( $(date +%s) - t0 )) + + # The BPE's own summary decides, not the exit code — consistent with how + # every other check in this project treats exit statuses. + if grep -qE 'Sessions finished: *OK: *1 +Error: *0' "$P/LUZON/BPE/$PCF.OUT" 2>/dev/null; then + printf 'OK (%02d:%02d)\n' $((dt / 60)) $((dt % 60)) | tee -a "$SUMMARY" + ok=$((ok + 1)) + else + err=$(grep -oE '[0-9]{3}_[0-9]{3} +\S+' "$P/LUZON/BPE/$PCF.OUT" 2>/dev/null \ + | tail -1 | awk '{print $2}') + printf 'FAIL (%02d:%02d) at %s — see %s\n' \ + $((dt / 60)) $((dt % 60)) "${err:-unknown}" "$(basename "$log")" | tee -a "$SUMMARY" + failed=$((failed + 1)); failed_days="$failed_days $d3" + fi +done + +{ + total=$(( $(date +%s) - start_all )) + echo + echo "=== finished $(date '+%F %T') ===" + printf ' OK %s FAILED %s SKIPPED %s EXCLUDED %s elapsed %02d:%02d:%02d\n' \ + "$ok" "$failed" "$skipped" "$excluded" \ + $((total/3600)) $(((total%3600)/60)) $((total%60)) + if [ -n "$failed_days" ]; then echo " failed days:$failed_days"; fi + if [ "$excluded" -gt 0 ]; then + echo " excluded days: $SKIP_DOYS (data gaps, not failures — see script header)" + fi + echo + echo " Solutions: \$S/LUZON/$YEAR/SOL/ (archived by R2S_SAV)" + echo + echo " REMINDER: this ran under I20. The coordinates are NOT comparable" + echo " with Abegail's I14 solutions — see runbook §1.4 and §4b.6." +} | tee -a "$SUMMARY" + +exit "$failed" From 35fbc56a0ae8a8f0a7d4fc5ef92284f4678f1894 Mon Sep 17 00:00:00 2001 From: alfieprojectsdev Date: Thu, 6 Aug 2026 14:40:39 +0800 Subject: [PATCH 09/17] fix(bernese): month driver scored every day OK from one stale file MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Two defects, both found by dry-running the loop with the BPE call stubbed out rather than by reading it. LOADGPS.setvar exports PCF, and the script set PCF=LUZON_DLY *above* the source, so the source clobbered it and the first launch died looking for `$U/PCF/$U/PCF.PCF`. This is the fourth such collision in this campaign after $SRC, $S and $P, so config names now carry a LUZON_ prefix and an assertion across the source fails loudly if a future one is added without it. Renaming the one variable would have left the trap armed. Worse: LUZON_DLY.OUT is rewritten in place each run, and the success test grepped it without checking whose run wrote it. A day whose BPE never started would be scored against the PREVIOUS day's summary and counted OK. The stubbed dry run scored all thirty days OK from one file left by DOY 121 — a full month of green with nothing executed. The summary must now be at least as new as the day's start, and a missing or stale one is a failure carrying the exit code. With the guard the same dry run reports 30 FAILED, 1 EXCLUDED. That is the defect this campaign keeps producing: a check that reports success without having inspected anything. --- scripts/run_luzon_month.sh | 86 ++++++++++++++++++++++++++------------ 1 file changed, 59 insertions(+), 27 deletions(-) diff --git a/scripts/run_luzon_month.sh b/scripts/run_luzon_month.sh index 8ab08c0..f497998 100755 --- a/scripts/run_luzon_month.sh +++ b/scripts/run_luzon_month.sh @@ -25,10 +25,18 @@ # halting at hour two and learning nothing about the rest. set -uo pipefail -YEAR=2025 -DOY_FROM=121 -DOY_TO=151 -PCF=LUZON_DLY +# NAMES ARE PREFIXED because LOADGPS.setvar exports a large set of short +# uppercase variables -- P D U C T S SRC XG XQ PCF PAN OPT SCR and more -- and +# it is sourced BELOW this block, so anything sharing a name is silently +# overwritten before first use. `PCF=LUZON_DLY` became `$U/PCF` that way and the +# first launch died looking for +# `/home/gps3/GPSUSER/PCF//home/gps3/BERN54/USER/PCF.PCF`. That was the fourth +# such collision in this campaign ($SRC, $S, $P before it), so the fix is the +# naming rule, not another rename. The assertion after the source enforces it. +LUZON_YEAR=2025 +LUZON_DOY_FROM=121 +LUZON_DOY_TO=151 +LUZON_PCF=LUZON_DLY # Days excluded for want of data, NOT for convenience. Each needs a reason. # @@ -40,21 +48,29 @@ PCF=LUZON_DLY # of the two stations the DOY 121 run dropped — and that file would # sit in SOL/ beside thirty proper ones, indistinguishable without # opening it. A missing day is honest; a degenerate one is a trap. -SKIP_DOYS="139" -LOG_DIR="$HOME/luzon-month-logs" -SUMMARY="$LOG_DIR/summary.txt" +LUZON_SKIP_DOYS="139" +LUZON_LOG_DIR="$HOME/luzon-month-logs" +LUZON_SUMMARY="$LUZON_LOG_DIR/summary.txt" die() { printf '\nERROR: %s\n' "$*" >&2; exit 1; } +# Snapshot the config so the source below cannot change it unnoticed. This +# catches the whole collision class rather than the one name that bit us. +_snap="$LUZON_YEAR|$LUZON_DOY_FROM|$LUZON_DOY_TO|$LUZON_PCF|$LUZON_SKIP_DOYS|$LUZON_LOG_DIR" + # shellcheck disable=SC1090 source "$HOME/BERN54/LOADGPS.setvar" >/dev/null 2>&1 || die "cannot source LOADGPS.setvar" + +if [ "$_snap" != "$LUZON_YEAR|$LUZON_DOY_FROM|$LUZON_DOY_TO|$LUZON_PCF|$LUZON_SKIP_DOYS|$LUZON_LOG_DIR" ]; then + die "LOADGPS.setvar overwrote a config variable — rename it with a LUZON_ prefix" +fi if [ -z "${P:-}" ] || [ -z "${U:-}" ] || [ -z "${S:-}" ]; then die "P/U/S unset after LOADGPS" fi -[ -f "$U/PCF/$PCF.PCF" ] || die "$U/PCF/$PCF.PCF not found" +[ -f "$U/PCF/$LUZON_PCF.PCF" ] || die "$U/PCF/$LUZON_PCF.PCF not found" [ -x "$U/SCRIPT/luzon_pcs.pl" ] || die "driver missing: $U/SCRIPT/luzon_pcs.pl" -mkdir -p "$LOG_DIR" +mkdir -p "$LUZON_LOG_DIR" # One instance only. A second concurrent BPE on the same campaign would corrupt # the working directories both are using. @@ -65,25 +81,25 @@ fi trap 'rmdir "$LOCK" 2>/dev/null' EXIT { - echo "=== LUZON month run: $YEAR DOY $DOY_FROM-$DOY_TO ===" + echo "=== LUZON month run: $LUZON_YEAR DOY $LUZON_DOY_FROM-$LUZON_DOY_TO ===" echo "started : $(date '+%F %T')" - echo "PCF : $PCF (V_PCV=$(grep -oE '^V_PCV *= *\S+' "$U/PCF/$PCF.PCF" | awk '{print $3}' | tr -d ';'))" + echo "PCF : $LUZON_PCF (V_PCV=$(grep -oE '^V_PCV *= *\S+' "$U/PCF/$LUZON_PCF.PCF" | awk '{print $3}' | tr -d ';'))" echo "campaign: $P/LUZON" echo -} | tee "$SUMMARY" +} | tee "$LUZON_SUMMARY" ok=0; failed=0; skipped=0; excluded=0; failed_days="" start_all=$(date +%s) -for doy in $(seq "$DOY_FROM" "$DOY_TO"); do +for doy in $(seq "$LUZON_DOY_FROM" "$LUZON_DOY_TO"); do d3=$(printf '%03d' "$doy") sess="${d3}0" - log="$LOG_DIR/luzon-$YEAR$sess.log" + log="$LUZON_LOG_DIR/luzon-$LUZON_YEAR$sess.log" - case " $SKIP_DOYS " in + case " $LUZON_SKIP_DOYS " in *" $doy "*) printf ' DOY %s EXCLUDED (insufficient observations — see SKIP_DOYS)\n' \ - "$d3" | tee -a "$SUMMARY" + "$d3" | tee -a "$LUZON_SUMMARY" excluded=$((excluded + 1)) continue ;; @@ -94,31 +110,47 @@ for doy in $(seq "$DOY_FROM" "$DOY_TO"); do # Archived name is FIN_0.SNX.gz — four-digit year, and gzipped # by R2S_SAV. An earlier version looked for a two-digit year in a directory # that does not exist, so every day would have re-run. - if ls "$S/LUZON/$YEAR/SOL/FIN_${YEAR}${sess}."* >/dev/null 2>&1; then - printf ' DOY %s SKIP (already done)\n' "$d3" | tee -a "$SUMMARY" + if ls "$S/LUZON/$LUZON_YEAR/SOL/FIN_${LUZON_YEAR}${sess}."* >/dev/null 2>&1; then + printf ' DOY %s SKIP (already done)\n' "$d3" | tee -a "$LUZON_SUMMARY" skipped=$((skipped + 1)) continue fi t0=$(date +%s) - printf ' DOY %s start %s ... ' "$d3" "$(date '+%H:%M:%S')" | tee -a "$SUMMARY" + printf ' DOY %s start %s ... ' "$d3" "$(date '+%H:%M:%S')" | tee -a "$LUZON_SUMMARY" echo "start: $(date '+%F %T')" > "$log" - perl "$U/SCRIPT/luzon_pcs.pl" "$YEAR" "$sess" "$PCF" >> "$log" 2>&1 + perl "$U/SCRIPT/luzon_pcs.pl" "$LUZON_YEAR" "$sess" "$LUZON_PCF" >> "$log" 2>&1 rc=$? echo "BPE_EXIT=$rc end: $(date '+%F %T')" >> "$log" dt=$(( $(date +%s) - t0 )) # The BPE's own summary decides, not the exit code — consistent with how # every other check in this project treats exit statuses. - if grep -qE 'Sessions finished: *OK: *1 +Error: *0' "$P/LUZON/BPE/$PCF.OUT" 2>/dev/null; then - printf 'OK (%02d:%02d)\n' $((dt / 60)) $((dt % 60)) | tee -a "$SUMMARY" + # + # But the summary must be THIS day's. LUZON_DLY.OUT is rewritten in place + # each run, so a day whose BPE never started leaves the PREVIOUS day's file + # sitting there saying "OK: 1 Error: 0" — and the day is scored OK having + # verified nothing. A dry run of this loop scored all thirty days OK from + # one stale file. So: the file must be at least as new as this day's start. + out="$P/LUZON/BPE/$LUZON_PCF.OUT" + fresh=no + if [ -f "$out" ] && [ "$(stat -c %Y "$out")" -ge "$t0" ]; then + fresh=yes + fi + + if [ "$fresh" = no ]; then + printf 'FAIL (%02d:%02d) — BPE produced no summary (exit %s); see %s\n' \ + $((dt / 60)) $((dt % 60)) "$rc" "$(basename "$log")" | tee -a "$LUZON_SUMMARY" + failed=$((failed + 1)); failed_days="$failed_days $d3" + elif grep -qE 'Sessions finished: *OK: *1 +Error: *0' "$out" 2>/dev/null; then + printf 'OK (%02d:%02d)\n' $((dt / 60)) $((dt % 60)) | tee -a "$LUZON_SUMMARY" ok=$((ok + 1)) else - err=$(grep -oE '[0-9]{3}_[0-9]{3} +\S+' "$P/LUZON/BPE/$PCF.OUT" 2>/dev/null \ + err=$(grep -oE "[0-9]{3}_[0-9]{3} +\S+" "$out" 2>/dev/null \ | tail -1 | awk '{print $2}') printf 'FAIL (%02d:%02d) at %s — see %s\n' \ - $((dt / 60)) $((dt % 60)) "${err:-unknown}" "$(basename "$log")" | tee -a "$SUMMARY" + $((dt / 60)) $((dt % 60)) "${err:-unknown}" "$(basename "$log")" | tee -a "$LUZON_SUMMARY" failed=$((failed + 1)); failed_days="$failed_days $d3" fi done @@ -132,13 +164,13 @@ done $((total/3600)) $(((total%3600)/60)) $((total%60)) if [ -n "$failed_days" ]; then echo " failed days:$failed_days"; fi if [ "$excluded" -gt 0 ]; then - echo " excluded days: $SKIP_DOYS (data gaps, not failures — see script header)" + echo " excluded days: $LUZON_SKIP_DOYS (data gaps, not failures — see script header)" fi echo - echo " Solutions: \$S/LUZON/$YEAR/SOL/ (archived by R2S_SAV)" + echo " Solutions: \$S/LUZON/$LUZON_YEAR/SOL/ (archived by R2S_SAV)" echo echo " REMINDER: this ran under I20. The coordinates are NOT comparable" echo " with Abegail's I14 solutions — see runbook §1.4 and §4b.6." -} | tee -a "$SUMMARY" +} | tee -a "$LUZON_SUMMARY" exit "$failed" From 4b2a1332ec27365a5ecac5ffa6818d3e9f479d57 Mon Sep 17 00:00:00 2001 From: alfieprojectsdev Date: Thu, 6 Aug 2026 14:43:49 +0800 Subject: [PATCH 10/17] =?UTF-8?q?docs(bernese):=20=C2=A74b.8=20=E2=80=94?= =?UTF-8?q?=20ocean=20loading=20closed,=20month=20run=20launched?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Records that the two failure modes §4b.7 left open are both resolved and that DOY 121 completed with 30 stations, matching her count. Also records a finding that is about our data holdings rather than about processing: DOY 139 in our copy of her DATAPOOL holds one RINEX2 station where neighbouring days hold 25, yet F1_251390.SNX exists in her SAVEDISK. The observations were there when she processed; our copy of that day is short. Nothing has yet checked for other such days outside this 31-day window, which is the part worth acting on. --- docs/bernese54_luzon_reprocessing_runbook.md | 35 ++++++++++++++++++++ 1 file changed, 35 insertions(+) diff --git a/docs/bernese54_luzon_reprocessing_runbook.md b/docs/bernese54_luzon_reprocessing_runbook.md index 6b43b36..e6765f8 100644 --- a/docs/bernese54_luzon_reprocessing_runbook.md +++ b/docs/bernese54_luzon_reprocessing_runbook.md @@ -713,6 +713,41 @@ produced would lack datum control. That closes the last known gap and makes both the datum verification and a multi-day run meaningful at the same time. +### 4b.8 Closed — one day complete, and the month launched (2026-08-06) + +The ocean-loading gap is closed. Coefficients for the nine fiducials came from +the Chalmers/Onsala service (FES2004, CMC:NO, Gutenberg-Bullen — chosen to match +the 135 existing stations, not for being newest) and merged into `LUZON.BLQ` via +`scripts/merge_blq.py`. **DOY 121 then completed cleanly:** `Sessions finished: +OK: 1 Error: 0`, 5m36s, 30 stations in `FIN_20251210.SNX` — the same count +Abegail's run produced — with `HELMCHK` and `COMPARF` both passing. + +So both failure modes in the table above are resolved, and the answer to "is a +multi-day batch worth running" flipped. `scripts/run_luzon_month.sh` runs +2025 DOY 121–151 and was launched on 2026-08-06. + +**Results go to `${S}/LUZON/$Y+0`, not the stock `${S}/RNX2SNX/$Y+0`.** Every +RNX2SNX-derived campaign shares the stock path, so `EXAMPLE`'s output would land +beside LUZON's with nothing in the filename to separate them. + +**DOY 139 is excluded, and this is a data-holdings finding rather than a +processing one.** Our copy of her `DATAPOOL/LUZON` holds exactly one RINEX2 +station for that day (`TGDN`) where every neighbouring day holds 25 — yet +`F1_251390.SNX` exists in her `SAVEDISK`, so the observations were present when +she processed. **Our copy of that day is short; the original was not.** Worth +raising alongside the I14 finding, because it means the transferred set is not a +complete mirror of what she worked from, and nothing else has yet checked for +other such days outside this 31-day window. + +Running it anyway would have produced a solution from the nine fiducials alone +(TGDN is one of the two stations DOY 121 dropped), sitting in `SOL/` beside +thirty proper ones and distinguishable only by opening it. + +**The month is a pipeline test under I20, not a comparison.** §4b.6 stands: I14 +cannot run on 5.4 at this epoch. Do not difference these coordinates against the +`F1_25*` series and attribute the residual to a Bernese version change — the +frame and antenna model both moved. + ## 5. Open questions — resolvable only by running it Most of the original list closed during the 2026-08-05 configuration survey From b6dfffa26ab501738a20b9091b90979a23cbc4c6 Mon Sep 17 00:00:00 2001 From: alfieprojectsdev Date: Thu, 6 Aug 2026 15:04:26 +0800 Subject: [PATCH 11/17] docs(bernese): verify the reproducible window by census, not by handover MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit §7's claim that only DOY 121-151 of 2025 has raw observations was carried forward from the transfer handover and had never been checked. It is correct: a find across both the array and the live GPSDATA tree returns exactly those 31 days and nothing else, against 365 solved days in 2025. Records the command so a successor can re-run it rather than trust this paragraph, and notes why the boundary falls there — her DATAPOOL/LUZON is a rolling staging area holding about a month, so the transfer captured a snapshot of that window, not a chosen interval. Within the window DOY 139 is short one day's stations, so the reproducible month is 30 days rather than 31. --- docs/bernese54_luzon_reprocessing_runbook.md | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/docs/bernese54_luzon_reprocessing_runbook.md b/docs/bernese54_luzon_reprocessing_runbook.md index e6765f8..bcebb6d 100644 --- a/docs/bernese54_luzon_reprocessing_runbook.md +++ b/docs/bernese54_luzon_reprocessing_runbook.md @@ -826,3 +826,22 @@ it is why `scripts/sudo/processed_transfer.sh` copies it first and alone. Capturing the raw archive from staff machines is a separate and larger piece of work, and it is the precondition for ever reprocessing the full LUZON history. + +**Confirmed by census, 2026-08-06.** The claim above was carried forward from the +transfer handover; it has now been checked directly against both trees: + +```bash +find /srv/gnss-archive /home/gps3/GPSDATA -name '????[0-3][0-9][0-9]0.25[oOdD]' \ + | sed 's|.*/....\([0-9]\{3\}\)0\.25.|\1|' | sort -u +``` + +Both return exactly DOY 121–151 and nothing else. Against 365 solved days in +2025 alone, **the reproducible fraction of that year is 8.5%**, and of the +sixteen-year series, well under 1%. + +Within the reproducible month, one day is itself short: **DOY 139 holds one +RINEX2 station where its neighbours hold 25**, though she solved it — so even our +"complete" month is 30 days, not 31 (§4b.8). Her `DATAPOOL/LUZON` is a rolling +staging area holding roughly a month, not an archive; what was transferred is a +snapshot of that window, which is why the boundary falls where it does rather +than at anything meaningful in the data. From 99c8c749440fef1da529a10a428a797923d3bbd6 Mon Sep 17 00:00:00 2001 From: alfieprojectsdev Date: Thu, 6 Aug 2026 15:18:51 +0800 Subject: [PATCH 12/17] =?UTF-8?q?docs(session):=20=C2=A719=20=E2=80=94=20o?= =?UTF-8?q?cean=20loading=20closed,=20month=20run=20launched,=20two=20find?= =?UTF-8?q?ings=20for=20Abegail?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .../SESSION_LOG_20260729_storage.md | 103 ++++++++++++++++++ 1 file changed, 103 insertions(+) diff --git a/docs/gps3-sessions/SESSION_LOG_20260729_storage.md b/docs/gps3-sessions/SESSION_LOG_20260729_storage.md index 1355f04..dfeae69 100644 --- a/docs/gps3-sessions/SESSION_LOG_20260729_storage.md +++ b/docs/gps3-sessions/SESSION_LOG_20260729_storage.md @@ -1482,3 +1482,106 @@ them silently clobbered locals, making every source file report as missing. One commit (`c4bc867`) was pushed **directly to `main`**, breaking Rule 1. After the PR #66 merge the branch had been deleted and work continued on `main` without creating a new one. The commit is sound; the route was not. + +--- + +## 19. Ocean loading closed, the month launched — 2026-08-06 + +### 19.1 What ran + +Ocean-loading coefficients for the nine fiducials arrived from the +Chalmers/Onsala service and were merged into `LUZON.BLQ`. **DOY 121 then +completed cleanly** — `Sessions finished: OK: 1 Error: 0`, 5m36s, 30 stations in +`FIN_20251210.SNX`, the same count Abegail's run produced, with `HELMCHK` and +`COMPARF` both passing. `scripts/run_luzon_month.sh` was launched over DOY +121–151 at 14:40 and is running at roughly 5m30s per day. + +The run is under **I20** and is a pipeline test, not a comparison. §18.5 stands: +I14 cannot run on 5.4 at this epoch. + +### 19.2 The BLQ merge took four attempts, each a different wrong assumption + +`scripts/merge_blq.py` was written once and corrected four times, and the +corrections are more instructive than the script: + +1. The station-name parser matched the **documented** single-token form and + found **zero** stations in a file holding 135. The service does not emit what + its documentation shows. +2. Rewritten to require two equal tokens — which held for all 135 local stations + and rejected **every fiducial**, because IGS sites carry a DOMES number in + that column instead of a repeated name. A rule generalised from the only + examples available. +3. New blocks were appended at end-of-file, i.e. **after `$$ END TABLE`**, where + Bernese never reads them. The file looked correct, the station was plainly + there on inspection, and `GTOCNL` still reported the coefficients missing. +4. Padding to position the key line was inserted **after** it rather than before, + so `GTOCNL`'s `FORMAT(//,2X,A10)` — where `//` skips *two* records — landed on + a blank comment. + +Every one of these produced a file that looked right. Three of them produced a +file that was silently wrong at read time rather than loudly wrong at parse time. + +### 19.3 The month driver scored thirty days OK without running anything + +Two defects in `run_luzon_month.sh`, **both found by dry-running the loop with +the BPE call stubbed out**, neither by reading it. + +`LOADGPS.setvar` exports `PCF`, and the script set `PCF=LUZON_DLY` *above* the +source, so the source clobbered it. The first launch died looking for +`$U/PCF/$U/PCF.PCF`. That is the **fourth** such collision after `$SRC`, `$S` and +`$P` (§18.6), so the fix is now the naming rule and not another rename: config +names carry a `LUZON_` prefix and an assertion across the source fails loudly if +a future one is added without it. + +The serious one: `LUZON_DLY.OUT` is rewritten in place each run, and the success +test grepped it **without checking whose run wrote it**. A day whose BPE never +started would be scored against the previous day's summary and counted OK. The +stubbed dry run reported **all thirty days OK from the single file DOY 121 left +behind** — a full month of green with nothing executed. The summary must now be +at least as new as the day's start; with the guard the same dry run reports 30 +FAILED, 1 EXCLUDED. + +This is the same defect the session has produced repeatedly and now for the +ninth time: **a check that reports success without having inspected anything.** +It would have been invisible in production — thirty OK lines and an empty `SOL/` +nobody opened until much later. + +Worth noting separately: `shellcheck` passed clean on a version of this script +that referenced three unset variables under `set -u`. A grep caught them. + +### 19.4 Two findings for Abegail, both about limits rather than bugs + +**The series cannot be continued under its original model.** I14's satellite +tables end in 2023, AIUB no longer publishes `SATELLIT_I14.SAT`, and the I14 +ANTEX fails 5.4's SVN/PRN consistency check (§18.5). + +**The series cannot be rebuilt from what we hold.** §7 of the runbook asserted +that only DOY 121–151 of 2025 has raw observations; that claim came from the +transfer handover and had never been checked. It is now verified by census +across both the array and the live tree — exactly those 31 days, nothing else, +against **365 solved days in 2025 alone**. The reproducible fraction is 8.5% of +that year and well under 1% of the sixteen-year series. + +The boundary falls there because her `DATAPOOL/LUZON` is a **rolling staging +area** holding about a month; the transfer captured a snapshot of it. The missing +observations were therefore probably never on the DOSTB drive, so re-transferring +will not recover them. Locating raw RINEX on staff machines is the only path, and +it is a much larger piece of work than this run. + +Within the window, **DOY 139 holds one RINEX2 station where its neighbours hold +25**, though she solved it — so our copy of that day is short and the +reproducible month is 30 days, not 31. It is excluded from the run rather than +processed into a fiducials-only solution that would sit in `SOL/` looking +legitimate. + +### 19.5 State at end of 2026-08-06 + +- Commits on `docs/luzon-i14-investigation` (PR #67): `f42a69d`, `69a16b4`, + `35fbc56`, `4b2a133`, `b6dfffa` — all pushed and verified on `origin`. +- Results land in `${S}/LUZON/$Y+0`, moved off the stock `${S}/RNX2SNX/$Y+0` + which every RNX2SNX-derived campaign shares. +- **Still open:** DOSTB unmount before unplugging + (`sudo scripts/sudo/mount_dostb.sh --umount`); no sha256 fixity on either + archive copy; T420 PRs #61–#63 now past the Rule 2 one-week limit; the leaked + `sk-ant-oat01-` token still needs revoking by a human at claude.ai; reboot + pending for kernel 6.8.0-136; iDRAC networking unconfigured. From f6b9b6e9ffc222d10ea8677b73341e965831d038 Mon Sep 17 00:00:00 2001 From: alfieprojectsdev Date: Thu, 6 Aug 2026 15:28:11 +0800 Subject: [PATCH 13/17] feat(bernese): check the solutions, not just the exit status MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Adds scripts/coord_repeatability.py — day-to-day scatter of each station about its own mean, from the FIN_*.SNX solutions a run produces. Every check this campaign has had so far asks whether the BPE finished. A BPE can finish cleanly and produce centimetre scatter, and nothing in the pipeline would have said so. Over the first eight days the run gives median N 2.9 mm, E 3.4 mm, U 7.4 mm across 31 stations, which is ordinary for daily double-difference solutions and is the first evidence that the derived PCF is not merely executable but correct. The docstring is emphatic that this is PRECISION and not accuracy: it would look just as good for a solution in the wrong reference frame, because every day would be wrong the same way. It is not evidence about I20 versus I14 and must not be quoted as such. TGDN stands out at 13/20 mm horizontal against a 3 mm median. That is the same station that is alone on DOY 139 and was dropped from DOY 121, so it reads as a station problem rather than a configuration one — if the configuration were wrong, every station would degrade together. --- scripts/coord_repeatability.py | 122 +++++++++++++++++++++++++++++++++ 1 file changed, 122 insertions(+) create mode 100755 scripts/coord_repeatability.py diff --git a/scripts/coord_repeatability.py b/scripts/coord_repeatability.py new file mode 100755 index 0000000..99a2e2c --- /dev/null +++ b/scripts/coord_repeatability.py @@ -0,0 +1,122 @@ +#!/usr/bin/env python3 +"""Day-to-day coordinate repeatability from a run's FIN_*.SNX solutions. + +WHY THIS EXISTS +`Sessions finished: OK: 1 Error: 0` says the pipeline executed. It says nothing +about whether the numbers are any good — a BPE can complete cleanly and produce +a solution with centimetre scatter. Given how many defects in this campaign were +"a check that reports success without having inspected anything", the run needed +one check that looks at the output. + +WHAT IT MEASURES, AND WHAT IT DOES NOT +Repeatability is PRECISION, not accuracy: the scatter of each station about its +own mean over the days processed. It will look excellent for a solution sitting +in the wrong reference frame, because every day is wrong in the same way. It +therefore does NOT validate the I20 frame and must not be cited as evidence that +these coordinates are comparable with the I14 series (runbook §1.4, §4b.6). + +Expect, for daily double-difference solutions: horizontal 2-4 mm, vertical +6-10 mm. Markedly worse at ONE station is usually that station — antenna, +multipath, short sessions — rather than the configuration. If the configuration +were wrong, everything would degrade together, which is what makes a single +outlier reassuring rather than alarming. + +Usage: + scripts/coord_repeatability.py # $S/LUZON/2025 by default + scripts/coord_repeatability.py '' +""" +from __future__ import annotations + +import collections +import glob +import gzip +import math +import os +import sys + +DEFAULT_PATTERN = "$S/LUZON/2025/SOL/FIN_2025*.SNX.gz" +MIN_DAYS = 3 # below this, an RMS about the mean is not worth printing + + +def read_positions(path: str) -> dict[str, tuple[float, float, float]]: + """Return {station: (X, Y, Z)} from a SINEX SOLUTION/ESTIMATE block.""" + vals: dict[str, dict[str, float]] = collections.defaultdict(dict) + in_block = False + with gzip.open(path, "rt", errors="replace") as fh: + for line in fh: + if line.startswith("+SOLUTION/ESTIMATE"): + in_block = True + continue + if line.startswith("-SOLUTION/ESTIMATE"): + break + if not in_block or line.startswith("*"): + continue + parts = line.split() + # index type code pt soln epoch unit s value stddev + if len(parts) >= 9 and parts[1] in ("STAX", "STAY", "STAZ"): + vals[parts[2]][parts[1]] = float(parts[8]) + return { + sta: (d["STAX"], d["STAY"], d["STAZ"]) + for sta, d in vals.items() + if len(d) == 3 + } + + +def to_neu(dx: float, dy: float, dz: float, lat: float, lon: float) -> tuple[float, ...]: + """Rotate an ECEF offset into local North/East/Up at the given geodetic point.""" + sl, cl = math.sin(lon), math.cos(lon) + sb, cb = math.sin(lat), math.cos(lat) + return ( + -sb * cl * dx - sb * sl * dy + cb * dz, + -sl * dx + cl * dy, + cb * cl * dx + cb * sl * dy + sb * dz, + ) + + +def main() -> int: + pattern = sys.argv[1] if len(sys.argv) > 1 else os.path.expandvars(DEFAULT_PATTERN) + files = sorted(glob.glob(pattern)) + if not files: + print(f"no solutions matched: {pattern}", file=sys.stderr) + return 1 + + # station -> {doy: (X, Y, Z)} + series: dict[str, dict[str, tuple[float, float, float]]] = collections.defaultdict(dict) + for path in files: + doy = os.path.basename(path).split("FIN_2025")[1][:3] + for sta, xyz in read_positions(path).items(): + series[sta][doy] = xyz + + days = sorted({d for byday in series.values() for d in byday}) + print(f"days: {len(days)} ({days[0]}-{days[-1]}) stations: {len(series)}\n") + print(f"{'STA':<6}{'n':>3} {'RMS N':>8}{'RMS E':>8}{'RMS U':>8} (mm about the mean)") + + rows = [] + for sta, byday in sorted(series.items()): + if len(byday) < MIN_DAYS: + continue + pos = list(byday.values()) + mean = [sum(p[i] for p in pos) / len(pos) for i in range(3)] + x, y, z = mean + lon = math.atan2(y, x) + lat = math.atan2(z, math.hypot(x, y)) + acc = [0.0, 0.0, 0.0] + for p in pos: + for i, v in enumerate(to_neu(*(p[i] - mean[i] for i in range(3)), lat, lon)): + acc[i] += v * v + rms = [math.sqrt(a / len(pos)) * 1000 for a in acc] + rows.append((sta, len(pos), rms)) + print(f"{sta:<6}{len(pos):>3} {rms[0]:8.1f}{rms[1]:8.1f}{rms[2]:8.1f}") + + if rows: + def median(i: int) -> float: + return sorted(r[2][i] for r in rows)[len(rows) // 2] + + print(f"\nmedian: N {median(0):.1f} mm E {median(1):.1f} mm U {median(2):.1f} mm") + print("\nPRECISION ONLY — says nothing about the reference frame. See the") + print("module docstring before quoting these numbers anywhere.") + return 0 + + +if __name__ == "__main__": + sys.exit(main()) From 523aa4405fba82f216546f96169298c1aaf43f70 Mon Sep 17 00:00:00 2001 From: alfieprojectsdev Date: Thu, 6 Aug 2026 15:41:02 +0800 Subject: [PATCH 14/17] =?UTF-8?q?docs(bernese):=20=C2=A74b.9/4b.10=20?= =?UTF-8?q?=E2=80=94=20a=20station=20missing=20from=20every=20solution,=20?= =?UTF-8?q?silently?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Station-for-station against her DOY 122: 30 of 31 agree. She has S01R and not PIMO; we have PIMO and not S01R. The counts match at 31, so a station-count check passes and the difference is invisible. PIMO is explained — it is one of the nine fiducials we stage and her set did not include it. S01R is not. It is absent from all ten of our solutions, present in hers, and its input IS staged, so this is a processing gap rather than a staging one. It appears in FIN_*.CRD without an estimation flag — the a priori value passed through, never solved — and no *** message anywhere in the BPE logs names it. The obvious cause does not survive checking: its RINEX samples at 15 s against everyone else's 30 s, but V_SAMPL is 180 and divides evenly by both. The section says explicitly not to record sampling as the explanation until someone traces the import step. Better an open question than a plausible wrong answer written down as fact. §4b.10 records repeatability — median N 2.9, E 3.4, U 7.4 mm over ten days — and is emphatic that this is precision and not accuracy, since a solution in the wrong frame would look identical. TGDN's 13/20 mm is explained by a 43%-length session, the only short one in the network. Also corrects the DOY 139 rationale: TGDN is dropped on DOY 121 only, not generally, so that day would yield ten stations rather than nine. Still degenerate, but the stated reasoning was wrong. --- docs/bernese54_luzon_reprocessing_runbook.md | 68 +++++++++++++++++++- 1 file changed, 65 insertions(+), 3 deletions(-) diff --git a/docs/bernese54_luzon_reprocessing_runbook.md b/docs/bernese54_luzon_reprocessing_runbook.md index bcebb6d..5600e46 100644 --- a/docs/bernese54_luzon_reprocessing_runbook.md +++ b/docs/bernese54_luzon_reprocessing_runbook.md @@ -739,9 +739,71 @@ raising alongside the I14 finding, because it means the transferred set is not a complete mirror of what she worked from, and nothing else has yet checked for other such days outside this 31-day window. -Running it anyway would have produced a solution from the nine fiducials alone -(TGDN is one of the two stations DOY 121 dropped), sitting in `SOL/` beside -thirty proper ones and distinguishable only by opening it. +Running it anyway would have produced a solution from the nine fiducials plus +TGDN — whose own session is 43% of a day (§4b.10) — sitting in `SOL/` beside +thirty proper ones and distinguishable only by opening it. Ten stations of which +nine are fiducials is not a Luzon network solution. + +*(An earlier version of this paragraph said TGDN was "one of the two stations +DOY 121 dropped" and would therefore be absent too. TGDN is dropped on DOY 121 +only; it appears in the other nine solutions. The conclusion is unchanged — +the day is degenerate either way — but the reasoning was wrong.)* + +### 4b.9 Two stations differ from her run, and neither is the frame + +Comparing DOY 122 station-for-station against her `F1_251220.SNX`: **30 of 31 +stations agree. She has `S01R` and not `PIMO`; we have `PIMO` and not `S01R`.** +The counts match at 31, which is why this went unnoticed — a station-count check +would have passed. + +`PIMO` is straightforward: it is one of the nine fiducials we stage from +RINEX3, and her fiducial set did not include it. + +**`S01R` is the one to look at.** It is absent from *all ten* of our solutions so +far, and present in hers. What is verified: + +- Its RINEX samples at **15 s** where every other station samples at 30 s + (5760 epochs against 2880). +- It is present in `LUZON.CRD`, `.STA`, `.BLQ`, `.CLU` and `.ABB`, and appears in + `FIN_*.CRD` — but **carries no estimation flag**, i.e. it is the a priori value + passed through, never solved. Estimated stations carry `G`. +- It has entries in `BSL_*.BSL` but produces **no observation files** under its + `S0` abbreviation. +- **No `***` message anywhere in the BPE logs names it.** It leaves the solution + silently. + +What is *not* established is the cause. The obvious suspect — the 15 s sampling +— does not survive contact: `V_SAMPL = 180`, and 180 divides evenly by both 15 +and 30, so decimation alone should not exclude it. **Do not record the sampling +interval as the explanation until someone has traced the import step.** + +The consequence is worth stating plainly: a station in the network contributes +to her results and nothing to ours, and the pipeline reports success either way. +That is the same defect class as §19.3 of the session log, this time in Bernese +rather than in our own tooling. + +### 4b.10 Repeatability: the solutions are good, and that is not the same as correct + +`scripts/coord_repeatability.py` over the first ten days gives **median N 2.9 mm, +E 3.4 mm, U 7.4 mm** across 31 stations — ordinary for daily double-difference +solutions, and the first evidence that the derived PCF is not merely executable +but sound. + +**This is precision, not accuracy.** A solution in the wrong reference frame +would show the same repeatability, because every day would be wrong identically. +It is not evidence about I20 versus I14 and must not be quoted as such. + +Two stations sit outside the band, both explained by their input rather than by +the configuration: + +- **TGDN** — 13.3 / 20.0 mm horizontal. Its session is **43% of a day** (1266 + epochs of 2880), the only short session in the network. It is also the only + station present on DOY 139. +- **ANTP** — 42.8 mm vertical with normal horizontals, on a full-length session. + Unexplained; older LEICA GRX1200GGPRO / LEIAT504 equipment. + +That a single station degrades while the rest hold at 3 mm is itself the useful +signal: a wrong configuration degrades everything together. **The month is a pipeline test under I20, not a comparison.** §4b.6 stands: I14 cannot run on 5.4 at this epoch. Do not difference these coordinates against the From 9f2e761fc0dd5884855556cd414a70bdcbcd5e51 Mon Sep 17 00:00:00 2001 From: alfieprojectsdev Date: Thu, 6 Aug 2026 18:29:42 +0800 Subject: [PATCH 15/17] =?UTF-8?q?docs(bernese):=20the=20month=20completed?= =?UTF-8?q?=20=E2=80=94=2030=20days,=20and=20what=20the=20numbers=20show?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 30 of 31 days processed, 0 failures, 2h47m wall. Repeatability over the full month: median N 2.8 mm, E 3.0 mm, U 10.9 mm. Horizontal held steady as the series lengthened from ten days to thirty, which is what a stable configuration looks like. The useful result is the shape of the bad days rather than their count. Scanning for stations more than 30 mm from their own mean: 25 days are completely clean, five have exactly one bad station, and only two stations are ever involved. A bad configuration degrades every station on every day; bad stations degrade themselves. TGDN is fully explained — its sessions run 112 to 1119 epochs against a 2880 full day, and its two worst days are its two shortest. The 43% figure recorded earlier came from one day and understated the variability. LGYE is NOT explained and is flagged as the open item: full 2880 epochs on every one of its bad days, yet a 200 mm excursion on DOY 137. Session length is ruled out. Recorded as open rather than attributed to a guess. --- docs/bernese54_luzon_reprocessing_runbook.md | 51 +++++++++++++++----- 1 file changed, 38 insertions(+), 13 deletions(-) diff --git a/docs/bernese54_luzon_reprocessing_runbook.md b/docs/bernese54_luzon_reprocessing_runbook.md index 5600e46..10194dd 100644 --- a/docs/bernese54_luzon_reprocessing_runbook.md +++ b/docs/bernese54_luzon_reprocessing_runbook.md @@ -784,26 +784,51 @@ rather than in our own tooling. ### 4b.10 Repeatability: the solutions are good, and that is not the same as correct -`scripts/coord_repeatability.py` over the first ten days gives **median N 2.9 mm, -E 3.4 mm, U 7.4 mm** across 31 stations — ordinary for daily double-difference -solutions, and the first evidence that the derived PCF is not merely executable -but sound. +`scripts/coord_repeatability.py` over the **completed 30 days** gives **median +N 2.8 mm, E 3.0 mm, U 10.9 mm** across 31 stations — ordinary for daily +double-difference solutions, and the first evidence that the derived PCF is not +merely executable but sound. Horizontal held steady as the series lengthened +(2.9/3.4 mm at ten days), which is what a stable configuration looks like. **This is precision, not accuracy.** A solution in the wrong reference frame would show the same repeatability, because every day would be wrong identically. It is not evidence about I20 versus I14 and must not be quoted as such. -Two stations sit outside the band, both explained by their input rather than by -the configuration: +**No day is bad network-wide.** Scanning all 30 days for stations more than +30 mm from their own mean: **25 days are completely clean**, and the other five +have **exactly one** bad station each. Only two stations are ever involved. -- **TGDN** — 13.3 / 20.0 mm horizontal. Its session is **43% of a day** (1266 - epochs of 2880), the only short session in the network. It is also the only - station present on DOY 139. -- **ANTP** — 42.8 mm vertical with normal horizontals, on a full-length session. - Unexplained; older LEICA GRX1200GGPRO / LEIAT504 equipment. +| DOY | stations >30 mm | worst | +|---|---|---| +| 124 | 1 | TGDN 66 mm | +| 137 | 1 | LGYE 200 mm | +| 138 | 1 | TGDN 67 mm | +| 140 | 1 | LGYE 35 mm | +| 151 | 1 | LGYE 111 mm | + +That distinction carries the weight. A bad *configuration* degrades every station +on every day; bad *stations* degrade themselves. This is the second pattern. + +**TGDN is fully explained by session length**, and the 43% figure quoted from +DOY 122 understated it — sessions vary enormously day to day: + +``` +DOY 123: 875 epochs DOY 124: 124 epochs (~1 hour) +DOY 137: 1119 epochs DOY 138: 112 epochs +``` + +Its two worst days are its two shortest. Nothing to fix in the pipeline. + +**LGYE is not explained.** It has a **full 2880 epochs on every one of its bad +days** (137, 140, 151), so session length is ruled out. A 200 mm excursion on +DOY 137 against a station that is otherwise unremarkable suggests ambiguity +resolution failing on that day, or a site-specific data problem. **This is open, +and it is the most worthwhile thing to look at next** — it is the only quality +finding in the month with no explanation attached. -That a single station degrades while the rest hold at 3 mm is itself the useful -signal: a wrong configuration degrades everything together. +**ANTP** — 30.1 mm vertical with normal horizontals across the month, on +full-length sessions. Elevated but never an outlier by the 30 mm horizontal test; +older LEICA GRX1200GGPRO / LEIAT504 equipment. **The month is a pipeline test under I20, not a comparison.** §4b.6 stands: I14 cannot run on 5.4 at this epoch. Do not difference these coordinates against the From fb1b58ce885e90721aae5feb9f946edb3917a036 Mon Sep 17 00:00:00 2001 From: alfieprojectsdev Date: Thu, 6 Aug 2026 18:30:25 +0800 Subject: [PATCH 16/17] =?UTF-8?q?docs:=20the=20month=20is=20done=20?= =?UTF-8?q?=E2=80=94=2030=20days,=20zero=20failures,=20and=20the=20header?= =?UTF-8?q?=20was=20lying?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Adds §19.6 (results) and §19.7 (measured capacity) to the session log. Also rewrites the runbook's status header, which still read "Status: not yet run" after the run had completed. Anyone opening that document would have taken its opening line as current. It now states plainly that the run succeeded, that it was under I20 and therefore does NOT achieve the document's own stated goal of reproducing her numbers, and that only 31 days of 2025 can be reprocessed at all. Saying so at the top matters more than saying it in §4b.6, because the person most likely to be misled is the one who reads the first screen and stops. --- docs/bernese54_luzon_reprocessing_runbook.md | 21 +++++++- .../SESSION_LOG_20260729_storage.md | 48 +++++++++++++++++++ 2 files changed, 67 insertions(+), 2 deletions(-) diff --git a/docs/bernese54_luzon_reprocessing_runbook.md b/docs/bernese54_luzon_reprocessing_runbook.md index 10194dd..89e9d8a 100644 --- a/docs/bernese54_luzon_reprocessing_runbook.md +++ b/docs/bernese54_luzon_reprocessing_runbook.md @@ -7,8 +7,25 @@ the DOSTB drive to `/srv/gnss-archive/processed/luzon-bern52/`. comparison ng results / adjustment (fine tuning) ng PCF."* Reprocess Abegail's LUZON network under 5.4, reproduce her 5.2 numbers, and only then tune. -**Status: not yet run.** Everything below about *inputs* is measured on gps3. -Everything about *execution* is untested — §5 lists what only a run will settle. +**Status as of 2026-08-06: RUN, and it completes.** 30 of the 31 days processed +end to end with **zero failures** in 2h47m (§4b.8–§4b.10). Repeatability is +median N 2.8 mm, E 3.0 mm, U 10.9 mm. + +Two things that section headings below will not tell you, so read them here: + +1. **It ran under I20, not I14.** I14 cannot run on 5.4 at this epoch — the + satellite tables end in 2023 and the ANTEX fails a consistency check 5.2 did + not perform (§4b.6). **These coordinates are therefore NOT comparable with + Abegail's `F1_25*` series**, and differencing them will show a frame and + antenna-model change, not a Bernese-version effect. The stated goal at the top + of this document — reproduce her numbers — is **not** what this run achieved. +2. **Only 31 days of 2025 can be reprocessed at all**, verified by census (§7). + She solved 365. The rest of the sixteen-year series has solutions and no + inputs. + +Sections below marked "not yet run" or "untested" predate 2026-08-06 and are +left in place because the reasoning that led to them is still worth reading; +§4b.7 onward records what actually happened. **Revised twice on 2026-08-05.** §1.1 replaces an earlier conclusion that the reference solutions did not exist — they are in `SAVEDISK/`. §1.1a replaces a diff --git a/docs/gps3-sessions/SESSION_LOG_20260729_storage.md b/docs/gps3-sessions/SESSION_LOG_20260729_storage.md index dfeae69..a61e111 100644 --- a/docs/gps3-sessions/SESSION_LOG_20260729_storage.md +++ b/docs/gps3-sessions/SESSION_LOG_20260729_storage.md @@ -1585,3 +1585,51 @@ legitimate. archive copy; T420 PRs #61–#63 now past the Rule 2 one-week limit; the leaked `sk-ant-oat01-` token still needs revoking by a human at claude.ai; reboot pending for kernel 6.8.0-136; iDRAC networking unconfigured. + +### 19.6 The month completed — 30 days, zero failures + +Finished 17:28, **2h47m** wall, **OK 30, FAILED 0, EXCLUDED 1**. Verified against +the filesystem rather than the summary: 30 `FIN_*.SNX.gz` and 30 `.NQ0.gz`, every +one readable under `gzip -t`, every one carrying ≥25 stations, no day missing +from the 121–151 range bar the excluded 139. + +**Repeatability over the full month: median N 2.8 mm, E 3.0 mm, U 10.9 mm.** +Horizontal held steady as the series grew from ten days to thirty (2.9/3.4 mm at +ten), which is what a stable configuration looks like. + +The result that carries weight is the *shape* of the bad days. Scanning all 30 +for stations more than 30 mm from their own mean: **25 days are completely +clean**, five have **exactly one** bad station, and only **two stations** are +ever involved. A bad configuration degrades every station on every day; bad +stations degrade themselves. Nothing here points at the PCF. + +**TGDN** is fully explained: its sessions run 112 to 1119 epochs against a 2880 +full day, and its two worst days are its two shortest. The "43%" recorded earlier +came from a single day and understated the variability considerably. + +**LGYE is not explained and is the open item.** Full 2880 epochs on every one of +its bad days, yet a **200 mm** excursion on DOY 137. Session length is ruled out. +It is recorded as open rather than attributed to a guess — the same discipline +applied to S01R in §4b.9. + +### 19.7 Capacity, measured rather than assumed + +Asked what a full year would cost. Measured during the run rather than +extrapolated from the per-day figure: + +- 333 s/day at 30 stations, consuming **3.94 of 24 cores** — 16% of the machine. +- A year at this configuration: **~34 h serial**, or **~7 h** at five concurrent + days. Storage is ~1.5 GB/year, so the whole sixteen-year series is ~23 GB. +- Parallelism needs **one campaign directory per worker**; `$P/LUZON` holds the + working files for whichever day is in flight, which is why the driver takes a + single-instance lock. +- Scaling to the ~135-station national network is **estimated** at 4–8× per day + and must be measured, not modelled. +- `MAXPAR` is 1000 and 135 stations need ~405 parameters, so the ceiling sits + near 330 stations — but only because `TROPEST 0` pre-eliminates the troposphere + parameters, which would otherwise dominate. Worth knowing before anyone points + this at the 417-entry CRD. + +**None of this is the binding constraint.** We hold 31 days of observations. The +compute budget for the full series is weeks, entirely tractable, and irrelevant +until the data exists. From ad44e2a7d7c28a3cbd18485cc557ed1f7b542103 Mon Sep 17 00:00:00 2001 From: alfieprojectsdev Date: Thu, 6 Aug 2026 18:37:23 +0800 Subject: [PATCH 17/17] feat(bernese): check for coherent multi-station motion, not just outliers MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit This project exists to detect earthquakes. §4b.10's quality check scanned for any single station exceeding 30 mm from its own mean and concluded "no day is bad network-wide" -- the wrong test for a seismic signal, which is several NEARBY stations moving TOGETHER, often below what would flag any one of them alone. scripts/network_coherence_scan.py checks pairs of stations within 120 km both exceeding 8 mm horizontal in the same direction. It found what the single-station scan missed: DOY 126 moved 14 stations 8-30 mm together across the whole southern-to-central Luzon cluster, invisible before because no single station individually cleared 30 mm by much. DOY 129 and 145 show smaller versions of the same pattern. Distinguishing real from artifact: a coseismic offset is a permanent step; a processing artifact specific to one day's products is a spike that reverts. DOY 126 reverts completely by DOY 127 -- a spike, not a step. Corroborated against the PHIVOLCS/USGS catalog: no earthquake recorded on any of the three flagged dates. The confirmed M4.6 near Quezon on DOY 147 produces no anomaly at the nearest stations either -- a useful negative control, since M4.6 at tens of km is below what daily static GNSS resolves, and the scan correctly stays quiet rather than manufacturing a signal from noise. The technical cause of the spikes is not identified. Orbit file size and the fiducial-fixing list were both checked and ruled out. Left open rather than assigned a plausible guess. Also closes the LGYE question raised in §4b.10 as a possible earthquake: its nearest neighbour (BLN2, 51 km) stays within 3 mm on both of LGYE's worst days, ruling out a real coseismic origin for that excursion. It remains open as a station or processing problem. --- docs/bernese54_luzon_reprocessing_runbook.md | 78 ++++++++- scripts/network_coherence_scan.py | 163 +++++++++++++++++++ 2 files changed, 232 insertions(+), 9 deletions(-) create mode 100755 scripts/network_coherence_scan.py diff --git a/docs/bernese54_luzon_reprocessing_runbook.md b/docs/bernese54_luzon_reprocessing_runbook.md index 89e9d8a..98b5898 100644 --- a/docs/bernese54_luzon_reprocessing_runbook.md +++ b/docs/bernese54_luzon_reprocessing_runbook.md @@ -811,9 +811,11 @@ merely executable but sound. Horizontal held steady as the series lengthened would show the same repeatability, because every day would be wrong identically. It is not evidence about I20 versus I14 and must not be quoted as such. -**No day is bad network-wide.** Scanning all 30 days for stations more than -30 mm from their own mean: **25 days are completely clean**, and the other five -have **exactly one** bad station each. Only two stations are ever involved. +**By a single-station threshold, no day is bad network-wide** — see §4b.11 for +why that qualifier matters and turns out to be wrong for a different, more +consequential reason. Scanning all 30 days for stations more than 30 mm from +their own mean: **25 days are completely clean**, and the other five have +**exactly one** bad station each. Only two stations are ever involved. | DOY | stations >30 mm | worst | |---|---|---| @@ -836,17 +838,75 @@ DOY 137: 1119 epochs DOY 138: 112 epochs Its two worst days are its two shortest. Nothing to fix in the pipeline. -**LGYE is not explained.** It has a **full 2880 epochs on every one of its bad -days** (137, 140, 151), so session length is ruled out. A 200 mm excursion on -DOY 137 against a station that is otherwise unremarkable suggests ambiguity -resolution failing on that day, or a site-specific data problem. **This is open, -and it is the most worthwhile thing to look at next** — it is the only quality -finding in the month with no explanation attached. +**LGYE is not explained by session length.** It has a **full 2880 epochs on +every one of its bad days** (137, 140, 151). It is, however, **ruled out as a +seismic event** by §4b.11's neighbour check: BLN2 sits 51 km away and stays +within 3 mm on both of LGYE's worst days. A real earthquake large enough to +move LGYE 200 mm would move BLN2 too. **This is open as a station/processing +problem** — ambiguity resolution failing on that day, or a site-specific issue +at LGYE — but closed as a possible earthquake. **ANTP** — 30.1 mm vertical with normal horizontals across the month, on full-length sessions. Elevated but never an outlier by the 30 mm horizontal test; older LEICA GRX1200GGPRO / LEIAT504 equipment. +### 4b.11 A single-station threshold misses the signal this project exists to +### detect — network coherence, checked properly + +§4b.10's "no day is bad network-wide" used a 30 mm **single-station** threshold. +That is the wrong test for a seismic event: a real earthquake displaces several +**nearby** stations **together**, often by amounts well under what would flag +any one of them alone. `scripts/network_coherence_scan.py` checks for that +directly — pairs of stations within 120 km both exceeding 8 mm horizontal +(≈2.5× the median repeatability) in the same direction (cosine similarity >0.5). + +**It found what the single-station scan missed. DOY 126 (2025-05-06): 14 +stations moved together, 8–30 mm, dozens of coherent pairs across the entire +southern-to-central Luzon cluster** — ALAB, ANTP, BLN2, CAC2, GUMA, GUNG, IBAZ, +MAUB, MLPA, PIMO, SAPN, TANY, TGDN, and more. Smaller versions of the same +pattern appear on DOY 129 (8 stations) and DOY 145 (13 stations). None of these +were visible in §4b.10 — no single station on DOY 126 individually cleared +30 mm by much (ANTP peaked at 29.8), so a network-wide 14-station shift hid +inside a check built to catch one bad station. + +**Distinguishing a real event from a processing artifact: is it a step or a +spike?** A coseismic offset is permanent — it persists in every subsequent +day's solution because the ground actually moved. Reading the day-by-day series +for the DOY 126 stations: DOY 125 is quiet (1–8 mm, ordinary), DOY 126 jumps to +9–30 mm across nearly the whole network — including BLN2, IBAZ and TGDN in the +north, so it is not confined to one geographic cluster — and **DOY 127 drops +straight back to 1–5 mm.** That is a spike, not a step, and a spike that +reverts completely in one day is the signature of something specific to that +day's processing, not of ground motion. + +**Corroborated against the catalog.** A web search against PHIVOLCS/USGS +reporting found **no earthquake recorded on 2025-05-06, 05-09, or 05-25** — the +three flagged dates. There **is** a confirmed M4.6 near General Nakar, Quezon on +**2025-05-27 (DOY 147)**, and checking the stations nearest that epicenter +(POLI, MAUB, GUMA, and others) on that date shows **no anomaly at all** — 0.5 to +6.8 mm, ordinary noise. That is a useful negative control: a real but small +(M4.6) event at tens of km from the nearest station is below what daily static +GNSS resolves, and the scan correctly stays quiet for it rather than +manufacturing a signal out of noise. Both halves — flagging three unexplained +network-wide days with no earthquake behind them, and staying silent for a +real one too small to see — say the method is behaving sensibly. + +**The technical cause of the DOY 126/129/145 spikes is not identified.** Two +candidates were checked and ruled out: the CODE SP3 orbit file for DOY 126 is a +normal size (no truncation), and the fiducial-fixing list in `HLM_20251260.FIX` +is identical to every ordinary day — just AIRA. Whatever produces a whole-day, +whole-network, fully-reverting shift remains open. Recorded as unexplained +rather than assigned a plausible-sounding cause, on the same principle as the +S01R and LGYE findings above. + +**What this means for anyone using this pipeline for actual event detection**: +a single-station outlier check is not sufficient and will miss a coordinated +multi-station shift unless it happens to also blow past the single-station +threshold. Any future monitoring built on this pipeline needs the coherence +check as a matter of course, not as an afterthought — and needs the step/spike +distinction made explicit, since an automated system that flags DOY 126 as +"earthquake" without checking DOY 127 would have been wrong. + **The month is a pipeline test under I20, not a comparison.** §4b.6 stands: I14 cannot run on 5.4 at this epoch. Do not difference these coordinates against the `F1_25*` series and attribute the residual to a Bernese version change — the diff --git a/scripts/network_coherence_scan.py b/scripts/network_coherence_scan.py new file mode 100755 index 0000000..f661f20 --- /dev/null +++ b/scripts/network_coherence_scan.py @@ -0,0 +1,163 @@ +#!/usr/bin/env python3 +"""Scan a run's FIN_*.SNX solutions for coherent multi-station motion. + +WHY THIS EXISTS, AND WHY coord_repeatability.py IS NOT ENOUGH +This project's purpose is earthquake monitoring, and a real coseismic offset +is not one bad station — it is SEVERAL NEARBY stations moving together, often +by amounts well under any single-station outlier threshold. The first pass at +quality-checking this run's output (coord_repeatability.py) scanned for +stations exceeding 30 mm horizontal individually and concluded "no day is bad +network-wide." That conclusion was an artifact of the threshold: on 2025 DOY +126, FOURTEEN stations moved 8-30 mm in the same direction on the same day +(runbook §4b.11), and it was invisible to a check built to catch one bad +station rather than a coordinated shift. + +WHAT DISTINGUISHES A REAL EVENT FROM A PROCESSING ARTIFACT +A real coseismic offset is a STEP: it persists in every subsequent daily +solution relative to the pre-event position, because the ground actually +moved. A processing artifact specific to one day's products or one day's +ambiguity resolution is a SPIKE: elevated on that day only, back to baseline +the next. This script flags coherent-motion days; TELLING THEM APART is a +manual read of the printed day-by-day series for the flagged stations, not +something this script decides for you. + +For 2025 DOY 121-151, DOY 126 (14 stations), 129 (8) and 145 (13) all show +this pattern and all reverted completely the next day -- read as processing +artifacts, not seismic, and corroborated by the PHIVOLCS/USGS catalog showing +no recorded event on any of those three dates. DOY 147 carries a confirmed +M4.6 near General Nakar, Quezon and shows NO anomaly at the nearest stations +(POLI, MAUB, GUMA) -- consistent with M4.6 being below what daily static GNSS +resolves at tens of km, and a useful negative control: this scan does not +manufacture events out of ordinary noise. + +WHAT THIS SCRIPT DOES NOT DO +It does not check the earthquake catalog itself, does not distinguish step +from spike automatically, and the NEIGHBOR_KM/THRESH_MM constants are tuned +for the ~135-station Luzon spacing, not validated for any other network. + +Usage: + scripts/network_coherence_scan.py # $S/LUZON/2025 by default + scripts/network_coherence_scan.py '' +""" +from __future__ import annotations + +import collections +import glob +import gzip +import math +import os +import sys + +DEFAULT_PATTERN = "$S/LUZON/2025/SOL/FIN_2025*.SNX.gz" +NEIGHBOR_KM = 120 # local network station spacing; excludes the fiducials by construction +THRESH_MM = 8.0 # roughly 2.5x median horizontal repeatability -- see coord_repeatability.py +COS_SIM_MIN = 0.5 # "same general direction", not opposite or orthogonal + + +def read_positions(path: str) -> dict[str, tuple[float, float, float]]: + vals: dict[str, dict[str, float]] = collections.defaultdict(dict) + in_block = False + with gzip.open(path, "rt", errors="replace") as fh: + for line in fh: + if line.startswith("+SOLUTION/ESTIMATE"): + in_block = True + continue + if line.startswith("-SOLUTION/ESTIMATE"): + break + if not in_block or line.startswith("*"): + continue + parts = line.split() + if len(parts) >= 9 and parts[1] in ("STAX", "STAY", "STAZ"): + vals[parts[2]][parts[1]] = float(parts[8]) + return {s: (d["STAX"], d["STAY"], d["STAZ"]) for s, d in vals.items() if len(d) == 3} + + +def to_neu(dx: float, dy: float, dz: float, lat: float, lon: float) -> tuple[float, float, float]: + sl, cl = math.sin(lon), math.cos(lon) + sb, cb = math.sin(lat), math.cos(lat) + return ( + -sb * cl * dx - sb * sl * dy + cb * dz, + -sl * dx + cl * dy, + cb * cl * dx + cb * sl * dy + sb * dz, + ) + + +def main() -> int: + pattern = sys.argv[1] if len(sys.argv) > 1 else os.path.expandvars(DEFAULT_PATTERN) + files = sorted(glob.glob(pattern)) + if not files: + print(f"no solutions matched: {pattern}", file=sys.stderr) + return 1 + + series: dict[str, dict[str, tuple[float, float, float]]] = collections.defaultdict(dict) + for path in files: + doy = os.path.basename(path).split("FIN_2025")[1][:3] + for sta, xyz in read_positions(path).items(): + series[sta][doy] = xyz + + mean = { + sta: tuple(sum(p[i] for p in byday.values()) / len(byday) for i in range(3)) + for sta, byday in series.items() + } + + def dist_km(a: str, b: str) -> float: + return math.dist(mean[a], mean[b]) / 1000 + + def offsets_for_day(doy: str) -> dict[str, tuple[float, float, float]]: + out = {} + for sta, byday in series.items(): + if doy not in byday: + continue + x, y, z = byday[doy] + mx, my, mz = mean[sta] + lon = math.atan2(my, mx) + lat = math.atan2(mz, math.hypot(mx, my)) + out[sta] = to_neu(x - mx, y - my, z - mz, lat, lon) + return out + + days = sorted({d for v in series.values() for d in v}) + print(f"{len(days)} days, {len(series)} stations, threshold {THRESH_MM} mm, " + f"neighbor radius {NEIGHBOR_KM} km\n") + + any_flagged = False + for doy in days: + offs = offsets_for_day(doy) + movers = { + s: (n * 1000, e * 1000) + for s, (n, e, u) in offs.items() + if math.hypot(n, e) * 1000 > THRESH_MM + } + pairs = [] + names = sorted(movers) + for i, a in enumerate(names): + for b in names[i + 1:]: + if dist_km(a, b) > NEIGHBOR_KM: + continue + na, ea = movers[a] + nb, eb = movers[b] + mag_a, mag_b = math.hypot(na, ea), math.hypot(nb, eb) + cos_sim = (na * nb + ea * eb) / (mag_a * mag_b) + if cos_sim > COS_SIM_MIN: + pairs.append((a, b, dist_km(a, b), cos_sim, mag_a, mag_b)) + if pairs: + any_flagged = True + stations_involved = sorted({s for a, b, *_ in pairs for s in (a, b)}) + print(f"DOY {doy}: {len(stations_involved)} station(s) moving coherently " + f"({', '.join(stations_involved)})") + for a, b, d, c, ma, mb in sorted(pairs, key=lambda r: -r[3]): + print(f" {a}({ma:.0f}mm) <-> {b}({mb:.0f}mm) {d:.0f} km apart, " + f"cos_sim={c:.2f}") + print() + + if not any_flagged: + print("No day shows two or more nearby stations moving coherently above threshold.") + + print("A flagged day is a SPIKE (processing artifact) if it reverts the next day, and a") + print("STEP (real, worth investigating) if it persists. This script does not tell them") + print("apart -- read the flagged stations' full day-by-day series before concluding either") + print("way, and check the PHIVOLCS/USGS catalog for the date before ruling seismic in or out.") + return 0 + + +if __name__ == "__main__": + sys.exit(main())