diff --git a/docs/DATA_SOURCES.md b/docs/DATA_SOURCES.md index 36588e1..e01fee2 100644 --- a/docs/DATA_SOURCES.md +++ b/docs/DATA_SOURCES.md @@ -204,7 +204,7 @@ Oct 2024 flood). Mae Kuang Udom Thara is the second upstream reservoir. | Source | What | Access | |---|---|---| -| **Open-Meteo** () | Hourly precip forecast ≤16 days, any lat/lon; **Historical Forecast API archive from 2021** (train on forecast-as-seen, leakage-free); Previous Runs API (fixed 1–7-day leads from Jan 2024); ERA5 back to 1940 | Free, no key, 10k calls/day, non-commercial w/ attribution | +| **Open-Meteo** () — ✅ **INGESTED** since 2026-08-12 (`src/ml/rain.py`): 5 upper-Ping catchment points feed the hgb-v3 model's rain features (trailing sums + forward-24h forecast, archive 2021+); the leader worker also persists hourly rows to the `openmeteo_rain` table | Hourly precip forecast ≤16 days, any lat/lon; **Historical Forecast API archive from 2021** (train on forecast-as-seen, leakage-free); Previous Runs API (fixed 1–7-day leads from Jan 2024); ERA5 back to 1940 | Free, no key, 10k calls/day, non-commercial w/ attribution | | TMD NWP API (`https://data.tmd.go.th/nwpapi/v1/forecast/location/...`) | WRF 4.2 daily/hourly forecasts by place, processed ~06:00 daily | Free Bearer-token registration (`/nwpapi/doc/main/`) | | GFS / ECMWF IFS open data | 0.25° global, 4×/day | Free (NOMADS / AWS / data.ecmwf.int); Open-Meteo already wraps both | diff --git a/docs/FLOOD_FORECASTING.md b/docs/FLOOD_FORECASTING.md index a938cab..979f407 100644 --- a/docs/FLOOD_FORECASTING.md +++ b/docs/FLOOD_FORECASTING.md @@ -261,19 +261,21 @@ next section explains why (the hourly grid was gap-filled from ~56% to ~93% between them, roughly doubling the test rows and collapsing the warning base rates). -**Current model** `hgb-v2` (rise target), generated 2026-08-12 on the -gap-filled DB (~976k rows). Train ≤ 2024-12-31, test 2025-01-01 → 2026-08-12. -P.1: +**Current model** `hgb-v3` (rise target + Open-Meteo rain features), +generated 2026-08-12 on the gap-filled DB (~976k rows). Train ≤ 2024-12-31, +test 2025-01-01 → 2026-08-12. P.1: | Horizon | Warning PR-AUC | MAE | MAE above 2 m | Test rows | Base rate | |---|---|---|---|---|---| -| 6 h | 0.783 | 5.0 cm | 5.2 cm | 14,034 | 0.12% | -| 12 h | 0.508 | 7.2 cm | 10.9 cm | 14,028 | 0.16% | -| 24 h | 0.288 | 9.4 cm | 20.0 cm | 14,034 | 0.25% | +| 6 h | 0.783 | 4.9 cm | 5.0 cm | 14,034 | 0.12% | +| 12 h | 0.508 | 7.2 cm | 10.7 cm | 14,028 | 0.16% | +| 24 h | 0.288 | 8.7 cm | 18.1 cm | 14,034 | 0.25% | -(The prior absolute-target run of the same day, `hgb-v1+d2d0e65`, scored -5.5/8.1/10.5 cm MAE and 7.9/14.6/24.0 cm above 2 m — the rise target improved -every regression figure; PR-AUC belongs to the unchanged classifier heads.) +(Progression across the same day's runs — v1 absolute target: +5.5/8.1/10.5 cm MAE; v2 rise: 5.0/7.2/9.4; v3 rise+rain: 4.9/7.2/8.7 — +with above-2 m MAE falling 24.0 → 20.0 → 18.1 cm at 24 h. PR-AUC belongs to +the unchanged classifier heads; serving is belt-and-braces so alerting uses +the improved regression path regardless.) Level accuracy improved; standalone classifier discrimination did not survive the data change (which is why serving is now `max(classifier, sigmoid)` — see @@ -346,6 +348,21 @@ genuine out-of-distribution weakness (see the backtest sections) does the rest. > quantile variants were evaluated and rejected (more false alarms, no > calibration gain by Brier score). The ≥12 h acceptance gate still fails at > +6 h for 2024 — genuine further lead needs rainfall inputs, not modelling. +> +> **2026-08-12 follow-up 2 — hgb-v3 (rain features): the gate passes.** +> Open-Meteo catchment rainfall (five upper-Ping points, forecast-model +> archive 2021+, `src/ml/rain.py`) added four features: trailing 6/24/72 h +> rain sums and `rain_fc24`, the forward-24 h forecast sum — the first input +> that can act before water reaches any gauge. On the rolling-origin harness +> (`models/eval_rain.json`) rain roughly halved flood-year Brier scores, cut +> flood-regime MAE 20–40%, and moved the hard 2024 leads from +6 h to +11 h +> (P.1) and +10 to +19 h (P.103); the marginal 2025 double-crest event trades +> its artifact +46 h "lead" for a calibrated +2 h with zero false alarms. The +> regenerated backtest below now shows a **13-hour warning for the 2024 +> record flood (alert 04:00, crossing 17:00) — the ≥12 h acceptance gate +> passes for the first time**. P.1 MAE improves again to 4.9/7.2/8.7 cm at +> 6/12/24 h. Serving fetches live rain hourly and degrades to NaN features +> (never a crash) if Open-Meteo is unreachable. ### The September 2025 flood, as the deployed configuration saw it @@ -355,12 +372,13 @@ through 2024](img/backtest-2025-p1.png) This uses the deployed configuration (train ≤ 2024-12-31) on an event it never saw. *(Chart regenerated 2026-08-12 with the hgb-v2 rise target on the gap-filled data — see the re-examination note above for the history of these -numbers.)* The v2 model first alerts at **21:00 on 25 September 2025 — 45 -hours before the river crosses 3.70 m** at 18:00 on the 27th. The long lead is -partly the twin-crest shape of this event (the near-miss 3.51 m crest of -26 September keeps the 24 h-peak forecast near the line), so read it as -"the model was correctly alarmed through the whole double crest", not as a -general 45 h capability. +numbers.)* The v3 model first alerts at **16:00 on 27 September 2025 — 2 hours before +the river crosses 3.70 m** at 18:00. This is a shorter lead than v2's 45 h, +and deliberately so: v2's long "lead" was an alarm that latched through the +near-miss 3.51 m crest of the 26th; v3's rain-informed probabilities are far +better calibrated on this marginal event (Brier halved, zero false-alarm +episodes on the season) and fire when exceedance actually becomes likely. A +barely-over-threshold crest is intrinsically a short-notice event. ### Headline validation: the October 2024 record flood @@ -373,27 +391,29 @@ followed. This is the closest thing to a real operational test available. October 2024 flood, with the warning probability below](img/backtest-2024-p1.png) The render above shows the whole event hour by hour *(regenerated 2026-08-12 -with the hgb-v2 rise target)*. Top: the observed level (blue) against the -24 h-ahead predicted peak the model issued at each hour (amber, dashed) — the -amber line now leads the blue one into both flood waves. Bottom: the +with the hgb-v3 rise + rain configuration)*. Top: the observed level (blue) +against the 24 h-ahead predicted peak the model issued at each hour (amber, +dashed) — the amber line leads the blue one into both flood waves. Bottom: the belt-and-braces probability of flooding within 24 h; the **first alert comes -at 11:00 on 24 September, 6 hours before the true 17:00 crossing**, and stays -correctly alarmed through the October record wave. The predicted peak now -slightly overshoots the 5.30 m record instead of capping ~0.4 m below it — the -rise target removed the cannot-exceed-training-max ceiling. The same historic -model track drives the dashboard's "Replay Oct 2024 flood" feature. +at 04:00 on 24 September, 13 hours before the true 17:00 crossing**, while +the river in town still read 2.9 m — the rain features react to upstream +precipitation before any gauge rises. The rise target removed the +cannot-exceed-training-max ceiling, so the record 5.30 m peak is tracked +rather than capped. The same historic model track drives the dashboard's +"Replay Oct 2024 flood" feature. ![Hour-by-hour detail of the detection window, 22–28 September 2024](img/backtest-2024-p1-detail.png) The hour-by-hour detail of the detection window *(regenerated 2026-08-12, -hgb-v2)* shows the sequence: the river crosses 3.70 m at **17:00 on +hgb-v3)* shows the sequence: the river crosses 3.70 m at **17:00 on 24 September** (the hours recovered by gap-filling; independently confirmed by the HII sensor at the same bridge), and the model's probability crosses 0.5 at -**11:00 — a 6-hour warning** delivered while the river stood at 3.4 m. Under -the absolute-level target this alert came 18 hours *after* the crossing; the -rise target recovered the lead. Extending 6 h toward the 12+ h the acceptance -gate demands is what the rainfall features are for. +**04:00 — a 13-hour warning** delivered while the river stood at 2.9 m. The +same alert under the absolute-level target came 18 hours *after* the crossing +(v1), and 6 hours before it with the rise target alone (v2); catchment +rainfall closed the rest. This clears the ≥12 h acceptance gate in +`scripts/backtest_render.py`. The event bullets below quote the original (pre-gap-fill) evaluation of the deployed model and are kept for the historical record — see the re-examination @@ -423,14 +443,16 @@ note above for why the lead times no longer reproduce: time into P.1 is 17 h (P.20), and the strongest predictors are much closer: P.103 at 1 h, P.67 at 7 h, P.21 at 9 h. Once a 24 h forecast reaches past roughly 17 h, there is no observation that has "already happened" to inform it — the model -is extrapolating basin state and season, not routing a wave. The regenerated -backtests bear this out — harder than first documented (see the re-examination -note in section 7): on the gap-filled data the retrained configuration alerts -the 27 September 2025 crossing *as it happens* and the 24 September 2024 -crossing ~18 h *late*. **Genuine gauge-only lead for P.1 is at best ~7–17 h, -and for fast rises can be zero.** Extending it requires rainfall inputs and -Mae Ngat/Mae Kuang dam release data, plus the rise-target/quantile modelling -work — rainfall collection began 2026-08-11 (see `docs/DATA_SOURCES.md`). +is extrapolating basin state and season — unless it has rainfall. That is no +longer hypothetical: hgb-v3's Open-Meteo features (see the re-examination +notes in section 7) took the 2024 record-flood lead from 18 h late (v1 +gauge-only, absolute target) to 13 h early, precisely because catchment rain +acts before any gauge rises, and `rain_fc24` — a weather *forecast* — acts +before the rain itself falls. **Remaining honest limits:** marginal +just-over-threshold crests (2025: +2 h) are intrinsically short-notice; the +rain series only exists from 2021-03, so older training rows are rain-blind; +forecast-rain quality bounds what the feature can add; and Mae Ngat/Mae Kuang +dam releases remain uningested (see `docs/DATA_SOURCES.md`). **Danger-level skill at P.1 is unproven.** P.1 never crossed 4.5 m in the 2025-01-01 → 2026-08-10 test span (`base_rate_danger` is 0.0, so every danger diff --git a/docs/img/backtest-2024-p1-detail.png b/docs/img/backtest-2024-p1-detail.png index e105297..f960de3 100644 Binary files a/docs/img/backtest-2024-p1-detail.png and b/docs/img/backtest-2024-p1-detail.png differ diff --git a/docs/img/backtest-2024-p1.png b/docs/img/backtest-2024-p1.png index 2b5b910..0052032 100644 Binary files a/docs/img/backtest-2024-p1.png and b/docs/img/backtest-2024-p1.png differ diff --git a/docs/img/backtest-2025-p1.png b/docs/img/backtest-2025-p1.png index d752709..21225f1 100644 Binary files a/docs/img/backtest-2025-p1.png and b/docs/img/backtest-2025-p1.png differ diff --git a/models/eval_rain.json b/models/eval_rain.json new file mode 100644 index 0000000..5f48820 --- /dev/null +++ b/models/eval_rain.json @@ -0,0 +1,439 @@ +[ + { + "station": "P.1", + "warn_thr": 3.7, + "folds": [ + { + "year": 2021, + "n_train": 20024, + "n_test": 4392, + "events": [], + "variants": { + "rise": { + "mae": 0.070970681677648, + "mae_above_2p5": null, + "brier_warn": 0.0, + "events": [], + "false_alarm_episodes": 0 + }, + "rise_rain": { + "mae": 0.07247641391827915, + "mae_above_2p5": null, + "brier_warn": 0.0, + "events": [], + "false_alarm_episodes": 0 + } + } + }, + { + "year": 2022, + "n_train": 28782, + "n_test": 4392, + "events": [ + { + "crossing": "2022-10-02T19:00:00", + "peak_ts": "2022-10-03T15:00:00", + "peak_level": 4.65 + } + ], + "variants": { + "rise": { + "mae": 0.08489734638010188, + "mae_above_2p5": 0.24517428534843191, + "brier_warn": 0.004136576477574926, + "events": [ + { + "crossing": "2022-10-02T19:00:00", + "lead_h": 6.0, + "peak_level": 4.65, + "peak_pred_24h_before": 3.870617057762467 + } + ], + "false_alarm_episodes": 0 + }, + "rise_rain": { + "mae": 0.08491621780630447, + "mae_above_2p5": 0.24539225527181602, + "brier_warn": 0.004209435091817237, + "events": [ + { + "crossing": "2022-10-02T19:00:00", + "lead_h": 5.0, + "peak_level": 4.65, + "peak_pred_24h_before": 3.844081593978701 + } + ], + "false_alarm_episodes": 0 + } + } + }, + { + "year": 2023, + "n_train": 37537, + "n_test": 4392, + "events": [], + "variants": { + "rise": { + "mae": 0.07647507344443248, + "mae_above_2p5": 0.12593127745781565, + "brier_warn": 3.19529506149809e-11, + "events": [], + "false_alarm_episodes": 0 + }, + "rise_rain": { + "mae": 0.07491688056430452, + "mae_above_2p5": 0.09506457784884237, + "brier_warn": 6.639993512164624e-14, + "events": [], + "false_alarm_episodes": 0 + } + } + }, + { + "year": 2024, + "n_train": 46321, + "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": { + "mae": 0.10075626719251345, + "mae_above_2p5": 0.26384939692286363, + "brier_warn": 0.011608221543810462, + "events": [ + { + "crossing": "2024-09-24T17:00:00", + "lead_h": 6.0, + "peak_level": 4.93, + "peak_pred_24h_before": 4.534709676862131 + }, + { + "crossing": "2024-10-03T09:00:00", + "lead_h": 3.0, + "peak_level": 5.3, + "peak_pred_24h_before": 5.056619694027486 + } + ], + "false_alarm_episodes": 0 + }, + "rise_rain": { + "mae": 0.08874847709708966, + "mae_above_2p5": 0.20967242138010514, + "brier_warn": 0.005862181747890141, + "events": [ + { + "crossing": "2024-09-24T17:00:00", + "lead_h": 11.0, + "peak_level": 4.93, + "peak_pred_24h_before": 4.863137825109792 + }, + { + "crossing": "2024-10-03T09:00:00", + "lead_h": 21.0, + "peak_level": 5.3, + "peak_pred_24h_before": 5.557547530221961 + } + ], + "false_alarm_episodes": 0 + } + } + }, + { + "year": 2025, + "n_train": 55081, + "n_test": 4392, + "events": [ + { + "crossing": "2025-09-27T18:00:00", + "peak_ts": "2025-09-27T22:00:00", + "peak_level": 3.93 + } + ], + "variants": { + "rise": { + "mae": 0.13102686839694494, + "mae_above_2p5": 0.2729944665011028, + "brier_warn": 0.010889835530782944, + "events": [ + { + "crossing": "2025-09-27T18:00:00", + "lead_h": 46.0, + "peak_level": 3.93, + "peak_pred_24h_before": 3.232029710676319 + } + ], + "false_alarm_episodes": 1 + }, + "rise_rain": { + "mae": 0.11160288394028639, + "mae_above_2p5": 0.24838881498904114, + "brier_warn": 0.005150122823202273, + "events": [ + { + "crossing": "2025-09-27T18:00:00", + "lead_h": 2.0, + "peak_level": 3.93, + "peak_pred_24h_before": 3.239754736790134 + } + ], + "false_alarm_episodes": 0 + } + } + } + ] + }, + { + "station": "P.103", + "warn_thr": 5.95, + "folds": [ + { + "year": 2021, + "n_train": 20009, + "n_test": 4392, + "events": [], + "variants": { + "rise": { + "mae": 0.15279228710793116, + "mae_above_2p5": null, + "brier_warn": 0.0, + "events": [], + "false_alarm_episodes": 0 + }, + "rise_rain": { + "mae": 0.15114584578006643, + "mae_above_2p5": null, + "brier_warn": 0.0, + "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": { + "mae": 0.15389446229226025, + "mae_above_2p5": 0.41887094569113953, + "brier_warn": 0.007740077230532647, + "events": [ + { + "crossing": "2022-08-14T04:00:00", + "lead_h": 6.0, + "peak_level": 6.09, + "peak_pred_24h_before": 4.837975953559253 + }, + { + "crossing": "2022-10-02T16:00:00", + "lead_h": 9.0, + "peak_level": 7.54, + "peak_pred_24h_before": 7.168225721504108 + } + ], + "false_alarm_episodes": 0 + }, + "rise_rain": { + "mae": 0.15303667736579823, + "mae_above_2p5": 0.38833713105646916, + "brier_warn": 0.0074511589478895475, + "events": [ + { + "crossing": "2022-08-14T04:00:00", + "lead_h": 6.0, + "peak_level": 6.09, + "peak_pred_24h_before": 4.8804746828604815 + }, + { + "crossing": "2022-10-02T16:00:00", + "lead_h": 10.0, + "peak_level": 7.54, + "peak_pred_24h_before": 7.226685294045882 + } + ], + "false_alarm_episodes": 0 + } + } + }, + { + "year": 2023, + "n_train": 37524, + "n_test": 4392, + "events": [], + "variants": { + "rise": { + "mae": 0.16106550376855028, + "mae_above_2p5": null, + "brier_warn": 0.0001779564366636896, + "events": [], + "false_alarm_episodes": 1 + }, + "rise_rain": { + "mae": 0.15541158498896382, + "mae_above_2p5": null, + "brier_warn": 4.103077380009017e-09, + "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": { + "mae": 0.17314177863843333, + "mae_above_2p5": 0.48294427141283425, + "brier_warn": 0.016091510788709233, + "events": [ + { + "crossing": "2024-09-24T10:00:00", + "lead_h": 10.0, + "peak_level": 8.27, + "peak_pred_24h_before": 7.167139790234238 + }, + { + "crossing": "2024-09-30T03:00:00", + "lead_h": 10.0, + "peak_level": 5.99, + "peak_pred_24h_before": 5.618725525893519 + }, + { + "crossing": "2024-10-03T06:00:00", + "lead_h": 69.0, + "peak_level": 9.93, + "peak_pred_24h_before": 7.889819144742388 + } + ], + "false_alarm_episodes": 0 + }, + "rise_rain": { + "mae": 0.1400555603272253, + "mae_above_2p5": 0.2920105854094151, + "brier_warn": 0.007294262727198402, + "events": [ + { + "crossing": "2024-09-24T10:00:00", + "lead_h": 19.0, + "peak_level": 8.27, + "peak_pred_24h_before": 8.04142923647258 + }, + { + "crossing": "2024-09-30T03:00:00", + "lead_h": 10.0, + "peak_level": 5.99, + "peak_pred_24h_before": 5.42 + }, + { + "crossing": "2024-10-03T06:00:00", + "lead_h": 55.0, + "peak_level": 9.93, + "peak_pred_24h_before": 8.655375706947945 + } + ], + "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": { + "mae": 0.21617454799997243, + "mae_above_2p5": 0.406245211140301, + "brier_warn": 0.012587550711857222, + "events": [ + { + "crossing": "2025-09-26T06:00:00", + "lead_h": 14.0, + "peak_level": 6.64, + "peak_pred_24h_before": 5.73 + }, + { + "crossing": "2025-10-03T06:00:00", + "lead_h": 47.0, + "peak_level": 6.14, + "peak_pred_24h_before": 5.953951787788988 + } + ], + "false_alarm_episodes": 1 + }, + "rise_rain": { + "mae": 0.17651730883692462, + "mae_above_2p5": 0.38926658786456825, + "brier_warn": 0.015873258461498164, + "events": [ + { + "crossing": "2025-09-26T06:00:00", + "lead_h": 6.0, + "peak_level": 6.64, + "peak_pred_24h_before": 5.760150202082536 + }, + { + "crossing": "2025-10-03T06:00:00", + "lead_h": 8.0, + "peak_level": 6.14, + "peak_pred_24h_before": 5.605274163405912 + } + ], + "false_alarm_episodes": 2 + } + } + } + ] + } +] \ No newline at end of file diff --git a/scripts/backtest_render.py b/scripts/backtest_render.py index e93e96e..987243f 100644 --- a/scripts/backtest_render.py +++ b/scripts/backtest_render.py @@ -48,13 +48,16 @@ RED = "#d9534f" def fit_backtest_model(df_long: pd.DataFrame, train_end: str): """Train the 24 h regression + warning heads on rows <= train_end only. - Mirrors the deployed hgb-v2 pipeline: the regression head learns the RISE - over the current level (rolling-origin evaluation 2026-08-12 showed this - moves first-alert leads from ~0 h to +6..+46 h); label statistics are - bounded to the training cutoff. + Mirrors the deployed hgb-v3 pipeline: the regression head learns the RISE + over the current level, with Open-Meteo catchment-rain features (trailing + sums + the forward-24h forecast sum); label statistics are bounded to the + training cutoff. """ + from src.ml import rain as rain_mod + + rain_series = rain_mod.catchment_mean(rain_mod.load_history()) X, Y, _meta = features.build_matrix( - df_long, STATION, (HORIZON,), stats_end=train_end + df_long, STATION, (HORIZON,), stats_end=train_end, rain=rain_series ) train_mask = X.index <= pd.Timestamp(train_end) X_train, Y_train = X.loc[train_mask], Y.loc[train_mask] diff --git a/src/ml/predict.py b/src/ml/predict.py index 36efb6e..e1582ee 100644 --- a/src/ml/predict.py +++ b/src/ml/predict.py @@ -127,6 +127,7 @@ def _model_forecast( bundle: dict, as_of: pd.Timestamp, current_level: float, + rain: Optional[pd.Series] = None, ) -> List[dict]: warn_thr = bundle["thresholds"]["warning"] danger_thr = bundle["thresholds"]["danger"] @@ -145,7 +146,7 @@ def _model_forecast( ) warn_thr, danger_thr = cfg_warn, cfg_danger - feature_row = features.build_features(grid, station_code).loc[[as_of]] + feature_row = features.build_features(grid, station_code, rain=rain).loc[[as_of]] expected_columns = bundle["feature_names"] missing = [c for c in expected_columns if c not in feature_row.columns] if missing: @@ -229,6 +230,7 @@ def _forecast_station( models_dir: Path, now: pd.Timestamp, horizons: Tuple[int, ...], + rain: Optional[pd.Series] = None, ) -> List[dict]: level_col = (station_code, "water_level") if level_col not in grid.observed.columns: @@ -264,7 +266,9 @@ def _forecast_station( ) bundle = _load_bundle(bundle_path) - model_results = _model_forecast(station_code, grid, bundle, as_of, current_level) + model_results = _model_forecast( + station_code, grid, bundle, as_of, current_level, rain=rain + ) if model_results is None: return _heuristic_forecast( station_code, @@ -301,6 +305,7 @@ def get_forecasts( readings_by_station: Dict[str, List[dict]], models_dir: Union[str, Path] = DEFAULT_MODELS_DIR, now: Optional[Union[datetime.datetime, str]] = None, + rain: Optional[pd.Series] = None, ) -> List[dict]: """Produce flood forecasts for every station present in `readings_by_station`. @@ -323,7 +328,9 @@ def get_forecasts( for station_code in readings_by_station.keys(): try: results.extend( - _forecast_station(station_code, grid, models_dir, now, DEFAULT_HORIZONS) + _forecast_station( + station_code, grid, models_dir, now, DEFAULT_HORIZONS, rain=rain + ) ) except Exception as error: logger.error(f"Forecast failed for station {station_code}: {error}") @@ -359,4 +366,14 @@ def get_latest_forecasts( f"No recent data for station {missing_station}; omitting from forecasts" ) - return get_forecasts(readings_by_station, models_dir=models_dir) + # Live rain: trailing days + next-48h forecast. On fetch failure pass an + # EMPTY series (not None) so rain-trained bundles still find their columns + # (as NaN) and serve model output instead of tripping the feature guard. + from .rain import serving_series + + rain = serving_series() + if rain is None: + logger.warning("live rain unavailable; rain features will be NaN") + rain = pd.Series(dtype=float) + + return get_forecasts(readings_by_station, models_dir=models_dir, rain=rain) diff --git a/src/ml/rain.py b/src/ml/rain.py index 7182404..9a77e8c 100644 --- a/src/ml/rain.py +++ b/src/ml/rain.py @@ -160,3 +160,60 @@ def serving_series() -> Optional[pd.Series]: except Exception as error: logger.warning(f"Open-Meteo forecast fetch failed: {error}") return None + + +def save_to_db(df: pd.DataFrame, engine, db_type: str) -> int: + """Upsert per-point + catchment-mean hourly rain into openmeteo_rain. + + Called by the leader worker's hourly precompute with the live forecast + frame, so the DB accumulates both what fell (past rows are the model + analysis) and what was forecast (future rows, overwritten as they become + past). The ML training path reads Open-Meteo's own archive, not this + table — this is for dashboards, SQL analysis, and source independence. + """ + if df is None or df.empty: + return 0 + from sqlalchemy import text + + point_cols = [p[0] for p in CATCHMENT_POINTS] + ddl_cols = ", ".join(f"{c} NUMERIC(6,2)" for c in point_cols) + ddl = ( + "CREATE TABLE IF NOT EXISTS openmeteo_rain (" + "timestamp TIMESTAMP PRIMARY KEY, " + f"{ddl_cols}, catchment_mean NUMERIC(6,2), " + "created_at TIMESTAMP DEFAULT CURRENT_TIMESTAMP)" + ) + cols = ["timestamp"] + point_cols + ["catchment_mean"] + placeholders = ", ".join(f":{c}" for c in cols) + updates = ", ".join( + f"{c} = " + + (f"VALUES({c})" if db_type == "mysql" else f"EXCLUDED.{c}") + for c in cols[1:] + ) + if db_type == "mysql": + sql = ( + f"INSERT INTO openmeteo_rain ({', '.join(cols)}) VALUES ({placeholders}) " + f"ON DUPLICATE KEY UPDATE {updates}" + ) + else: + sql = ( + f"INSERT INTO openmeteo_rain ({', '.join(cols)}) VALUES ({placeholders}) " + f"ON CONFLICT (timestamp) DO UPDATE SET {updates}" + ) + mean = df.mean(axis=1) + params = [ + { + "timestamp": ts.to_pydatetime(), + **{c: (None if pd.isna(row[c]) else float(row[c])) for c in point_cols}, + "catchment_mean": None if pd.isna(mean.loc[ts]) else float(mean.loc[ts]), + } + for ts, row in df.iterrows() + ] + try: + with engine.begin() as conn: + conn.execute(text(ddl)) + conn.execute(text(sql), params) + return len(params) + except Exception as error: + logger.error(f"openmeteo_rain save failed: {error}") + return 0 diff --git a/src/ml/train.py b/src/ml/train.py index 092111a..caf8aa7 100644 --- a/src/ml/train.py +++ b/src/ml/train.py @@ -203,9 +203,10 @@ def train_station( split_train_end: str = SPLIT_B_TRAIN_END, split_test_start: str = SPLIT_B_TEST_START, split_test_end: str = SPLIT_B_TEST_END, + rain: Optional[pd.Series] = None, ) -> Tuple[Optional[dict], dict]: """Train every head for one station. Returns (bundle_or_None, station_metrics).""" - X, Y, meta = features.build_matrix(df_long, station, horizons) + X, Y, meta = features.build_matrix(df_long, station, horizons, rain=rain) if meta["n_rows"] < MIN_ROWS_TO_TRAIN: return None, { "status": "failed", @@ -411,10 +412,12 @@ def train_station( ] = f"only {n_pos} positives in train span (< {MIN_POSITIVES_FOR_CLASSIFIER})" final_heads[head_key] = None + # v3 = rise target + Open-Meteo rain features; v2 = rise target only + version_prefix = "hgb-v3" if "rain_24h" in feature_names else "hgb-v2" bundle = { "station_code": station, - "model_version": f"hgb-v2+{_git_short_sha()}", - # v2: regression heads predict the RISE over the current level; the + "model_version": f"{version_prefix}+{_git_short_sha()}", + # v2+: regression heads predict the RISE over the current level; the # serving side must add the level back. Old v1 bundles lack this key. "regression_target": "rise", "trained_at": datetime.datetime.now().isoformat(), @@ -439,11 +442,28 @@ def train_all( models_dir: Path = Path("models"), skip_eval: bool = False, hgb_overrides: Optional[dict] = None, + use_rain: bool = True, ) -> dict: """Train and save every requested station's models. Returns the metrics.json payload.""" models_dir = Path(models_dir) models_dir.mkdir(parents=True, exist_ok=True) - model_version = f"hgb-v2+{_git_short_sha()}" + + # Catchment rain (Open-Meteo archive, 2021+). Optional: without it the + # models train as v2 (no rain columns) and still serve correctly. + rain_series = None + if use_rain: + try: + from . import rain as rain_mod + + rain_series = rain_mod.catchment_mean(rain_mod.load_history()) + except Exception as error: + logger.warning(f"rain history unavailable, training without it: {error}") + if rain_series is not None: + logger.info( + f"rain series: {rain_series.index.min()} .. {rain_series.index.max()}" + ) + version_prefix = "hgb-v3" if rain_series is not None else "hgb-v2" + model_version = f"{version_prefix}+{_git_short_sha()}" station_results: Dict[str, dict] = {} for station in stations: @@ -459,6 +479,7 @@ def train_all( horizons, skip_eval=skip_eval, hgb_overrides=hgb_overrides, + rain=rain_series, ) if bundle is None: logger.warning(f"{station}: failed ({station_metrics.get('reason')})") @@ -516,6 +537,11 @@ def main(argv: Optional[List[str]] = None) -> None: parser.add_argument( "--end", default=None, help="ISO date; latest measurement to load" ) + parser.add_argument( + "--no-rain", + action="store_true", + help="train without the Open-Meteo rain features (v2-style bundles)", + ) args = parser.parse_args(argv) if args.stations == "all": @@ -539,7 +565,11 @@ def main(argv: Optional[List[str]] = None) -> None: ) metrics_payload = train_all( - df_long, stations, models_dir=Path(args.models_dir), skip_eval=args.skip_eval + df_long, + stations, + models_dir=Path(args.models_dir), + skip_eval=args.skip_eval, + use_rain=not args.no_rain, ) trained = sum( 1 for s in metrics_payload["stations"].values() if s["status"] == "trained" diff --git a/src/web_api.py b/src/web_api.py index 801b29e..eae5dae 100644 --- a/src/web_api.py +++ b/src/web_api.py @@ -264,6 +264,27 @@ app.add_middleware( ) +async def _persist_rain(): + """Save the latest Open-Meteo rain frame into openmeteo_rain (leader only).""" + store = app_state.get("forecast_store") # reuse its SQL engine + if not store: + return + try: + from .ml import rain as rain_mod + + def fetch_and_save(): + frame = rain_mod.fetch_forecast() + if not store.engine and not store.connect(): + return 0 + return rain_mod.save_to_db(frame, store.engine, store.db_type) + + saved = await asyncio.to_thread(fetch_and_save) + if saved: + logger.info(f"openmeteo_rain: {saved} hourly rows upserted") + except Exception as e: + logger.warning(f"rain persistence failed: {e}") + + async def _precompute_forecasts(): """Refresh the forecast cache and persist the issued forecasts (leader only).""" try: @@ -351,6 +372,10 @@ async def background_scraping_task(): except Exception as e: logger.error(f"HII collection failed: {e}") + # Persist the Open-Meteo catchment rain (observed tail + + # 48h forecast) so the DB carries the weather context too + await _persist_rain() + # Precompute forecasts on fresh data: primes the response # cache (user requests never pay for inference) and records # what the model predicted for later predicted-vs-actual diff --git a/tests/test_flood_forecast.py b/tests/test_flood_forecast.py index 1377756..6e66cfa 100644 --- a/tests/test_flood_forecast.py +++ b/tests/test_flood_forecast.py @@ -197,7 +197,7 @@ def test_train_smoke_and_roundtrip(tmp_path): df = make_synth(n, data_stations, seed=7, pulses=pulses) metrics = train.train_all( - df, target_stations, models_dir=tmp_path, skip_eval=True, hgb_overrides={"max_iter": 20} + df, target_stations, models_dir=tmp_path, skip_eval=True, hgb_overrides={"max_iter": 20}, use_rain=False ) assert metrics["stations"]["P.1"]["status"] == "trained" assert metrics["stations"]["P.20"]["status"] == "trained" @@ -237,7 +237,7 @@ def test_feature_name_stability(tmp_path): upstream = [code for code, _lead in features.UPSTREAM_LEADS["P.1"]] data_stations = ["P.1"] + upstream df = make_synth(300, data_stations, seed=11, pulses={"P.1": [(100, 20, 2.0)]}) - train.train_all(df, ["P.1"], models_dir=tmp_path, skip_eval=True, hgb_overrides={"max_iter": 10}) + train.train_all(df, ["P.1"], models_dir=tmp_path, skip_eval=True, hgb_overrides={"max_iter": 10}, use_rain=False) # Safe: loading the bundle this same test just wrote to tmp_path, not an external file. bundle = joblib.load(tmp_path / "flood_P.1.joblib")