diff --git a/config/bernese/otl/README.md b/config/bernese/otl/README.md new file mode 100644 index 0000000..58f760e --- /dev/null +++ b/config/bernese/otl/README.md @@ -0,0 +1,104 @@ +# Ocean-loading coefficients for the seven IGS fiducials + +`LUZON.BLQ` covers the 135 local PHIVOLCS stations and **none of the seven IGS +fiducials** the network uses for 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 service, +which delivers **by email** — so one manual step is unavoidable. Everything +either side of it is prepared here. + +## ⚠ Match the existing settings — do not use the service defaults + +The 135 stations already in `LUZON.BLQ` were computed with a specific +configuration, recorded in the file's own header: + +``` +$$ Ocean tide model: FES2004 +$$ CMC: NO (corr.tide centre of mass) +$$ Gutenberg-Bullen Greens function is used +``` + +**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. + +## The request + +Paste `REQUEST_fiducials.txt`: + +``` +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 + +| 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 | + +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 +``` + +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 + +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 new file mode 100644 index 0000000..98c54f4 --- /dev/null +++ b/config/bernese/otl/REQUEST_fiducials.txt @@ -0,0 +1,7 @@ +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 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/docs/bernese54_luzon_reprocessing_runbook.md b/docs/bernese54_luzon_reprocessing_runbook.md index 6f8beba..98b5898 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 @@ -580,6 +597,321 @@ 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. + +### 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. + +### 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 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 **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. + +**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 | +|---|---|---| +| 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 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 +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 @@ -658,3 +990,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. diff --git a/docs/gps3-sessions/SESSION_LOG_20260729_storage.md b/docs/gps3-sessions/SESSION_LOG_20260729_storage.md index 8b32f80..a61e111 100644 --- a/docs/gps3-sessions/SESSION_LOG_20260729_storage.md +++ b/docs/gps3-sessions/SESSION_LOG_20260729_storage.md @@ -1381,3 +1381,255 @@ 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. + +--- + +## 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. + +### 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. 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()) diff --git a/scripts/derive_luzon_pcf.py b/scripts/derive_luzon_pcf.py index 070bd15..d7ca38a 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,27 @@ # 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". + # + # 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", @@ -116,12 +138,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") diff --git a/scripts/merge_blq.py b/scripts/merge_blq.py new file mode 100755 index 0000000..68efa2e --- /dev/null +++ b/scripts/merge_blq.py @@ -0,0 +1,252 @@ +#!/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 +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 +""" +from __future__ import annotations + +import argparse +import re +import shutil +import sys +from datetime import datetime +from pathlib import Path + +# 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 + # 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() + + +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: + name = _station_name(ln) + if name: + if current: + blocks[current] = buf + current = name + 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)") + + # 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 []): + rows = [ln for ln in incoming[name] + if ln.strip() and not ln.strip().startswith("$$") + 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)") + 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) + + # 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: + # 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, 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(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}") + return 0 + + +if __name__ == "__main__": + sys.exit(main()) 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()) diff --git a/scripts/run_luzon_month.sh b/scripts/run_luzon_month.sh new file mode 100755 index 0000000..f497998 --- /dev/null +++ b/scripts/run_luzon_month.sh @@ -0,0 +1,176 @@ +#!/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 + +# 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. +# +# 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. +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/$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 "$LUZON_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: $LUZON_YEAR DOY $LUZON_DOY_FROM-$LUZON_DOY_TO ===" + echo "started : $(date '+%F %T')" + 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 "$LUZON_SUMMARY" + +ok=0; failed=0; skipped=0; excluded=0; failed_days="" +start_all=$(date +%s) + +for doy in $(seq "$LUZON_DOY_FROM" "$LUZON_DOY_TO"); do + d3=$(printf '%03d' "$doy") + sess="${d3}0" + log="$LUZON_LOG_DIR/luzon-$LUZON_YEAR$sess.log" + + case " $LUZON_SKIP_DOYS " in + *" $doy "*) + printf ' DOY %s EXCLUDED (insufficient observations — see SKIP_DOYS)\n' \ + "$d3" | tee -a "$LUZON_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/$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 "$LUZON_SUMMARY" + + echo "start: $(date '+%F %T')" > "$log" + 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. + # + # 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+" "$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 "$LUZON_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: $LUZON_SKIP_DOYS (data gaps, not failures — see script header)" + fi + echo + 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 "$LUZON_SUMMARY" + +exit "$failed"