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)	max observed level in (t, t+h]
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. Classifier head, when the training span contains at least MIN_POSITIVES_FOR_CLASSIFIER = 30 positive examples. Below that, a classifier would be fitting noise, and the head is recorded in skipped_heads with its reason.
2. Sigmoid on the regression head, when the classifier is absent. p = 1/(1 + exp(−(predicted_max − threshold)/σ)), where σ is the standard deviation of the regressor's test residuals (floor MIN_SIGMA = 0.15 m). This turns the peak-level prediction into a calibrated-ish probability that widens correctly when the regressor is less accurate at that horizon — at P.1, σ is 0.15 m at 6 and 12 h but 0.166 m at 24 h.
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)

Model version hgb-v1+49a3de0, generated 2026-08-10. Train ≤ 2024-12-31, test 2025-01-01 → 2026-08-10 — the test span is entirely unseen future data relative to training.

P.1 (Nawarat Bridge), the station that matters most:

Horizon	Warning PR-AUC	Recall @1% FAR	Recall @5% FAR	MAE	MAE above 2 m	Test rows	Base rate
6 h	0.974	98.3%	100%	6.1 cm	9.2 cm	8,536	1.36%
12 h	0.904	93.8%	97.7%	9.0 cm	15.0 cm	7,932	1.61%
24 h	0.900	90.1%	93.4%	11.3 cm	24.5 cm	8,572	1.77%

Read PR-AUC against the base rate — 0.974 versus a 1.36% positive rate is a wide margin over chance. "Recall at 1% false-alarm rate" is the operationally honest number: at a threshold that fires on 1% of quiet hours, the 6 h model still catches 98.3% of warning exceedances.

P.103 (Ring Bridge 3) is the only station with enough danger-level events to evaluate a danger head on the 2025–26 span (base rate 5.7–7.4%): PR-AUC 0.979 / 0.953 / 0.892 and recall at 1% FAR of 97.9% / 89.9% / 79.5% at 6 / 12 / 24 h.

Across the other stations the 6 h warning PR-AUC spans 0.996 (P.5) down to 0.302 (P.82), and tracks almost exactly with how many exceedances that station saw. The strong ones are the frequently-flooded gauges — P.5 0.996, P.81 0.992, P.77 0.968, P.85 0.953, P.75 0.927 — and the weak ones are un-routed western tributaries with almost no positives (P.84 0.570, P.82 0.302 on 0.22% of test hours). P.92 and P.20 have no evaluable warning metric at all: neither crossed 3.0 m often enough in the test span (P.92 not once, P.20 in 0.09% of hours) to score.

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. 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, which is the entire point of the system. Bottom: the model's probability of flooding within 24 h; it fires its first alert at 01:00 on 24 September, a full day before the river crossed the 3.70 m flooding line, stays pinned near 1.0 through both waves, and stands down between and after them. Also visible, honestly: the predicted peak tops out ~0.4 m short of the actual 5.30 m record (the extreme-compression limitation discussed below), and the prediction is noisier on the recession limbs. The same model track drives the dashboard's "Replay Oct 2024 flood" feature, so this chart can be watched live on the map.

The hour-by-hour detail of the detection window shows the sequence exactly: the predicted 24 h peak (amber) starts pulling away from the observed level late on 23 September as upstream gauges rise, the warning probability snaps from ~0 to 1.0 at 01:00 on 24 September, and the river crosses 3.70 m at 01:00 on 25 September — a clean 24-hour warning, delivered while the river in town still looked normal at 2.8 m.

• 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, not routing a wave. The 2025–26 test events bear this out: the 25 September 2025 cold-start crossing was called 7 h ahead by the 12 h model and 9 h ahead by the 24 h model. Practical lead for P.1 is ~7–17 h. Extending it requires rainfall forecasts and Mae Ngat/Mae Kuang dam release data, neither of which this system currently ingests.

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.

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. If P.1's 6 h warning PR-AUC has dropped materially below ~0.97 or its MAE has risen well above ~6 cm, investigate before deploying — that 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'])"