diff --git a/docs/FLOOD_FORECASTING.md b/docs/FLOOD_FORECASTING.md index 5622377..04b642e 100644 --- a/docs/FLOOD_FORECASTING.md +++ b/docs/FLOOD_FORECASTING.md @@ -536,6 +536,54 @@ exists alongside its flood events. (`fill_from_hii`, +9,341 h at P.81, +682 h at P.92, +810 h at P.20) is lead-neutral — the gate holds at 13 h with fill on — and ships enabled. +### 2026-09-12: three candidates on top of hgb-v3 — two rejected, one deferred + +Same rolling-origin harness (`src/ml/evaluate.py`, five monsoon folds +2021–2025, P.1 and P.103), all variants run from the identical +`models/cache/` snapshot (`--from-cache`), results in +`models/eval_2026-09-12*.json`, tables via `scripts/summarize_eval.py`. +Baseline is `rise_rain`, the deployed configuration. + +**Quantile regression heads (`rise_rain_quantile`, `_uw`) — rejected.** The +August result that quantile loss beat L2 on MAE held with rain in the model +(P.1 0.083/0.081 vs 0.087; P.103 0.143/0.124 vs 0.152), and Brier improved a +hair, but the operational numbers went the wrong way: at P.103 the 2022-08-14 +crossing dropped from +6 h to +1 h lead, 2022-10-02 from +9 h to +5/+3 h, and +the 2024-09-30 event from +9 h to +4 h; at P.1 2022 dropped +5 → +3/+2 h and +2025 +2 → +1 h, with one false-alarm episode where the baseline had none. A +median predicts the *typical* rise, and on the run-up to a crossing the typical +rise is not the one that matters. MAE is not the objective; lead is. + +**Quantile heads for sigma only (`rise_rain_qsigma`) — no effect.** The +hybrid keeps the L2 point prediction (so every lead is identical to the +baseline by construction — p≥0.5 alerts are sigma-independent) and derives a +per-row sigma from q90−q50. Brier moved 0.0031 → 0.0029 at P.1 and +0.0061 → 0.0060 at P.103, i.e. within noise, at the cost of three fitted +heads per horizon instead of one. Per-row uncertainty from this family of +models is not informative enough here to be worth the training time; the +0.15 m floor stays. + +**Forward-48 h forecast rain (`rise_rain_fc48`) — deferred.** Adding the +`(t, t+48]` Open-Meteo sum alongside `rain_fc24` left MAE, Brier and false +alarms unchanged and every event lead within ±1 h of baseline, *except* the +2024-10-03 P.1 record crossing, which went from +21 h to +72 h (and +55 → +69 h +at P.103). That is one event with the highest stakes in the record, on the +same feature family that already produced the 2024 gain, but n=1 is not +evidence: the P.103 2025-09-26 event lost 2 h in the same run. Rerun after the +2026 season adds events; if the 48 h window still moves only the biggest +onsets, promote it. Serving would need no new data source (`fetch_forecast` +already pulls `forecast_days=2`). + +**HII gauge rain — not evaluable yet.** `hii_rainfall` (~130 gauges in the +upper-Ping box, DWR/FOP/HII/RID/TMD) is the obvious independent rain source, +but the table only exists since 2026-08-11 and the api-v3 archive endpoint +ignores its date range (see `docs/DATA_SOURCES.md` §2.1), so every training +row before that is NaN and no fold in the harness has gauge data in its test +span. `src/ml/hii_rain.py` builds the catchment mean and +`GET /api/hii/rainfall/catchment` exposes it next to the Open-Meteo series with +a 24 h-sum bias/MAE/correlation, so the two sources' relationship is on record +by the time the 2027 fold (train ≤ 2027-04-30, test Jun–Nov 2027) can test it. + ## 6. Deployment ### API diff --git a/models/eval_2026-09-12.json b/models/eval_2026-09-12.json new file mode 100644 index 0000000..0dc7d68 --- /dev/null +++ b/models/eval_2026-09-12.json @@ -0,0 +1,723 @@ +[ + { + "station": "P.1", + "warn_thr": 3.7, + "folds": [ + { + "year": 2021, + "n_train": 20024, + "n_test": 4392, + "events": [], + "variants": { + "rise_rain": { + "mae": 0.07375926701460789, + "mae_above_2p5": null, + "brier_warn": 0.0, + "events": [], + "false_alarm_episodes": 0 + }, + "rise_rain_fc48": { + "mae": 0.07332598842242062, + "mae_above_2p5": null, + "brier_warn": 0.0, + "events": [], + "false_alarm_episodes": 0 + }, + "rise_rain_quantile": { + "mae": 0.07383022350644125, + "mae_above_2p5": null, + "brier_warn": 1.0850721383440065e-12, + "events": [], + "false_alarm_episodes": 0 + }, + "rise_rain_quantile_uw": { + "mae": 0.06977345537342049, + "mae_above_2p5": null, + "brier_warn": 7.128994064266462e-16, + "events": [], + "false_alarm_episodes": 0 + } + } + }, + { + "year": 2022, + "n_train": 28784, + "n_test": 4392, + "events": [ + { + "crossing": "2022-10-02T19:00:00", + "peak_ts": "2022-10-03T15:00:00", + "peak_level": 4.65 + } + ], + "variants": { + "rise_rain": { + "mae": 0.08452847354970178, + "mae_above_2p5": 0.25082252888260664, + "brier_warn": 0.004300908725927739, + "events": [ + { + "crossing": "2022-10-02T19:00:00", + "lead_h": 5.0, + "peak_level": 4.65, + "peak_pred_24h_before": 3.8173954245046406 + } + ], + "false_alarm_episodes": 0 + }, + "rise_rain_fc48": { + "mae": 0.08369773912557228, + "mae_above_2p5": 0.24702357745371217, + "brier_warn": 0.004325741658698973, + "events": [ + { + "crossing": "2022-10-02T19:00:00", + "lead_h": 5.0, + "peak_level": 4.65, + "peak_pred_24h_before": 3.8114190118860223 + } + ], + "false_alarm_episodes": 0 + }, + "rise_rain_quantile": { + "mae": 0.08230165237819607, + "mae_above_2p5": 0.2599241552494541, + "brier_warn": 0.004191550822623699, + "events": [ + { + "crossing": "2022-10-02T19:00:00", + "lead_h": 3.0, + "peak_level": 4.65, + "peak_pred_24h_before": 3.693734826616603 + } + ], + "false_alarm_episodes": 0 + }, + "rise_rain_quantile_uw": { + "mae": 0.08181109208140971, + "mae_above_2p5": 0.2639358812546371, + "brier_warn": 0.004581181754314636, + "events": [ + { + "crossing": "2022-10-02T19:00:00", + "lead_h": 2.0, + "peak_level": 4.65, + "peak_pred_24h_before": 3.620470606696475 + } + ], + "false_alarm_episodes": 0 + } + } + }, + { + "year": 2023, + "n_train": 37539, + "n_test": 4392, + "events": [], + "variants": { + "rise_rain": { + "mae": 0.07696637416933236, + "mae_above_2p5": 0.10643655855133666, + "brier_warn": 1.919860722404625e-15, + "events": [], + "false_alarm_episodes": 0 + }, + "rise_rain_fc48": { + "mae": 0.07708101918232377, + "mae_above_2p5": 0.09974894450126857, + "brier_warn": 1.7028444765622288e-15, + "events": [], + "false_alarm_episodes": 0 + }, + "rise_rain_quantile": { + "mae": 0.06968496247786034, + "mae_above_2p5": 0.10210094973854984, + "brier_warn": 1.0987601008721297e-09, + "events": [], + "false_alarm_episodes": 0 + }, + "rise_rain_quantile_uw": { + "mae": 0.06726603345661021, + "mae_above_2p5": 0.09716644572204487, + "brier_warn": 2.7085590803510675e-09, + "events": [], + "false_alarm_episodes": 0 + } + } + }, + { + "year": 2024, + "n_train": 46323, + "n_test": 4392, + "events": [ + { + "crossing": "2024-09-24T17:00:00", + "peak_ts": "2024-09-26T02:00:00", + "peak_level": 4.93 + }, + { + "crossing": "2024-10-03T09:00:00", + "peak_ts": "2024-10-05T12:00:00", + "peak_level": 5.3 + } + ], + "variants": { + "rise_rain": { + "mae": 0.0890019157543212, + "mae_above_2p5": 0.2093245927883516, + "brier_warn": 0.005826169840715695, + "events": [ + { + "crossing": "2024-09-24T17:00:00", + "lead_h": 10.0, + "peak_level": 4.93, + "peak_pred_24h_before": 4.817519939833057 + }, + { + "crossing": "2024-10-03T09:00:00", + "lead_h": 21.0, + "peak_level": 5.3, + "peak_pred_24h_before": 5.546588884631041 + } + ], + "false_alarm_episodes": 0 + }, + "rise_rain_fc48": { + "mae": 0.08875916089292855, + "mae_above_2p5": 0.20975114218621593, + "brier_warn": 0.006061154593275248, + "events": [ + { + "crossing": "2024-09-24T17:00:00", + "lead_h": 11.0, + "peak_level": 4.93, + "peak_pred_24h_before": 4.800924141216692 + }, + { + "crossing": "2024-10-03T09:00:00", + "lead_h": 72.0, + "peak_level": 5.3, + "peak_pred_24h_before": 5.577164984770105 + } + ], + "false_alarm_episodes": 0 + }, + "rise_rain_quantile": { + "mae": 0.08341835741803395, + "mae_above_2p5": 0.18688839374651373, + "brier_warn": 0.003188229521816099, + "events": [ + { + "crossing": "2024-09-24T17:00:00", + "lead_h": 17.0, + "peak_level": 4.93, + "peak_pred_24h_before": 5.058029430632501 + }, + { + "crossing": "2024-10-03T09:00:00", + "lead_h": 21.0, + "peak_level": 5.3, + "peak_pred_24h_before": 5.3311787370709975 + } + ], + "false_alarm_episodes": 0 + }, + "rise_rain_quantile_uw": { + "mae": 0.08412336465267245, + "mae_above_2p5": 0.19173218504592401, + "brier_warn": 0.0034722638888286116, + "events": [ + { + "crossing": "2024-09-24T17:00:00", + "lead_h": 15.0, + "peak_level": 4.93, + "peak_pred_24h_before": 5.0118284217314 + }, + { + "crossing": "2024-10-03T09:00:00", + "lead_h": 21.0, + "peak_level": 5.3, + "peak_pred_24h_before": 5.3520431553190155 + } + ], + "false_alarm_episodes": 0 + } + } + }, + { + "year": 2025, + "n_train": 55083, + "n_test": 4392, + "events": [ + { + "crossing": "2025-09-27T18:00:00", + "peak_ts": "2025-09-27T22:00:00", + "peak_level": 3.93 + } + ], + "variants": { + "rise_rain": { + "mae": 0.11202835461695825, + "mae_above_2p5": 0.2456782496767171, + "brier_warn": 0.005178356039745929, + "events": [ + { + "crossing": "2025-09-27T18:00:00", + "lead_h": 2.0, + "peak_level": 3.93, + "peak_pred_24h_before": 3.23 + } + ], + "false_alarm_episodes": 0 + }, + "rise_rain_fc48": { + "mae": 0.11083068059655521, + "mae_above_2p5": 0.23787013498709667, + "brier_warn": 0.005228043825289021, + "events": [ + { + "crossing": "2025-09-27T18:00:00", + "lead_h": 2.0, + "peak_level": 3.93, + "peak_pred_24h_before": 3.23 + } + ], + "false_alarm_episodes": 0 + }, + "rise_rain_quantile": { + "mae": 0.10372526043263834, + "mae_above_2p5": 0.23846363852091013, + "brier_warn": 0.005898259026876986, + "events": [ + { + "crossing": "2025-09-27T18:00:00", + "lead_h": 1.0, + "peak_level": 3.93, + "peak_pred_24h_before": 3.23 + } + ], + "false_alarm_episodes": 1 + }, + "rise_rain_quantile_uw": { + "mae": 0.10267267834487567, + "mae_above_2p5": 0.2507951319537275, + "brier_warn": 0.005020467408392599, + "events": [ + { + "crossing": "2025-09-27T18:00:00", + "lead_h": 2.0, + "peak_level": 3.93, + "peak_pred_24h_before": 3.2417205711942434 + } + ], + "false_alarm_episodes": 1 + } + } + } + ] + }, + { + "station": "P.103", + "warn_thr": 5.95, + "folds": [ + { + "year": 2021, + "n_train": 20009, + "n_test": 4392, + "events": [], + "variants": { + "rise_rain": { + "mae": 0.14438683486790602, + "mae_above_2p5": null, + "brier_warn": 0.0, + "events": [], + "false_alarm_episodes": 0 + }, + "rise_rain_fc48": { + "mae": 0.14305613024043953, + "mae_above_2p5": null, + "brier_warn": 0.0, + "events": [], + "false_alarm_episodes": 0 + }, + "rise_rain_quantile": { + "mae": 0.13512418754716052, + "mae_above_2p5": null, + "brier_warn": 6.937198832251013e-10, + "events": [], + "false_alarm_episodes": 0 + }, + "rise_rain_quantile_uw": { + "mae": 0.10270057685061172, + "mae_above_2p5": null, + "brier_warn": 6.565945930999852e-14, + "events": [], + "false_alarm_episodes": 0 + } + } + }, + { + "year": 2022, + "n_train": 28769, + "n_test": 4392, + "events": [ + { + "crossing": "2022-08-14T04:00:00", + "peak_ts": "2022-08-14T08:00:00", + "peak_level": 6.09 + }, + { + "crossing": "2022-10-02T16:00:00", + "peak_ts": "2022-10-03T16:00:00", + "peak_level": 7.54 + } + ], + "variants": { + "rise_rain": { + "mae": 0.14392334606606189, + "mae_above_2p5": 0.41075382386412596, + "brier_warn": 0.00764679988213082, + "events": [ + { + "crossing": "2022-08-14T04:00:00", + "lead_h": 6.0, + "peak_level": 6.09, + "peak_pred_24h_before": 4.825270553204425 + }, + { + "crossing": "2022-10-02T16:00:00", + "lead_h": 9.0, + "peak_level": 7.54, + "peak_pred_24h_before": 7.194064117976157 + } + ], + "false_alarm_episodes": 0 + }, + "rise_rain_fc48": { + "mae": 0.14284722696431765, + "mae_above_2p5": 0.39401132538096667, + "brier_warn": 0.007427375799450148, + "events": [ + { + "crossing": "2022-08-14T04:00:00", + "lead_h": 6.0, + "peak_level": 6.09, + "peak_pred_24h_before": 4.845896993132048 + }, + { + "crossing": "2022-10-02T16:00:00", + "lead_h": 9.0, + "peak_level": 7.54, + "peak_pred_24h_before": 7.194661123502819 + } + ], + "false_alarm_episodes": 0 + }, + "rise_rain_quantile": { + "mae": 0.1458954690776336, + "mae_above_2p5": 0.4674051188369517, + "brier_warn": 0.00899991317044724, + "events": [ + { + "crossing": "2022-08-14T04:00:00", + "lead_h": 1.0, + "peak_level": 6.09, + "peak_pred_24h_before": 4.82714042795344 + }, + { + "crossing": "2022-10-02T16:00:00", + "lead_h": 5.0, + "peak_level": 7.54, + "peak_pred_24h_before": 6.59871451008115 + } + ], + "false_alarm_episodes": 0 + }, + "rise_rain_quantile_uw": { + "mae": 0.13203958422483053, + "mae_above_2p5": 0.4902773906613983, + "brier_warn": 0.008514527559601331, + "events": [ + { + "crossing": "2022-08-14T04:00:00", + "lead_h": 4.0, + "peak_level": 6.09, + "peak_pred_24h_before": 4.895369771408423 + }, + { + "crossing": "2022-10-02T16:00:00", + "lead_h": 3.0, + "peak_level": 7.54, + "peak_pred_24h_before": 6.294082735853337 + } + ], + "false_alarm_episodes": 0 + } + } + }, + { + "year": 2023, + "n_train": 37524, + "n_test": 4392, + "events": [], + "variants": { + "rise_rain": { + "mae": 0.16195301512514512, + "mae_above_2p5": null, + "brier_warn": 3.0208441147560167e-07, + "events": [], + "false_alarm_episodes": 0 + }, + "rise_rain_fc48": { + "mae": 0.16334224973870387, + "mae_above_2p5": null, + "brier_warn": 1.1626187140381066e-08, + "events": [], + "false_alarm_episodes": 0 + }, + "rise_rain_quantile": { + "mae": 0.14755024879522577, + "mae_above_2p5": null, + "brier_warn": 3.850104714388072e-05, + "events": [], + "false_alarm_episodes": 0 + }, + "rise_rain_quantile_uw": { + "mae": 0.10589805327026759, + "mae_above_2p5": null, + "brier_warn": 1.8896757611081001e-07, + "events": [], + "false_alarm_episodes": 0 + } + } + }, + { + "year": 2024, + "n_train": 46308, + "n_test": 4058, + "events": [ + { + "crossing": "2024-09-24T10:00:00", + "peak_ts": "2024-09-26T00:00:00", + "peak_level": 8.27 + }, + { + "crossing": "2024-09-30T03:00:00", + "peak_ts": "2024-09-30T06:00:00", + "peak_level": 5.99 + }, + { + "crossing": "2024-10-03T06:00:00", + "peak_ts": "2024-10-05T07:00:00", + "peak_level": 9.93 + } + ], + "variants": { + "rise_rain": { + "mae": 0.13901842325128816, + "mae_above_2p5": 0.2887416310966684, + "brier_warn": 0.0072040857753137046, + "events": [ + { + "crossing": "2024-09-24T10:00:00", + "lead_h": 19.0, + "peak_level": 8.27, + "peak_pred_24h_before": 8.212734363860193 + }, + { + "crossing": "2024-09-30T03:00:00", + "lead_h": 9.0, + "peak_level": 5.99, + "peak_pred_24h_before": 5.47565489914091 + }, + { + "crossing": "2024-10-03T06:00:00", + "lead_h": 55.0, + "peak_level": 9.93, + "peak_pred_24h_before": 8.780278464915938 + } + ], + "false_alarm_episodes": 0 + }, + "rise_rain_fc48": { + "mae": 0.13850733053801131, + "mae_above_2p5": 0.29817392926900865, + "brier_warn": 0.007928846574280278, + "events": [ + { + "crossing": "2024-09-24T10:00:00", + "lead_h": 19.0, + "peak_level": 8.27, + "peak_pred_24h_before": 8.166020251020889 + }, + { + "crossing": "2024-09-30T03:00:00", + "lead_h": 10.0, + "peak_level": 5.99, + "peak_pred_24h_before": 5.524584037259288 + }, + { + "crossing": "2024-10-03T06:00:00", + "lead_h": 69.0, + "peak_level": 9.93, + "peak_pred_24h_before": 8.595120917488185 + } + ], + "false_alarm_episodes": 0 + }, + "rise_rain_quantile": { + "mae": 0.12640546558686208, + "mae_above_2p5": 0.27944115421093924, + "brier_warn": 0.007730268539879654, + "events": [ + { + "crossing": "2024-09-24T10:00:00", + "lead_h": 20.0, + "peak_level": 8.27, + "peak_pred_24h_before": 8.336487732683672 + }, + { + "crossing": "2024-09-30T03:00:00", + "lead_h": 4.0, + "peak_level": 5.99, + "peak_pred_24h_before": 5.45284915024346 + }, + { + "crossing": "2024-10-03T06:00:00", + "lead_h": 69.0, + "peak_level": 9.93, + "peak_pred_24h_before": 8.575915226697406 + } + ], + "false_alarm_episodes": 0 + }, + "rise_rain_quantile_uw": { + "mae": 0.12480975852168202, + "mae_above_2p5": 0.2857575557415695, + "brier_warn": 0.005894928998647136, + "events": [ + { + "crossing": "2024-09-24T10:00:00", + "lead_h": 21.0, + "peak_level": 8.27, + "peak_pred_24h_before": 8.444720326882825 + }, + { + "crossing": "2024-09-30T03:00:00", + "lead_h": 9.0, + "peak_level": 5.99, + "peak_pred_24h_before": 5.42 + }, + { + "crossing": "2024-10-03T06:00:00", + "lead_h": 32.0, + "peak_level": 9.93, + "peak_pred_24h_before": 8.781001847167515 + } + ], + "false_alarm_episodes": 0 + } + } + }, + { + "year": 2025, + "n_train": 54577, + "n_test": 4392, + "events": [ + { + "crossing": "2025-09-26T06:00:00", + "peak_ts": "2025-09-27T21:00:00", + "peak_level": 6.64 + }, + { + "crossing": "2025-10-03T06:00:00", + "peak_ts": "2025-10-03T12:00:00", + "peak_level": 6.14 + } + ], + "variants": { + "rise_rain": { + "mae": 0.17206685418009837, + "mae_above_2p5": 0.3840520085986099, + "brier_warn": 0.015824852891785323, + "events": [ + { + "crossing": "2025-09-26T06:00:00", + "lead_h": 6.0, + "peak_level": 6.64, + "peak_pred_24h_before": 5.73043665401637 + }, + { + "crossing": "2025-10-03T06:00:00", + "lead_h": 8.0, + "peak_level": 6.14, + "peak_pred_24h_before": 5.609919760117381 + } + ], + "false_alarm_episodes": 2 + }, + "rise_rain_fc48": { + "mae": 0.1722274451362234, + "mae_above_2p5": 0.3872216974630654, + "brier_warn": 0.016245727130063177, + "events": [ + { + "crossing": "2025-09-26T06:00:00", + "lead_h": 4.0, + "peak_level": 6.64, + "peak_pred_24h_before": 5.7312263707793685 + }, + { + "crossing": "2025-10-03T06:00:00", + "lead_h": 8.0, + "peak_level": 6.14, + "peak_pred_24h_before": 5.685511824758637 + } + ], + "false_alarm_episodes": 2 + }, + "rise_rain_quantile": { + "mae": 0.15912275848686555, + "mae_above_2p5": 0.3845695612674703, + "brier_warn": 0.014797606329811499, + "events": [ + { + "crossing": "2025-09-26T06:00:00", + "lead_h": 11.0, + "peak_level": 6.64, + "peak_pred_24h_before": 5.73 + }, + { + "crossing": "2025-10-03T06:00:00", + "lead_h": 10.0, + "peak_level": 6.14, + "peak_pred_24h_before": 5.58133012298031 + } + ], + "false_alarm_episodes": 2 + }, + "rise_rain_quantile_uw": { + "mae": 0.15385355845103685, + "mae_above_2p5": 0.41314645854725585, + "brier_warn": 0.016528116336109958, + "events": [ + { + "crossing": "2025-09-26T06:00:00", + "lead_h": 9.0, + "peak_level": 6.64, + "peak_pred_24h_before": 5.741254234340573 + }, + { + "crossing": "2025-10-03T06:00:00", + "lead_h": 8.0, + "peak_level": 6.14, + "peak_pred_24h_before": 5.418105405306207 + } + ], + "false_alarm_episodes": 2 + } + } + } + ] + } +] \ No newline at end of file diff --git a/models/eval_2026-09-12_qsigma.json b/models/eval_2026-09-12_qsigma.json new file mode 100644 index 0000000..58e7f00 --- /dev/null +++ b/models/eval_2026-09-12_qsigma.json @@ -0,0 +1,439 @@ +[ + { + "station": "P.1", + "warn_thr": 3.7, + "folds": [ + { + "year": 2021, + "n_train": 20024, + "n_test": 4392, + "events": [], + "variants": { + "rise_rain": { + "mae": 0.07375926701460789, + "mae_above_2p5": null, + "brier_warn": 0.0, + "events": [], + "false_alarm_episodes": 0 + }, + "rise_rain_qsigma": { + "mae": 0.07375926701460789, + "mae_above_2p5": null, + "brier_warn": 7.852802408474157e-14, + "events": [], + "false_alarm_episodes": 0 + } + } + }, + { + "year": 2022, + "n_train": 28784, + "n_test": 4392, + "events": [ + { + "crossing": "2022-10-02T19:00:00", + "peak_ts": "2022-10-03T15:00:00", + "peak_level": 4.65 + } + ], + "variants": { + "rise_rain": { + "mae": 0.08452847354970178, + "mae_above_2p5": 0.25082252888260664, + "brier_warn": 0.004300908725927739, + "events": [ + { + "crossing": "2022-10-02T19:00:00", + "lead_h": 5.0, + "peak_level": 4.65, + "peak_pred_24h_before": 3.8173954245046406 + } + ], + "false_alarm_episodes": 0 + }, + "rise_rain_qsigma": { + "mae": 0.08452847354970178, + "mae_above_2p5": 0.25082252888260664, + "brier_warn": 0.0039815091186836665, + "events": [ + { + "crossing": "2022-10-02T19:00:00", + "lead_h": 5.0, + "peak_level": 4.65, + "peak_pred_24h_before": 3.8173954245046406 + } + ], + "false_alarm_episodes": 0 + } + } + }, + { + "year": 2023, + "n_train": 37539, + "n_test": 4392, + "events": [], + "variants": { + "rise_rain": { + "mae": 0.07696637416933236, + "mae_above_2p5": 0.10643655855133666, + "brier_warn": 1.919860722404625e-15, + "events": [], + "false_alarm_episodes": 0 + }, + "rise_rain_qsigma": { + "mae": 0.07696637416933236, + "mae_above_2p5": 0.10643655855133666, + "brier_warn": 4.0844243581898366e-08, + "events": [], + "false_alarm_episodes": 0 + } + } + }, + { + "year": 2024, + "n_train": 46323, + "n_test": 4392, + "events": [ + { + "crossing": "2024-09-24T17:00:00", + "peak_ts": "2024-09-26T02:00:00", + "peak_level": 4.93 + }, + { + "crossing": "2024-10-03T09:00:00", + "peak_ts": "2024-10-05T12:00:00", + "peak_level": 5.3 + } + ], + "variants": { + "rise_rain": { + "mae": 0.0890019157543212, + "mae_above_2p5": 0.2093245927883516, + "brier_warn": 0.005826169840715695, + "events": [ + { + "crossing": "2024-09-24T17:00:00", + "lead_h": 10.0, + "peak_level": 4.93, + "peak_pred_24h_before": 4.817519939833057 + }, + { + "crossing": "2024-10-03T09:00:00", + "lead_h": 21.0, + "peak_level": 5.3, + "peak_pred_24h_before": 5.546588884631041 + } + ], + "false_alarm_episodes": 0 + }, + "rise_rain_qsigma": { + "mae": 0.0890019157543212, + "mae_above_2p5": 0.2093245927883516, + "brier_warn": 0.005475787026695418, + "events": [ + { + "crossing": "2024-09-24T17:00:00", + "lead_h": 10.0, + "peak_level": 4.93, + "peak_pred_24h_before": 4.817519939833057 + }, + { + "crossing": "2024-10-03T09:00:00", + "lead_h": 21.0, + "peak_level": 5.3, + "peak_pred_24h_before": 5.546588884631041 + } + ], + "false_alarm_episodes": 0 + } + } + }, + { + "year": 2025, + "n_train": 55083, + "n_test": 4392, + "events": [ + { + "crossing": "2025-09-27T18:00:00", + "peak_ts": "2025-09-27T22:00:00", + "peak_level": 3.93 + } + ], + "variants": { + "rise_rain": { + "mae": 0.11202835461695825, + "mae_above_2p5": 0.2456782496767171, + "brier_warn": 0.005178356039745929, + "events": [ + { + "crossing": "2025-09-27T18:00:00", + "lead_h": 2.0, + "peak_level": 3.93, + "peak_pred_24h_before": 3.23 + } + ], + "false_alarm_episodes": 0 + }, + "rise_rain_qsigma": { + "mae": 0.11202835461695825, + "mae_above_2p5": 0.2456782496767171, + "brier_warn": 0.00491346243876589, + "events": [ + { + "crossing": "2025-09-27T18:00:00", + "lead_h": 2.0, + "peak_level": 3.93, + "peak_pred_24h_before": 3.23 + } + ], + "false_alarm_episodes": 0 + } + } + } + ] + }, + { + "station": "P.103", + "warn_thr": 5.95, + "folds": [ + { + "year": 2021, + "n_train": 20009, + "n_test": 4392, + "events": [], + "variants": { + "rise_rain": { + "mae": 0.14438683486790602, + "mae_above_2p5": null, + "brier_warn": 0.0, + "events": [], + "false_alarm_episodes": 0 + }, + "rise_rain_qsigma": { + "mae": 0.14438683486790602, + "mae_above_2p5": null, + "brier_warn": 3.5349670949872053e-13, + "events": [], + "false_alarm_episodes": 0 + } + } + }, + { + "year": 2022, + "n_train": 28769, + "n_test": 4392, + "events": [ + { + "crossing": "2022-08-14T04:00:00", + "peak_ts": "2022-08-14T08:00:00", + "peak_level": 6.09 + }, + { + "crossing": "2022-10-02T16:00:00", + "peak_ts": "2022-10-03T16:00:00", + "peak_level": 7.54 + } + ], + "variants": { + "rise_rain": { + "mae": 0.14392334606606189, + "mae_above_2p5": 0.41075382386412596, + "brier_warn": 0.00764679988213082, + "events": [ + { + "crossing": "2022-08-14T04:00:00", + "lead_h": 6.0, + "peak_level": 6.09, + "peak_pred_24h_before": 4.825270553204425 + }, + { + "crossing": "2022-10-02T16:00:00", + "lead_h": 9.0, + "peak_level": 7.54, + "peak_pred_24h_before": 7.194064117976157 + } + ], + "false_alarm_episodes": 0 + }, + "rise_rain_qsigma": { + "mae": 0.14392334606606189, + "mae_above_2p5": 0.41075382386412596, + "brier_warn": 0.006975493312169236, + "events": [ + { + "crossing": "2022-08-14T04:00:00", + "lead_h": 6.0, + "peak_level": 6.09, + "peak_pred_24h_before": 4.825270553204425 + }, + { + "crossing": "2022-10-02T16:00:00", + "lead_h": 9.0, + "peak_level": 7.54, + "peak_pred_24h_before": 7.194064117976157 + } + ], + "false_alarm_episodes": 0 + } + } + }, + { + "year": 2023, + "n_train": 37524, + "n_test": 4392, + "events": [], + "variants": { + "rise_rain": { + "mae": 0.16195301512514512, + "mae_above_2p5": null, + "brier_warn": 3.0208441147560167e-07, + "events": [], + "false_alarm_episodes": 0 + }, + "rise_rain_qsigma": { + "mae": 0.16195301512514512, + "mae_above_2p5": null, + "brier_warn": 3.2176039819636275e-05, + "events": [], + "false_alarm_episodes": 0 + } + } + }, + { + "year": 2024, + "n_train": 46308, + "n_test": 4058, + "events": [ + { + "crossing": "2024-09-24T10:00:00", + "peak_ts": "2024-09-26T00:00:00", + "peak_level": 8.27 + }, + { + "crossing": "2024-09-30T03:00:00", + "peak_ts": "2024-09-30T06:00:00", + "peak_level": 5.99 + }, + { + "crossing": "2024-10-03T06:00:00", + "peak_ts": "2024-10-05T07:00:00", + "peak_level": 9.93 + } + ], + "variants": { + "rise_rain": { + "mae": 0.13901842325128816, + "mae_above_2p5": 0.2887416310966684, + "brier_warn": 0.0072040857753137046, + "events": [ + { + "crossing": "2024-09-24T10:00:00", + "lead_h": 19.0, + "peak_level": 8.27, + "peak_pred_24h_before": 8.212734363860193 + }, + { + "crossing": "2024-09-30T03:00:00", + "lead_h": 9.0, + "peak_level": 5.99, + "peak_pred_24h_before": 5.47565489914091 + }, + { + "crossing": "2024-10-03T06:00:00", + "lead_h": 55.0, + "peak_level": 9.93, + "peak_pred_24h_before": 8.780278464915938 + } + ], + "false_alarm_episodes": 0 + }, + "rise_rain_qsigma": { + "mae": 0.13901842325128816, + "mae_above_2p5": 0.2887416310966684, + "brier_warn": 0.006910847607098417, + "events": [ + { + "crossing": "2024-09-24T10:00:00", + "lead_h": 19.0, + "peak_level": 8.27, + "peak_pred_24h_before": 8.212734363860193 + }, + { + "crossing": "2024-09-30T03:00:00", + "lead_h": 9.0, + "peak_level": 5.99, + "peak_pred_24h_before": 5.47565489914091 + }, + { + "crossing": "2024-10-03T06:00:00", + "lead_h": 55.0, + "peak_level": 9.93, + "peak_pred_24h_before": 8.780278464915938 + } + ], + "false_alarm_episodes": 0 + } + } + }, + { + "year": 2025, + "n_train": 54577, + "n_test": 4392, + "events": [ + { + "crossing": "2025-09-26T06:00:00", + "peak_ts": "2025-09-27T21:00:00", + "peak_level": 6.64 + }, + { + "crossing": "2025-10-03T06:00:00", + "peak_ts": "2025-10-03T12:00:00", + "peak_level": 6.14 + } + ], + "variants": { + "rise_rain": { + "mae": 0.17206685418009837, + "mae_above_2p5": 0.3840520085986099, + "brier_warn": 0.015824852891785323, + "events": [ + { + "crossing": "2025-09-26T06:00:00", + "lead_h": 6.0, + "peak_level": 6.64, + "peak_pred_24h_before": 5.73043665401637 + }, + { + "crossing": "2025-10-03T06:00:00", + "lead_h": 8.0, + "peak_level": 6.14, + "peak_pred_24h_before": 5.609919760117381 + } + ], + "false_alarm_episodes": 2 + }, + "rise_rain_qsigma": { + "mae": 0.17206685418009837, + "mae_above_2p5": 0.3840520085986099, + "brier_warn": 0.015889933586185904, + "events": [ + { + "crossing": "2025-09-26T06:00:00", + "lead_h": 6.0, + "peak_level": 6.64, + "peak_pred_24h_before": 5.73043665401637 + }, + { + "crossing": "2025-10-03T06:00:00", + "lead_h": 8.0, + "peak_level": 6.14, + "peak_pred_24h_before": 5.609919760117381 + } + ], + "false_alarm_episodes": 2 + } + } + } + ] + } +] \ No newline at end of file diff --git a/scripts/summarize_eval.py b/scripts/summarize_eval.py new file mode 100644 index 0000000..995e398 --- /dev/null +++ b/scripts/summarize_eval.py @@ -0,0 +1,62 @@ +"""Summarise rolling-origin harness output side by side. + +Usage: + uv run python scripts/summarize_eval.py models/eval_2026-09-12.json [more.json ...] + +Aggregates each (station, variant) across folds: mean MAE, mean flood-regime +MAE, mean Brier, total false-alarm episodes, and every warning event with its +first-alert lead and the 24 h-ahead peak error -- the operational numbers that +decide whether a variant ships. +""" + +import json +import statistics +import sys +from collections import OrderedDict + + +def summarize(paths): + for path in paths: + results = json.load(open(path, encoding="utf-8")) + print(f"\n##### {path}") + for station in results: + print(f"\n=== {station['station']} (warn {station['warn_thr']:.2f} m) ===") + agg = OrderedDict() + for fold in station["folds"]: + for name, m in fold["variants"].items(): + a = agg.setdefault( + name, {"mae": [], "mae_hi": [], "brier": [], "fa": 0, "events": []} + ) + a["mae"].append(m["mae"]) + if m.get("mae_above_2p5") is not None: + a["mae_hi"].append(m["mae_above_2p5"]) + if m.get("brier_warn") is not None: + a["brier"].append(m["brier_warn"]) + a["fa"] += m["false_alarm_episodes"] + for e in m["events"]: + err = ( + None + if e["peak_pred_24h_before"] is None + else e["peak_pred_24h_before"] - e["peak_level"] + ) + a["events"].append((fold["year"], e["crossing"][:10], e["lead_h"], e["peak_level"], err)) + print(f"{'variant':22} {'MAE':>6} {'MAE_hi':>7} {'Brier':>7} {'FA':>3} events: year crossing lead_h peak(err24h)") + for name, a in agg.items(): + ev = " ".join( + f"{y} {d} {'—' if l is None else format(l, '+.0f')}h {p:.2f}({'—' if err is None else format(err, '+.2f')})" + for y, d, l, p, err in a["events"] + ) + leads = [l for *_, l, _, _ in a["events"] if l is not None] + print( + f"{name:22} {statistics.mean(a['mae']):6.3f} " + f"{statistics.mean(a['mae_hi']) if a['mae_hi'] else float('nan'):7.3f} " + f"{statistics.mean(a['brier']) if a['brier'] else float('nan'):7.4f} " + f"{a['fa']:>3} {ev}" + ) + if leads: + print(f"{'':22} lead: mean {statistics.mean(leads):+.1f} h, min {min(leads):+.0f} h, " + f"missed {sum(1 for *_, l, _, _ in a['events'] if l is None)}/{len(a['events'])}") + + +if __name__ == "__main__": + summarize(sys.argv[1:] or ["models/eval_variants.json"]) diff --git a/src/ml/data.py b/src/ml/data.py index 3699a12..5d43a16 100644 --- a/src/ml/data.py +++ b/src/ml/data.py @@ -282,6 +282,10 @@ def _read_cache(cache_dir: Path, stations: Optional[List[str]]) -> pd.DataFrame: frames = [] for path in sorted(cache_dir.glob("*.csv.gz")): code = path.name[: -len(".csv.gz")] + # The dir is shared with rain.py / dam.py caches (rain_openmeteo, + # dam_): only station files (P.) are measurements. + if not code.startswith("P."): + continue if stations and code not in stations: continue with gzip.open(path, "rt", encoding="utf-8") as handle: diff --git a/src/ml/evaluate.py b/src/ml/evaluate.py index 57e9e2a..1043faa 100644 --- a/src/ml/evaluate.py +++ b/src/ml/evaluate.py @@ -52,18 +52,31 @@ def _flood_weights(y_abs: pd.Series) -> np.ndarray: return 1.0 + 4.0 * np.clip((y_abs.to_numpy() - 2.5) / 1.2, 0.0, 1.0) +# Experimental forward-48h rain sum, built in evaluate_station (not in +# features.build_features) so the served feature set is untouched until the +# harness says it helps. Serving could supply it: fetch_forecast() already +# pulls forecast_days=2. +EXTRA_RAIN_FEATURES = ("rain_fc48",) + + class Variant: """A trainable candidate producing (pred_abs, sigma_per_row) on test rows.""" def __init__(self, name: str, target: str, weighted: bool = False, quantile: bool = False, use_rain: bool = False, - use_dam: bool = False): + use_dam: bool = False, use_fc48: bool = False, + qsigma: bool = False): self.name = name self.target = target # 'abs' or 'rise' self.weighted = weighted self.quantile = quantile self.use_rain = use_rain self.use_dam = use_dam + self.use_fc48 = use_fc48 + # Hybrid: L2 head for the point prediction (keeps the lead-time + # behaviour of the deployed model exactly, since p>=0.5 alerts are + # sigma-independent) and quantile heads ONLY for a per-row sigma. + self.qsigma = qsigma def fit_predict( self, X_tr, y_abs_tr, X_te @@ -84,6 +97,12 @@ class Variant: raise ValueError( f"{self.name} requires the dam series (rid_reservoir_daily backfilled)" ) + if not self.use_fc48: + drop = [c for c in EXTRA_RAIN_FEATURES if c in X_tr.columns] + X_tr = X_tr.drop(columns=drop) + X_te = X_te.drop(columns=drop) + elif "rain_fc48" not in X_tr.columns: + raise ValueError(f"{self.name} requires the rain series") level_tr = X_tr["level"] level_te = X_te["level"].to_numpy() y_tr = (y_abs_tr - level_tr) if self.target == "rise" else y_abs_tr @@ -99,7 +118,13 @@ class Variant: else: reg = _make_regressor().fit(X_tr, y_tr, sample_weight=weights) pred = reg.predict(X_te) - sigma = np.full(len(X_te), FIXED_SIGMA) + if self.qsigma: + q50 = _quantile_regressor(0.5).fit(X_tr, y_tr, sample_weight=weights) + q90 = _quantile_regressor(0.9).fit(X_tr, y_tr, sample_weight=weights) + spread = np.maximum(q90.predict(X_te) - q50.predict(X_te), 0.0) + sigma = np.maximum(spread / 1.2816, 0.05) + else: + sigma = np.full(len(X_te), FIXED_SIGMA) pred_abs = pred + level_te if self.target == "rise" else pred pred_abs = np.maximum(pred_abs, level_te) # peak >= current, as served @@ -116,12 +141,29 @@ VARIANTS: Dict[str, Variant] = { "rise_rain_dam": Variant("rise_rain_dam", target="rise", use_rain=True, use_dam=True), "rise_dam": Variant("rise_dam", target="rise", use_dam=True), + # 2026-09-12 experiments on top of the deployed rise_rain configuration: + # per-row sigma from quantile heads (the served sigma sits on the 0.15 + # floor at every P.1 horizon, so stage probabilities are constant- + # calibrated), and a longer forecast-rain window for the 24 h horizon. + "rise_rain_quantile": Variant("rise_rain_quantile", target="rise", + weighted=True, quantile=True, use_rain=True), + "rise_rain_quantile_uw": Variant("rise_rain_quantile_uw", target="rise", + quantile=True, use_rain=True), + "rise_rain_fc48": Variant("rise_rain_fc48", target="rise", use_rain=True, + use_fc48=True), + "rise_rain_qsigma": Variant("rise_rain_qsigma", target="rise", use_rain=True, + qsigma=True), } # Dam variants are opt-in by name: they require dam columns that only exist # for features.DAM_STATIONS and only when the reservoir series loaded, and -# the 2026-08-13 ablation concluded them a negative result. -DEFAULT_VARIANTS = [k for k, v in VARIANTS.items() if not v.use_dam] +# the 2026-08-13 ablation concluded them a negative result. The 2026-09-12 +# experiments are opt-in too (see their results in docs/FLOOD_FORECASTING.md). +DEFAULT_VARIANTS = [ + k for k, v in VARIANTS.items() + if not v.use_dam and not v.use_fc48 and not v.qsigma + and not (v.quantile and v.use_rain) +] def _find_events(observed: pd.Series, thr: float) -> List[dict]: @@ -222,6 +264,12 @@ def evaluate_station( warn_thr, _ = features.get_thresholds(station) grid = features.make_hourly_grid(df_long) X_all = features.build_features(grid, station, rain=rain, dam=dam) + if rain is not None: + # forward sum over (t, t+48]; same construction as rain_fc24 + r = rain.reindex(X_all.index) + X_all["rain_fc48"] = ( + r.shift(-1).iloc[::-1].rolling(48, min_periods=1).sum().iloc[::-1] + ) observed = grid.observed[(station, "water_level")] keep = X_all["obs_age_h"].notna() @@ -380,12 +428,19 @@ def main(argv=None) -> int: help="skip loading the Open-Meteo rain series") parser.add_argument("--no-dam", action="store_true", help="skip loading the Mae Ngat reservoir series") + parser.add_argument("--from-cache", action="store_true", + help="offline: read models/cache/ only (no DB, no API, " + "no Open-Meteo refresh) -- reproducible reruns") args = parser.parse_args(argv) logging.basicConfig( level=logging.INFO, format="%(asctime)s %(levelname)s %(message)s" ) - df = data.load_measurements(db_url=args.db_url) + if args.from_cache: + df = data._read_cache(data.CACHE_DIR, None) + logger.info(f"measurements from cache: {len(df)} rows") + else: + df = data.load_measurements(db_url=args.db_url) if df.empty: logger.error("no measurement data") return 1 @@ -394,7 +449,9 @@ def main(argv=None) -> int: if not args.no_rain: from . import rain as rain_mod - rain_series = rain_mod.catchment_mean(rain_mod.load_history()) + rain_series = rain_mod.catchment_mean( + rain_mod.load_history(refresh=not args.from_cache) + ) if rain_series is None: logger.warning("rain history unavailable; rain features will be NaN") else: @@ -404,7 +461,7 @@ def main(argv=None) -> int: ) dam_frame = None - if not args.no_dam: + if not args.no_dam and not args.from_cache: from . import dam as dam_mod dam_frame = dam_mod.load_history(db_url=args.db_url) diff --git a/src/ml/hii_rain.py b/src/ml/hii_rain.py new file mode 100644 index 0000000..92a81d9 --- /dev/null +++ b/src/ml/hii_rain.py @@ -0,0 +1,119 @@ +"""Catchment-mean hourly rain from the HII/ThaiWater gauge network. + +Independent of Open-Meteo (src/ml/rain.py): those are model-analysis values, +these are what the gauges measured. The `hii_rainfall` table has been filled +by the hourly collector since 2026-08-11 and there is NO archive behind it +(the api-v3 rain_24h_graph endpoint ignores its date range, see +docs/DATA_SOURCES.md 2.1), so this series cannot yet be a training feature: +every training row before 2026-08 would be NaN and HistGradientBoosting +would learn nothing from the column. It becomes a candidate once a full +monsoon season of gauge rows exists in the rolling-origin harness's test +span -- the 2027 fold (train through 2027-04-30, test Jun-Nov 2027) is the +first that could show anything. + +Until then it serves two purposes: + * a live cross-check of the Open-Meteo catchment mean (/api/hii/rainfall + already exposes the raw gauges; this gives the comparable aggregate); + * accumulating the comparison so the eventual feature evaluation has a + documented bias/variance relationship between the two sources. +""" + +import logging +from typing import Optional, Sequence, Tuple + +import pandas as pd + +from .data import resolve_db_url + +logger = logging.getLogger(__name__) + +# Same footprint as rain.CATCHMENT_POINTS: the upper Ping above P.1. Gauges +# inside this box are averaged; there are ~130 with recent data (DWR, FOP, +# HII, RID, TMD), far denser than the five Open-Meteo points. +CATCHMENT_BOX: Tuple[float, float, float, float] = (18.75, 19.60, 98.60, 99.30) +# A gauge that reports the same rain_24h for many hours is stuck; drop hours +# where fewer than this many gauges reported at all. +MIN_GAUGES_PER_HOUR = 5 + + +def load_gauge_mean( + db_url: Optional[str] = None, + start: Optional[pd.Timestamp] = None, + end: Optional[pd.Timestamp] = None, + box: Sequence[float] = CATCHMENT_BOX, + engine=None, +) -> Optional[pd.Series]: + """Hourly catchment-mean rain_1h (mm) across HII gauges in `box`. + + Pass `engine` (the API's HII store engine) to reuse a pool; otherwise a + connection is resolved from db_url / config. Returns None if the DB is + unavailable or the table is empty. Hours with fewer than + MIN_GAUGES_PER_HOUR reporting gauges are NaN. + """ + if engine is None: + resolved = resolve_db_url(db_url) + if not resolved: + return None + lat_lo, lat_hi, lon_lo, lon_hi = box + try: + from sqlalchemy import create_engine, text + + query = ( + "SELECT m.timestamp, COUNT(m.rain_1h) AS n, AVG(m.rain_1h) AS rain_1h " + "FROM hii_rainfall m JOIN hii_rain_stations s ON s.id = m.station_id " + "WHERE s.latitude BETWEEN :lat_lo AND :lat_hi " + "AND s.longitude BETWEEN :lon_lo AND :lon_hi " + "AND m.rain_1h IS NOT NULL" + ) + params = {"lat_lo": lat_lo, "lat_hi": lat_hi, "lon_lo": lon_lo, "lon_hi": lon_hi} + if start is not None: + query += " AND m.timestamp >= :start" + params["start"] = pd.Timestamp(start).to_pydatetime() + if end is not None: + query += " AND m.timestamp <= :end" + params["end"] = pd.Timestamp(end).to_pydatetime() + query += " GROUP BY m.timestamp ORDER BY m.timestamp" + if engine is None: + engine = create_engine(resolved, pool_pre_ping=True) + with engine.connect() as conn: + frame = pd.read_sql(text(query), conn, params=params) + except Exception as error: + logger.warning(f"HII gauge rain load failed: {error}") + return None + if frame.empty: + return None + frame["timestamp"] = pd.to_datetime(frame["timestamp"]).dt.floor("h") + frame = frame.groupby("timestamp").agg(n=("n", "sum"), rain_1h=("rain_1h", "mean")) + series = pd.to_numeric(frame["rain_1h"], errors="coerce") + series[frame["n"] < MIN_GAUGES_PER_HOUR] = float("nan") + series.name = "hii_gauge_mean" + return series + + +def compare_with_openmeteo( + gauge: pd.Series, openmeteo: pd.Series, window_h: int = 24 +) -> dict: + """Bias/correlation of Open-Meteo against the gauges over the overlap. + + Both are summed over trailing `window_h` so single-hour timing offsets + (gauges report at :00, the model's hour is an interval) do not dominate. + """ + joined = pd.concat( + {"gauge": gauge, "openmeteo": openmeteo}, axis=1 + ).dropna() + if joined.empty: + return {"overlap_hours": 0} + g = joined["gauge"].rolling(window_h, min_periods=window_h).sum() + o = joined["openmeteo"].rolling(window_h, min_periods=window_h).sum() + both = pd.concat({"g": g, "o": o}, axis=1).dropna() + if both.empty: + return {"overlap_hours": int(len(joined))} + return { + "overlap_hours": int(len(joined)), + "window_h": window_h, + "gauge_mean_mm": float(both["g"].mean()), + "openmeteo_mean_mm": float(both["o"].mean()), + "bias_mm": float((both["o"] - both["g"]).mean()), + "mae_mm": float((both["o"] - both["g"]).abs().mean()), + "corr": float(both["g"].corr(both["o"])), + } diff --git a/src/web_api.py b/src/web_api.py index 731caac..2e3239d 100644 --- a/src/web_api.py +++ b/src/web_api.py @@ -1001,6 +1001,84 @@ async def get_hii_rainfall_latest( return rows +@app.get("/api/hii/rainfall/catchment") +async def get_hii_rainfall_catchment( + response: Response, days: int = Query(14, ge=1, le=60) +): + """Upper-Ping catchment-mean hourly rain: HII gauges vs the Open-Meteo + series the flood model actually uses, plus their agreement over the window. + + Evidence-gathering endpoint (docs/FLOOD_FORECASTING.md, HII gauge rain): + the gauge table only exists since 2026-08 so it cannot be a training + feature yet; this makes the two sources' relationship observable meanwhile. + """ + increment_counter("api_requests", labels={"endpoint": "hii_rain_catchment"}) + start = datetime.now() - timedelta(days=days) + + def compute(): + import pandas as pd + + from .ml import hii_rain + + engine = _hii_engine() + if engine is None: + return {"box": hii_rain.CATCHMENT_BOX, "gauge": [], "openmeteo": [], + "comparison_24h_sums": {"overlap_hours": 0}} + gauge = hii_rain.load_gauge_mean(start=pd.Timestamp(start), engine=engine) + openmeteo = None + try: + from sqlalchemy import text + + with engine.connect() as conn: + frame = pd.read_sql( + text( + "SELECT timestamp, catchment_mean FROM openmeteo_rain " + "WHERE timestamp >= :start ORDER BY timestamp" + ), + conn, + params={"start": start}, + ) + if not frame.empty: + frame["timestamp"] = pd.to_datetime(frame["timestamp"]) + openmeteo = pd.to_numeric( + frame.set_index("timestamp")["catchment_mean"], errors="coerce" + ) + except Exception as error: # openmeteo_rain may not exist yet + logger.warning(f"openmeteo_rain read failed: {error}") + + def series_rows(s): + if s is None: + return [] + return [ + {"timestamp": ts.isoformat(), "rain_mm": None if pd.isna(v) else round(float(v), 2)} + for ts, v in s.items() + ] + + comparison = ( + hii_rain.compare_with_openmeteo(gauge, openmeteo) + if gauge is not None and openmeteo is not None + else {"overlap_hours": 0} + ) + return { + "box": hii_rain.CATCHMENT_BOX, + "gauge": series_rows(gauge), + "openmeteo": series_rows(openmeteo), + "comparison_24h_sums": comparison, + } + + payload, stale = await _cached_swr( + HII_CACHE, + HII_CACHE_LOCK, + _HII_COMPUTE_LOCKS["rain"], + f"rain_catchment:{days}", + Config.HII_CACHE_TTL_SECONDS, + compute, + ) + if stale: + response.headers["X-Data-Stale"] = "true" + return payload + + @app.get("/api/hii/waterlevel/latest") async def get_hii_waterlevel_latest( response: Response, hours: int = Query(26, ge=1, le=168) diff --git a/tests/test_hii_collector.py b/tests/test_hii_collector.py index 214dd21..1a53df4 100644 --- a/tests/test_hii_collector.py +++ b/tests/test_hii_collector.py @@ -353,6 +353,29 @@ class TestHiiApiEndpoints: self._get(web_api, "get_hii_rainfall_latest", hours=48) assert calls["n"] == 2 + def test_rainfall_catchment(self, web_api): + """One gauge in the box (CHM005, 19.12N 98.94E) is below the + MIN_GAUGES_PER_HOUR floor, so the catchment mean is NaN -> null, the + openmeteo_rain table does not exist in this store, and the comparison + reports no overlap. Shape is what matters: the endpoint must not 500 + on a fresh database.""" + payload, response = self._get(web_api, "get_hii_rainfall_catchment", days=7) + assert "x-data-stale" not in response.headers + assert list(payload) == ["box", "gauge", "openmeteo", "comparison_24h_sums"] + assert payload["openmeteo"] == [] + assert payload["comparison_24h_sums"] == {"overlap_hours": 0} + assert len(payload["gauge"]) == 1 + assert payload["gauge"][0]["rain_mm"] is None # < MIN_GAUGES_PER_HOUR + + def test_rainfall_catchment_disabled(self, monkeypatch): + from src import web_api + + monkeypatch.setitem(web_api.app_state, "hii_collector", None) + web_api.HII_CACHE.clear() + web_api._REFRESH_IN_FLIGHT.clear() + payload, _ = self._get(web_api, "get_hii_rainfall_catchment", days=7) + assert payload["gauge"] == [] and payload["openmeteo"] == [] + def test_stale_served_on_recompute_failure(self, web_api, monkeypatch): # Prime the cache, expire it, break the DB: the stale copy is served # and flagged via the X-Data-Stale header. diff --git a/tests/test_hii_rain.py b/tests/test_hii_rain.py new file mode 100644 index 0000000..45dd541 --- /dev/null +++ b/tests/test_hii_rain.py @@ -0,0 +1,47 @@ +"""HII gauge-rain aggregate: pure-function tests (no DB).""" + +import numpy as np +import pandas as pd + +from src.ml import hii_rain + + +def _hourly(start, n): + return pd.date_range(start, periods=n, freq="h") + + +def test_compare_identical_series_has_zero_bias(): + idx = _hourly("2026-08-12", 200) + rng = np.random.default_rng(1) + rain = pd.Series(rng.exponential(0.5, len(idx)), index=idx) + out = hii_rain.compare_with_openmeteo(rain, rain.copy(), window_h=24) + assert out["overlap_hours"] == 200 + assert out["bias_mm"] == 0.0 + assert out["mae_mm"] == 0.0 + assert out["corr"] > 0.999 + + +def test_compare_reports_constant_bias(): + idx = _hourly("2026-08-12", 100) + gauge = pd.Series(1.0, index=idx) + model = pd.Series(1.5, index=idx) # model wetter by 0.5 mm/h + out = hii_rain.compare_with_openmeteo(gauge, model, window_h=24) + assert abs(out["bias_mm"] - 12.0) < 1e-9 # 0.5 mm/h x 24 h + + +def test_compare_uses_overlap_only(): + gauge = pd.Series(1.0, index=_hourly("2026-08-12", 100)) + model = pd.Series(1.0, index=_hourly("2026-08-14", 100)) # 52 h overlap + out = hii_rain.compare_with_openmeteo(gauge, model, window_h=24) + assert out["overlap_hours"] == 52 + + +def test_compare_no_overlap(): + gauge = pd.Series(1.0, index=_hourly("2026-01-01", 10)) + model = pd.Series(1.0, index=_hourly("2026-06-01", 10)) + assert hii_rain.compare_with_openmeteo(gauge, model) == {"overlap_hours": 0} + + +def test_load_gauge_mean_without_db_returns_none(monkeypatch): + monkeypatch.setattr(hii_rain, "resolve_db_url", lambda *a, **k: None) + assert hii_rain.load_gauge_mean() is None