Northern-Thailand-Ping-River-Monitor/docs/FLOOD_FORECASTING.md

Flood forecasting

Short-range flood-risk forecasts for the Ping River gauge network, trained on the monitor's own PostgreSQL history. This document covers what the system predicts, what it is built from, how well it actually performs, how to run it on the server, and when to retrain.

Code lives in src/ml/ (data.py, features.py, train.py, predict.py), the training entry point is scripts/train_flood_model.py, tests are in tests/test_flood_forecast.py, and trained artifacts land in models/.

1. Overview

For every station the system answers three questions at three lead times (6, 12 and 24 hours):

• p_warning — probability the water level reaches or exceeds the warning threshold (3.0 m) at any point within the horizon.
• p_danger — same for the danger threshold (4.5 m).
• predicted_max_level — the expected peak level within the horizon, in metres on the station's own datum.

The window is open-ended forward: exceed_warn_6 at 09:00 asks whether the level touches 3.0 m anywhere in (09:00, 15:00], not what it will be at 15:00 — the question an operator actually has.

Fifteen of the sixteen stations have trained models. P.4A (Ban Mae Taeng) is excluded by features.NOT_TRAINABLE (17.2% hourly fill, effectively dead 2019–2024: 289 rows in 2019, 769 in 2024) and is served by the persistence heuristic instead. It still feeds downstream stations as an input, where HistGradientBoosting's native NaN handling copes with the gaps.

Every forecast row carries source (model or heuristic), model_version and trained_at, so a stale or degraded forecast is visible in the payload rather than silently indistinguishable from a good one.

2. Data

Source of record. PostgreSQL table water_measurements joined to stations. As verified on 2026-08-10 (inventory.json, db_cross_check): 592,240 rows, 16 stations, 2018-08-01 through 2026-08-10, zero mismatches against the HTTP API, and a status column that is uniformly active (there are no quality flags to filter on — bad readings must be caught by the feature pipeline, not by the database).

Coverage is the dominant data constraint. Readings are nominally hourly (modal interval 1 h, ~95% of gaps), but only about 56% of hours on the complete hourly grid have a reading: P.1, P.67, P.76 and P.84 at 56.1%, P.103 at 55.8%, P.21 at 55.4%, P.20 at 53.7%, P.87 at 53.3%, P.5 at 50.7%, P.4A at 17.2%.

The missingness is systematic, not random. Measured over the 587 days from 2025-01-01 in models/cache/P.1.csv.gz, the fraction of days with a reading at each hour is roughly 0.80 for 01:00–12:00, 0.68–0.69 for 13:00–16:00, 0.46–0.49 for 17:00–21:00, 0.40–0.42 for 22:00–23:00, and 0.32 at midnight. That is a scrape-schedule fingerprint, not hydrology — and it is why hour-of-day is deliberately not a feature (see section 3).

Long outages would poison training if used naively: P.87 lost 3,961 hours (165 days) in 2023, P.20 lost 2,681 hours in 2021, P.77 lost 2,522 hours in early 2022, P.5 lost 2,429 hours over the 2020–21 turn. features.TRAIN_START excludes P.5 before 2022-01-01; the rest are handled by the per-row coverage gates.

How src/ml/data.py loads it. resolve_db_url() picks a connection string in priority order: an explicit --db-url argument, then the FLOOD_ML_DB_URL environment variable, then Config.get_database_config() when DB_TYPE is postgresql, else None. load_measurements() then tries three tiers:

1. PostgreSQL (_fetch_from_db) — the primary path. NULL discharge stays NULL, which matters because the models must learn from the real missingness pattern.
2. HTTP API (_fetch_from_api, default http://100.81.167.42:8000) — a fallback for running off-server. Caveat: the public history endpoint backfills missing discharge with a synthetic rating-curve estimate, so this path is not equivalent to the DB path. It is flagged as discharge_maybe_synthetic: true in the cache metadata.
3. On-disk cache (models/cache/{station}.csv.gz plus meta.json) — last resort only. A successful DB or API fetch refreshes the cache; the cache is never treated as a source of fresh data.

load_latest() (used by the API) pulls the trailing 336 hours and never writes the cache.

3. Physics and features

Upstream routing

The Ping mainstem gives real forecast skill for free: a flood wave takes hours to travel downstream, so an upstream gauge reading now is information about a downstream gauge later. data-scout measured these travel times by cross-correlating water-level anomalies against the basin anchor P.1. The peak correlation lags, which are hard-coded in features.UPSTREAM_LEADS:

Station	Lead vs P.1	Peak anomaly correlation	Distance to P.1 (km)
P.20 (Ban Chiang Dao)	17 h	0.59	84.5
P.92 (Ban Muang Aut)	15 h	0.66	63.8
P.75 (Ban Chai Lat)	12 h	0.61	45.0
P.4A (Ban Mae Taeng)	12 h	0.73	37.9
P.67 (Ban Tae)	7 h	0.74	25.3
P.21 (Ban Rim Tai)	9 h	0.56	15.0
P.103 (Ring Bridge 3)	1 h	0.88	9.3

(Distances are cumulative straight-line gauge-to-gauge, from the inventory's spatial_order_north_to_south, not channel length — the real river is longer.)

The lags are broadly consistent with distance, with one exception worth knowing about: P.21 is 10 km closer to P.1 than P.67 yet lags by 9 h rather than 7 h, and it has the weakest correlation of the mainstem set (0.56). Whatever the cause, the table encodes the measured lag rather than the one distance would predict — which is the point of measuring instead of assuming.

Two stations are downstream of P.1 (P.5 at −12 h, P.81 at −4 h). For those, UPSTREAM_LEADS routes P.1 and P.103 forward as their inputs, which is the same physics running the other direction.

Six western-tributary stations — P.82, P.84, P.87, P.77, P.85, P.76 — have empty UPSTREAM_LEADS and are un-routed. Their anomaly correlations with P.1 are 0.22–0.32, low enough that routing them would inject noise rather than signal. They are forecast from their own history plus the P.1 basin-state features. This is a known gap: those catchments have no upstream gauge of their own in this network.

Feature set (features.build_features)

Everything is computed on an hourly grid built by make_hourly_grid(), which keeps three aligned frames: observed (raw, NaN where nothing was recorded), filled (forward-filled with FFILL_LIMIT_H = 3), and mask (True where a real reading exists). Features read filled; labels read observed only.

Per target station:

• Self level: current level, lags at 1/2/3/6/12/24/48/72 h.
• Rate of rise: level minus its own value 1/3/6/12/24 h ago — a river at 2.5 m and falling is a different situation from one at 2.5 m rising 30 cm/h.
• Rolling statistics: 6/24 h means, 6/24/72 h maxima, 24 h minimum.
• Discharge: current, lags at 6/24 h, 6 h rise (read from observed, so NULL discharge stays NULL).
• Observation health: obs_age_h (hours since the last real reading, capped at FFILL_LIMIT_H) and cov_24h (fraction of the last 24 h actually observed), so the model can learn to hedge when a gauge is going quiet.
• Routed upstream, per (upstream, lead) pair: the upstream level at lead−3, lead, and lead+3 hours ago, its 6 h rise at lead, and its 24 h rolling max at lead−3. The three-point bracket absorbs error in the measured travel time rather than depending on it being exact.
• Basin state (non-P.1 stations only): P.1 level, its 24 h rolling max, and its 24 h rise.
• Seasonality: doy_sin, doy_cos and an is_monsoon flag for June–October.

Hour-of-day is deliberately excluded. Given the availability profile in section 2, an hour-of-day feature would let the model learn "readings at 03:00 are more likely to exist" and route that through to the label — an artefact of when the scraper runs, with no hydrological content, that would evaporate the moment the scrape schedule changed.

No-leakage guarantees

• Only shift(), backward rolling() and forward-fill are used — nothing interpolates, and no row can read a value timestamped after itself.
• test_no_future_leakage enforces this empirically: it adds +50 m to every reading after time t, rebuilds the features, and asserts the rows at or before t are bit-identical.
• Labels come from observed, never filled, so a forward-filled value can never become its own target.
• The split is strictly temporal, and early_stopping is disabled in HGB_PARAMS specifically because scikit-learn's internal validation split is random and would leak across time.

Coverage gating

A label is only trusted if enough of its forward window was actually observed. build_labels requires MIN_WINDOW_COVERAGE = 0.5 — at least half the horizon's hours present — otherwise the label is NaN and the row is dropped from that head's training set. The one exception is deliberate: an observed exceedance always produces a positive label regardless of coverage, because a confirmed 3.5 m reading inside a sparse window is not ambiguous. Rows whose own features are stale (obs_age_h is NaN, i.e. the last real reading is more than 3 h old) are dropped entirely in build_matrix.

4. Models

Architecture

One HistGradientBoosting model per station × horizon × head:

Head	Type	Target
max_{h}	HistGradientBoostingRegressor (squared error)	rise: max observed level in (t, t+h] minus level at t (v2; serving adds the level back)
warn_{h}	HistGradientBoostingClassifier	level ≥ 3.0 m anywhere in (t, t+h]
danger_{h}	HistGradientBoostingClassifier	level ≥ 4.5 m anywhere in (t, t+h]

Nine heads per station, three horizons (6/12/24 h), fifteen trained stations. Hyperparameters are fixed (HGB_PARAMS: 300 iterations, learning rate 0.06, 31 leaf nodes, minimum 50 samples per leaf, L2 1.0, random_state=42), chosen in an earlier sweep and not re-searched per run — training is deterministic and repeatable.

HistGradientBoosting was chosen for three concrete reasons: it handles NaN natively (essential given ~44% missing hours), it needs no feature scaling, and it trains on CPU alone — no GPU anywhere in this pipeline (measured cost in section 6).

Head gating and fallbacks

The system degrades in tiers rather than failing:

1. Belt-and-braces probability (since 2026-08-11 — see the re-examination note in section 7): the sigmoid-of-regression probability p = 1/(1 + exp(−(predicted_max − threshold)/σ)) is always computed (σ = the regressor's test-residual std, floor MIN_SIGMA = 0.15 m), and when a classifier head exists — trained only if the span had at least MIN_POSITIVES_FOR_CLASSIFIER = 30 positives — the served probability is max(classifier, sigmoid). The classifier can raise the alarm but never silence it: on the gap-filled data a trained classifier stayed near zero through the 2024 record crossing while the regression tracked it.
2. Sigmoid only, when the classifier head is absent or skipped (recorded in skipped_heads with its reason).
3. Persistence heuristic (predict._heuristic_forecast), when there is no model file at all, or the station's newest reading is more than STALE_AFTER_H = 6 hours old. It extrapolates the last 3 h rate of rise forward with a 0.7 damping factor and a fixed σ of 0.3 m. It is not skilful; it exists so the endpoint always returns something structurally valid.

A station is skipped entirely if it has fewer than MIN_ROWS_TO_TRAIN = 200 usable rows; an individual head is skipped below MIN_ROWS_FOR_HEAD = 50 labeled rows. _safe_fit converts any fit failure (typically HistGradientBoosting's binning step rejecting an all-NaN or constant column) into a recorded skip rather than a station-killing exception.

Training procedure

train_station runs two passes. First it evaluates on the strict temporal holdout (train ≤ 2024-12-31, test 2025-01-01 → 2026-08-10) to produce the metrics and the σ calibration. Then it refits every head on the entire record for the deployed artifact, so the shipped model has seen the most recent data. Because the full record has more labeled rows than the training half, the head-gating decisions can differ between the two passes — skipped_heads is therefore re-derived during the refit so it always describes what is actually in the saved bundle, not what the evaluation pass decided.

Bundle format

models/flood_{station}.joblib contains: station_code, model_version (hgb-v1+<git short SHA>), trained_at, sklearn_version, feature_names, horizons, thresholds, heads, sigma, skipped_heads, train_span, and n_train_rows.

feature_names is the important one. At prediction time _model_forecast rebuilds the feature row from live data and checks it against the bundle's stored list; if any expected column is missing it logs an error and falls back to the heuristic rather than feeding scikit-learn silently misaligned columns. test_feature_name_stability guards the same invariant at build time. Bundles are cached in memory keyed by (path, mtime), so dropping in a retrained file invalidates the cache without a restart.

5. Measured performance

Holdout metrics (models/metrics.json)

Two evaluations exist and they differ sharply — the re-examination note in the 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-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	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%

(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 "Head gating"). Recall-at-FAR is null at all horizons on this run. Across stations the 6 h warning PR-AUC now spans 0.987 (P.77) / 0.982 (P.5) / 0.956 (P.85) / 0.950 (P.67) down to 0.436 (P.84), and danger heads are now evaluable at nine stations — strongest P.5 (0.958/0.883/0.811 at 6/12/24 h) and P.77 (0.942/0.863/0.786); P.103's danger metrics, previously the highlight, are null on this span.

Historical evaluation (hgb-v1+49a3de0, 2026-08-10, pre-gap-fill DB — kept for the record; these numbers described the sparser 56%-filled grid and do not reproduce on today's data):

Horizon	Warning PR-AUC	Recall @1% FAR	MAE	Test rows	Base rate
6 h	0.974	98.3%	6.1 cm	8,536	1.36%
12 h	0.904	93.8%	9.0 cm	7,932	1.61%
24 h	0.900	90.1%	11.3 cm	8,572	1.77%

The dramatic PR-AUC difference is mostly the base rate: the filled grid adds ~5,500 quiet test hours per horizon while the number of positive hours barely changes, so the same ranking quality scores far lower — and the classifier's genuine out-of-distribution weakness (see the backtest sections) does the rest.

2026-08-11 re-examination: fuller data changes the backtest story

Read this before the two backtest sections below. On 2026-08-11 the backtests were codified into scripts/backtest_render.py (previously they were one-off runs) and re-run after the database grew from 592k to ~976k rows (a --fill-gaps all pass repaired most of the missing 44% of the hourly grid). Three things changed:

1. The 2024 crossing was 8 hours earlier than documented. The recovered hours show P.1 crossing 3.70 m at 17:00 on 24 September 2024, not 01:00 on 25 September — confirmed independently by the HII sensor at Nawarat Bridge (hii_waterlevel, station 3226: 3.73 m at 17:00). The originally celebrated "24-hour warning" was therefore ~16 hours measured against the real river.
2. Retraining on the fuller data improves level accuracy but degrades the warning classifiers. P.1 24 h MAE improved (11.3 → 10.5 cm), but the warning-head PR-AUC collapsed (0.900 → 0.288 at 24 h): with the filled grid the classifier trains on many more dry-season rows and now stays silent through the September 2024 record crossing while the regression head tracks it. Serving was changed to belt-and-braces — max(classifier, sigmoid(regression)) — so alerting can never be worse than the regression path.
3. Honest current lead times, from the regenerated charts below: the retrained configuration first alerts ~18 h after the true 24 Sep 2024 crossing and roughly at the 27 Sep 2025 crossing. The earlier, better numbers came from models trained and evaluated on the sparser data. The conclusion is not that the old system was better — it is that gauge-only features fundamentally lack lead time for fast rises, which is exactly the rainfall-input and rise-target work now queued (see "Honest limits").

scripts/backtest_render.py regenerates all three charts and fails its acceptance gate while the 2024 lead stays under 12 h — keeping this page honest is now automatic.

2026-08-12 follow-up — hgb-v2 (rise target). A rolling-origin, event-aware evaluation (scripts/evaluate_variants.py, one fold per monsoon 2021-2025) compared the absolute-level target against rise-target variants. The rise target — regression predicts future max minus current level, the level is added back at serving — won decisively and is now deployed as hgb-v2: the regenerated charts below show the 2024 first alert moving from 18 h late to 6 h early (11:00 vs the 17:00 crossing), the 2025 alert from at-crossing to 45 h early, the record-peak underprediction eliminated (the model now slightly overshoots 5.30 m rather than capping ~0.4 m below it), and P.1 MAE improving ~11% at every horizon. Weighted and 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

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 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

The holdout above never sees a true extreme, because the 2024 flood is in the training half. So the model was retrained on data ending 2024-08-31 and asked to forecast September–November 2024 cold, with no knowledge of the event that followed. This is the closest thing to a real operational test available.

The render above shows the whole event hour by hour (regenerated 2026-08-12 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 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.

The hour-by-hour detail of the detection window (regenerated 2026-08-12, 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 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 note above for why the lead times no longer reproduce:

• 25 September cold start. P.1's first warning crossing of the episode was alerted 24–26 hours ahead. This is the genuinely impressive case: the river was in normal state, and the alert came from upstream routing alone.
• 5 October record peak (P.1 5.30 m, P.103 9.93 m — the highest levels in the eight-year record). Alerted 48 hours ahead. Read this one carefully: the river was already in sustained flood by then, so "48 hours" is the _first_alert_at lookback window (lookback_h = 48) saturating, not a measurement of true lead time. The model was correctly alarmed throughout; the metric simply cannot express how much earlier than 48 h that started.
• P.103 danger head over the same window: PR-AUC 0.98–0.99, recall at 1% FAR 87–95%. It called the danger-level crossings, not just the warning ones.
• 8 November re-flood. Caught 26–31 hours ahead by the 12 and 24 h models — a second, independent event in the same test window.
• P.103's 1 September "miss" is a test-boundary artefact: the event begins in the first hours of the test span, before the feature window has enough test-side history to have produced a sustained alert. It is not a model failure, but it is also not evidence of skill.

Honest limits

Genuine lead time is capped by gauge-only physics. The longest upstream travel 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 — 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 reservoir state, though now ingested daily (see docs/DATA_SOURCES.md), measurably hurts alert lead as a model feature — see the 2026-08-13 experiment below.

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 metric is null). The danger head exists and is trained on the full record — the river has spent 57 hours above 4.5 m historically, 0.144% of all hours — but no out-of-sample number backs it. Treat p_danger at P.1 as indicative, not validated.

Thresholds are per-station as of 2026-08-10. THRESHOLDS now carries calibrated (warning, danger) pairs for all 16 stations, derived from the DB's discharge_percent (RID % of channel capacity): warning = median level at 75–85% capacity, danger = median level at 95–105%. P.1 instead uses the official Chiang Mai inundation map (P1_FLOOD_STAGES): warning 3.70 m (city flooding begins, stage 1) and danger 4.20 m (stage 5). The prior single default of (3.0, 4.5) m made P.103 badly over-alert (its bank-full level is ~6.75 m) and P.67 under-alert (overflow at ~2.9 m, 1.6 m below the old danger line). A retrain is required after any threshold change — classifier labels depend on them; until then, model rows report the thresholds baked into their bundle. P.1 additionally reports stages: exceedance probability for each of the seven official inundation stages (3.70–4.60 m), computed from the regression head and its calibration sigma, so they need no retrain and no per-stage classifiers.

2026-08-13: Mae Ngat dam features — a documented negative result

With rid_reservoir_daily backfilled to 2018 (daily Mae Ngat storage/inflow/ outflow, src/ml/dam.py), the obvious v4 experiment was to feed reservoir state to the mainstem models: during the Oct 2024 flood the dam hit 113% of usable capacity with 19–22 MCM/day inflow spikes on the crossing days.

It fails the acceptance gate. On the 2024 record-flood backtest (train < 1 Sep 2024, belt-and-braces alerting, identical to the deployed pipeline):

dam features	first-alert lead	record-peak err (24 h ahead)
none (deployed v3 config)	+13 h (PASS)	+0.24 m
all four	+10 h (FAIL)	+0.22 m
storage % + 3-day delta	+12 h	+0.21…+0.27 m
inflow + outflow	+10 h (FAIL)	+0.35 m
outflow only	+12 h	+0.20 m

Every subset costs 1–3 h of warning for at most a ~3 cm peak-error gain. The mechanism is the publication lag: RID posts the daily report on the morning of its own date (features apply it from 07:00, dam.py's leakage rule), so at the 04:00 first-alert hour of 24 Sep 2024 the freshest dam row still described 23 Sep — a benign reservoir quietly absorbing inflow (outflow 0.13 MCM/day). The columns therefore argue against imminent flooding exactly when the rain features are (correctly) raising the alarm. The rolling-origin harness agrees: rise_rain_dam matches rise_rain on leads and false alarms, only nudging event-peak amplitude (−0.11 → −0.03 m on the Sep 2024 event), and rise_dam (dam without rain) is strictly worse with alarm-latch artifacts.

Disposition: dam features are OFF by default (train_all(use_dam=False); opt-in via --dam on the training CLI, scripts/backtest_render.py --dam, and the rise_rain_dam / rise_dam harness variants). The collector keeps accruing daily rows; revisit post-monsoon when the 2026 season adds dam-era flood events.

Why the lag is probably not the whole story — P.75 already is the dam signal. A 2026-08-13 source sweep put the negative result on firmer ground: P.75 "บ้านช่อแล" sits 3.8 km downstream of the Mae Ngat dam on the Mae Ngat river (nearest other station: P.4A at 10.9 km), it reports hourly, and it has been a model input since v1 with a 12 h routed lead into P.1. Whatever the reservoir releases flows past P.75 within the hour and the model already reads it. The daily reservoir table therefore offers a stale, coarser proxy of a signal the features capture hourly and directly — which is the more likely reason it adds nothing and costs alarm responsiveness.

That reframes what a future intraday source would have to beat: not "no dam information", but "hourly observed dam outflow". Genuine intraday reservoir-state feeds do exist and are open (bigdata-api.rid.go.th SWOC telemetry, and ThaiWater station ridhydro_TUP.16 at the dam), but both are snapshot-only — no archive (verified: the history endpoint returns empty grids for them at every era, including the current one). They can only be accumulated forward, so they cannot retrain against 2024/2025 events. The HII collector already captures both hourly as of 2026-08-11; revisit after the 2026 monsoon, when a season of true intraday reservoir state exists alongside its flood events.

Shipped from the same work: the HII gap-fill merge in the data loader (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.

6. Deployment

API

GET /forecast (src/web_api.py) returns one JSON row per station × horizon with the fields listed in section 1. Results are cached in-process for FORECAST_TTL = 900 seconds (15 minutes), which matches the data cadence — the underlying readings do not update faster than hourly. Inference runs in a thread via asyncio.to_thread so it never blocks the event loop.

Failure modes: 503 if the src.ml package cannot be imported (missing scikit-learn, say), 502 on any other exception.

One behaviour worth knowing, because the code comments suggest otherwise: the endpoint's FileNotFoundError ("No trained flood models found") and RuntimeError handlers are unreachable — nothing in src/ml/ raises either, and predict._forecast_station checks bundle_path.exists() and falls back to the heuristic instead. So before the first training run /forecast returns 200 with an all-heuristic payload, not a 503, provided there is recent gauge data; you get an empty 200 [] only when there is no recent data at all. Judge deployment state by the source field, not the status code.

Dashboard

The "Flood risk outlook" panel (src/static/dashboard.html, loadForecasts()) loads non-blocking after the map renders and stays hidden unless /forecast returns a non-empty array — a non-OK response, an empty array, or a thrown fetch all just leave the panel hidden, and the rest of the dashboard is unaffected. Per the note above, this means the panel appears with heuristic-only content once data is flowing but before any model is trained; the per-chip tooltip is what tells you so. Stations are sorted worst-risk first, each showing three chips (6/12/24 h) coloured by risk band, with the tooltip carrying the exact warning and danger percentages, the predicted peak level, and a "heuristic fallback" note when the row did not come from a model. The panel is labelled experimental.

Training on the server

The server already has the PostgreSQL connection configured, so no host override is needed:

cd /path/to/Northern-Thailand-Ping-River-Monitor
python scripts/train_flood_model.py --stations all

resolve_db_url() picks up Config.get_database_config() automatically when DB_TYPE=postgresql. The run writes fifteen models/flood_{station}.joblib bundles plus models/metrics.json.

Artifacts and dependencies

The fifteen bundles total 101.4 MB — mean 6.76 MB, from 4.04 MB (P.20) to 9.35 MB (P.103 and P.87, with P.5 next at 8.99 MB) — plus metrics.json at 0.34 MB and a 2.1 MB models/cache/. These are not in git, and they should stay that way — artifacts are produced on the server, not shipped. .gitignore excludes models/*.joblib, models/cache/ and models/metrics.json for exactly this reason.

Two pins matter and are already in requirements.txt / pyproject.toml: scikit-learn==1.9.0 and numpy>=1.24,<2 (pandas 2.0.3 wheels are ABI incompatible with numpy 2.x). Bundles record sklearn_version; unpickling a bundle under a different scikit-learn version is not guaranteed to work, so retrain after any scikit-learn upgrade rather than assuming the artifacts carry over.

Measured resource use

All figures below were measured on 2026-08-10 on a development workstation — 24 physical / 32 logical cores at 2.20 GHz, 32 GiB RAM (Python 3.11.9, scikit-learn 1.9.0, joblib 1.5.3, numpy 1.26.4, pandas 2.0.3) — not on the production server. They come from two independent benchmark runs on that same machine, which is why a couple of figures below are quoted as narrow ranges. Treat the CPU times as a floor and the memory figures as representative, since RSS barely depends on core count. Training read the models/cache/ csv.gz files (592,240 rows load in 0.5 s); loading the same history from PostgreSQL was not measured and will be slower.

Training (train_all, all 15 stations, evaluation pass plus full refit):

Measurement	Value
Full 15-station run, unrestricted threads	199 s (3.3 min)
Peak RSS during the full run	209 MB
Single station, unrestricted (P.1 / P.103)	16.5 s / 18.7 s

HistGradientBoosting threads through OpenMP, and it scales only modestly. Timing P.1 alone under OMP_NUM_THREADS:

Threads	1	2	4	unrestricted (32)
P.1 train time	36.5 s	23.5 s	14.7 s	16.5 s

Two things follow. Four threads is the sweet spot — 32 threads was marginally slower than 4, so oversubscription costs you a little. And even one core is enough: at 36.5 s per station, a single-core box retrains all fifteen in roughly 9 minutes (extrapolated, not measured end-to-end).

Per station the fit costs 8–17 s, and P.1 is the worst case at 16.9 s — it is the basin anchor, so it carries 64 features against 32 for stations with fewer upstream inputs (P.85 9.6 s, P.20 8.2 s). Two things are not the cost driver. Evaluation isn't: P.1 with skip_eval=True took 17.3 s, no faster than the full path. Nor is feature engineering — build_matrix over P.1's whole 8-year history is 256 ms against 8–17 s of fitting. The fit is the cost.

One honest caveat about the run that produced the current artifacts. By file mtime it wrote all fifteen models between 11:55:29 and 12:01:15 — 5 min 46 s, averaging 25 s/station including joblib serialization, which lines up with the measured fits. But models/cache/meta.json records the data fetch finishing at 11:45:49, so end to end that run spanned about 15.5 minutes, and the 9 min 40 s gap between fetch and first model could not be reconstructed from the surviving artifacts. Do not attribute it to per-station training cost. Either way the conclusion holds: retraining is minutes, not tens of minutes.

Inference (15 bundles, 16 stations × 3 horizons = 48 rows):

Measurement	Value
Cold call — every bundle unpickled from disk	6.7 s
Warm call — bundles in _MODEL_CACHE	0.72 s median (0.63–0.84 s)
RSS after imports, before any model	71 MB
RSS with all 15 bundles resident	288 MB

The 101.4 MB of on-disk pickles expand to roughly 203–211 MB resident — about 2× — and they stay there: _MODEL_CACHE replaces an entry when the file's mtime changes but never drops one to reclaim memory. That is the single largest memory cost of the whole feature.

Where the time goes: cold start is 5.30 s, of which 0.46 s is the import and 4.84 s is unpickling, and 2.6 s of that is the first bundle alone paying a one-time lazy sklearn.ensemble import — the remaining fourteen average 159 ms. Of the ~640 ms warm compute, model prediction is ~525 ms, the hourly grid 47 ms, and feature building 71 ms across all sixteen stations.

Live-endpoint measurements (a second, independent benchmark run against a real uvicorn instance of the app, same day, same workstation, RSS summed over the process tree):

Measurement	Value
/forecast cache hit (15-min TTL)	2.4 ms median
/forecast cache miss, default threads	15.2 s (≈7 s of that was the HTTP data fallback; a local DB replaces it)
/forecast cache miss, OMP_NUM_THREADS=1	10.9 s
API process RSS, idle → models resident	76 MB → 335 MB

The endpoint-level RSS (335 MB) is higher than the models-only figure above because the live process also retains the pandas frames from the data pull and the HTTP/JSON machinery — use 335 MB as the sizing number.

One threading subtlety cuts the other way in serving: inference is ~135 single-row predicts, and at one row OpenMP thread dispatch costs more than the math — OMP_NUM_THREADS=1 makes the warm compute 2.6× faster (642 ms → 252 ms). Training shows the opposite (2.3× slower single-threaded), so set the variable per process, never globally.

Server sizing, in plain terms: this is a small workload and almost any server runs it. RAM is the binding constraint, not CPU. Budget about 1 GB for the API process so the ~335 MB steady state has headroom on top of the rest of the app; training peaks at only ~210–315 MB and can share the same box. No GPU anywhere. Pin thread counts per process — OMP_NUM_THREADS=1 in the serving unit, OMP_NUM_THREADS=4 for retraining so it cannot monopolise every core while the API is serving. And since a cold call costs seconds against a 2.4 ms cache hit, consider warming /forecast once at startup rather than letting a user absorb it.

7. Retraining policy

Why it matters here specifically. This is not a generic "models go stale" argument:

• Channel geometry changes after every major flood. Scour, deposition and bank failure shift the level-to-discharge relationship at a gauge, and RID revises rating curves after big events. A model trained on the pre-2024 channel is predicting levels for a cross-section that no longer exists.
• Extreme events extend the label range. The highest P.103 reading before the 2024 season was 7.54 m (October 2022); the 2024 event pushed it to 8.27 m on 26 September and 9.93 m on 5 October. Gradient boosting cannot extrapolate past its training range — predictions saturate at the largest value it has seen — so every new record is what makes the next one predictable.
• Station outages change feature availability. P.87's 165-day gap in 2023 and P.4A's five dead years mean the set of populated features drifts over time. Retraining lets head gating and NaN handling re-adapt to the current sensors.

Recommended schedule:

When	Why
Every year, May–June (pre-monsoon)	The minimum. Ensures the model entering the flood season has seen last season in full.
Monthly, July–November	Cheap insurance during the season — a full retrain costs minutes, not hours (section 6), so nice it and forget it.
After any major flood event	Non-negotiable. Channel geometry and rating curves have changed, and the new extreme extends the trainable label range.

Staleness is auditable without guesswork: model_version embeds the git short SHA of the code that trained the bundle (hgb-v1+49a3de0), and trained_at is a timestamp in every bundle. Both are echoed in every /forecast row, so you can tell from the API response alone which code produced a forecast and how old the model is.

8. Operations runbook

All commands assume the project virtualenv is active (.venv locally).

Train (all stations, with evaluation):

python scripts/train_flood_model.py --stations all

Train a subset, refit-only (skips the holdout evaluation — much faster, but produces no metrics and leaves σ at the MIN_SIGMA floor):

python scripts/train_flood_model.py --stations P.1,P.103 --skip-eval

Train from a workstation against the server's database:

export FLOOD_ML_DB_URL='postgresql://user:pass@host:5432/dbname'
python scripts/train_flood_model.py --stations all

Do not commit that URL anywhere. If the DB is unreachable the loader silently falls back to the HTTP API, whose discharge values are partly synthetic — check the log line PostgreSQL fetch failed, falling back to HTTP API before trusting a run.

Verify before promoting. Training writes models/metrics.json alongside the bundles. Check it before treating a run as good:

python -c "import json; m=json.load(open('models/metrics.json')); \
print(m['model_version'], m['split']); \
print({s: v['status'] for s, v in m['stations'].items()}); \
print({h: (d.get('pr_auc_warn'), d.get('mae')) for h, d in m['stations']['P.1']['per_horizon'].items()})"

Expect fifteen trained and one heuristic (P.4A). A station that reports failed names its reason in the same payload. Compare against the previous run's metrics.json, not an absolute bar: after the 2026-08-11 gap-fill the expected baseline is P.1 6 h warning PR-AUC ≈ 0.78 and MAE ≈ 5.5 cm (the historical ~0.97 figure belonged to the sparse pre-fill grid — see section 5). A material drop from the previous run usually means a data problem (a gauge that went quiet, or a bad backfill) rather than a modelling one.

Run the tests (synthetic data only, no database or network required):

python -m pytest tests/test_flood_forecast.py -v

Seven tests covering leakage, label alignment, the coverage gate, forward-fill and staleness, a train/predict round trip, the heuristic fallback, and feature-name stability. The whole suite runs in about 8 seconds, so there is no excuse for skipping it before a deploy.

Understanding graceful degradation. Three things can make a forecast row non-model-backed, and all of them are visible in the payload:

• source: "heuristic", model_version: "heuristic-v1" — either no bundle exists for that station (P.4A always, every station before the first training run), or the station's newest reading is more than 6 hours old.
• A single horizon coming back heuristic while others are model-backed — that horizon's head is in the bundle's skipped_heads, almost always because the station had fewer than 30 positive examples for that threshold.
• A whole station flipping to heuristic after a code change — the feature-name check in _model_forecast caught a mismatch between the live feature builder and the stored feature_names. The fix is to retrain; the log line names the missing columns.

A live example from the 2026-08-10 cache: of 48 forecast rows, 42 came from models and 6 were heuristic — three for P.4A, which has no bundle by design, and three for P.92, whose newest reading was 02:00 while the basin's newest was 09:00. That 7 hours of staleness crossed STALE_AFTER_H = 6, so P.92 correctly dropped to persistence. Both fallback triggers, working as intended, in one ordinary call.

Inspect a bundle's skipped heads directly (joblib.load unpickles, so only ever point it at a bundle this pipeline's own train.py wrote — never a file from elsewhere):

python -c "import joblib; b=joblib.load('models/flood_P.1.joblib'); \
print(b['model_version'], b['trained_at'], b['n_train_rows']); print(b['skipped_heads'])"