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

![Observed vs predicted through the September 2025 flood — model trained only
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 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.

![Observed P.1 level vs the model's 24 h-ahead predicted peak through the
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-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.

![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-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. Removing the publication-lag objection needs an *intraday*
reservoir source; the obvious candidate, `lsim.rid.go.th`, is unreachable
even from a Thai consumer ISP (see `docs/DATA_SOURCES.md`), so that remains
an open search.

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

```bash
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):**

```bash
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):**

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

**Train from a workstation against the server's database:**

```bash
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:

```bash
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):

```bash
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):

```bash
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'])"
```