grabowski
6d4beae65a
CI / Format & lint (push) Successful in 24s
CI / Test suite (push) Successful in 48s
Docs / Validate documentation (push) Successful in 24s
Security / Dependency vulnerabilities (push) Successful in 2m20s
Security / Static analysis (push) Successful in 21s
Security / License report (push) Successful in 1m40s
build_features cut the rain series to the gauge grid before the forward sum, and the grid ends at the newest reading, so the row every live forecast is made from had rain_fc24 = NaN. Training and all evaluations had the real Open-Meteo forecast there. The forward sum now runs over the grid extended by 24 h; earlier rows are unchanged. Harness (new opt-in variant rise_rain_fc24_missing = rise_rain with rain_fc24 blanked on test rows, five folds, --from-cache): P.1 and P.103 nearly unaffected, P.21 false alarms 10 -> 5 and MAE 0.232 -> 0.215 with the fix. Write-up in docs/FLOOD_FORECASTING.md section 5; results in models/eval_2026-10-18_fc24.json. Also: HII mirror rows in one hour now resolve to the latest reading via a stable sort on the raw time (was an unstable sort on the floored hour, i.e. row-order dependent). Needed for an exact Rust port.
443 lines
17 KiB
Python
443 lines
17 KiB
Python
"""Static config and feature/label engineering for the Ping River flood forecast models.
|
|
|
|
All feature computation is strictly causal (no row uses information timestamped after
|
|
itself) so it is safe to run identically at training time and at prediction time.
|
|
"""
|
|
|
|
import datetime
|
|
import logging
|
|
from dataclasses import dataclass
|
|
from typing import Dict, List, Optional, Tuple
|
|
|
|
import numpy as np
|
|
import pandas as pd
|
|
|
|
logger = logging.getLogger(__name__)
|
|
|
|
# ---------------------------------------------------------------------------
|
|
# Static configuration
|
|
# ---------------------------------------------------------------------------
|
|
|
|
# Per-station (warning, danger) levels in metres on each gauge's own datum.
|
|
# "*" is the default applied to any station without an explicit override.
|
|
# Calibrated 2026-08-10 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 keyed to
|
|
# the P.1 gauge: city flooding begins at 3.70 m (stage 1) and reaches most
|
|
# districts by 4.20 m (stage 5) — see P1_FLOOD_STAGES.
|
|
THRESHOLDS: Dict[str, Tuple[float, float]] = {
|
|
"*": (3.0, 4.5),
|
|
"P.1": (3.70, 4.20),
|
|
"P.103": (5.95, 6.75),
|
|
"P.20": (2.35, 2.80),
|
|
"P.21": (3.20, 3.60),
|
|
# P.4A: LOW CONFIDENCE — its sensor was dead 2019-2024 and only 11 readings
|
|
# ever reached 3.40 m, so the capacity calibration rests on very few points.
|
|
# It only affects the heuristic sigmoid (P.4A is NOT_TRAINABLE).
|
|
"P.4A": (3.40, 3.90),
|
|
"P.5": (4.55, 4.95),
|
|
"P.67": (2.45, 2.90),
|
|
# P.75: 2024 (the only year with a full flood record, 191% capacity peak)
|
|
# puts 75-85% at 3.45 m and 95-105% at 3.72 m; 2018/2022 agree within
|
|
# 0.15 m. The 2026-08 value (2.75) alerted on 15 quiet-season hours.
|
|
"P.75": (3.20, 3.65),
|
|
"P.76": (5.35, 5.45),
|
|
# P.77: recalibrated 2026-09-12. The 2026-08 value (2.85) sat below the
|
|
# gauge's own dry-season baseline (2.6-2.7 m at 8-14% capacity), so the
|
|
# first ntfy cycle fired a "warning" at 22% capacity. Across 2018-2024,
|
|
# 75-85% capacity reads 3.35-4.57 m and 95-105% 4.27-5.08 m; 2024 (the
|
|
# best-sampled flood year) gives 4.57 / 5.08. Slightly conservative:
|
|
"P.77": (4.30, 4.90),
|
|
"P.81": (5.15, 6.30),
|
|
# P.82 never reached 100% capacity in the record (max level 3.78, max 96.4%);
|
|
# danger sits just below the observed maximum so the head can actually train.
|
|
"P.82": (3.40, 3.75),
|
|
"P.84": (3.45, 3.90),
|
|
"P.85": (2.90, 3.35),
|
|
"P.87": (3.75, 4.05),
|
|
"P.92": (2.95, 3.60),
|
|
}
|
|
|
|
# Official Chiang Mai flood-onset stages at the P.1 gauge (Nawarat Bridge),
|
|
# from the municipal inundation map (พื้นที่ท่วมตัวเมืองเชียงใหม่, events of
|
|
# 2548/2554/2565 BE): gauge level in m, RID discharge in m³/s. Each stage
|
|
# floods progressively more city zones.
|
|
P1_FLOOD_STAGES: List[Dict[str, float]] = [
|
|
{"stage": 1, "level": 3.70, "discharge_cms": 405},
|
|
{"stage": 2, "level": 3.90, "discharge_cms": 438},
|
|
{"stage": 3, "level": 4.00, "discharge_cms": 458},
|
|
{"stage": 4, "level": 4.10, "discharge_cms": 478},
|
|
{"stage": 5, "level": 4.20, "discharge_cms": 493},
|
|
{"stage": 6, "level": 4.30, "discharge_cms": 508},
|
|
{"stage": 7, "level": 4.60, "discharge_cms": 558},
|
|
]
|
|
|
|
FLOOD_STAGES: Dict[str, List[Dict[str, float]]] = {"P.1": P1_FLOOD_STAGES}
|
|
|
|
MONSOON_MONTHS = {6, 7, 8, 9, 10}
|
|
FFILL_LIMIT_H = 3
|
|
MIN_WINDOW_COVERAGE = 0.5
|
|
|
|
BASIN_ANCHOR = "P.1"
|
|
|
|
# Empirical hours a station's water-level anomaly leads the basin anchor (P.1),
|
|
# derived from data-scout cross-correlation analysis. UPSTREAM_LEADS[station]
|
|
# lists, for each station, the (upstream_code, lead_hours) pairs to use as
|
|
# routed-upstream input features when forecasting `station`.
|
|
UPSTREAM_LEADS: Dict[str, List[Tuple[str, int]]] = {
|
|
"P.1": [
|
|
("P.103", 1),
|
|
("P.67", 7),
|
|
("P.21", 9),
|
|
("P.75", 12),
|
|
("P.4A", 12),
|
|
("P.92", 15),
|
|
("P.20", 17),
|
|
],
|
|
"P.103": [
|
|
("P.67", 6),
|
|
("P.21", 8),
|
|
("P.75", 11),
|
|
("P.4A", 11),
|
|
("P.92", 14),
|
|
("P.20", 16),
|
|
],
|
|
"P.21": [("P.67", 1), ("P.75", 3), ("P.4A", 3), ("P.92", 6), ("P.20", 8)],
|
|
"P.67": [("P.75", 5), ("P.4A", 5), ("P.92", 8), ("P.20", 10)],
|
|
"P.75": [("P.92", 3), ("P.20", 5)],
|
|
"P.4A": [("P.92", 3), ("P.20", 5)],
|
|
"P.92": [("P.20", 2)],
|
|
"P.20": [],
|
|
"P.5": [("P.1", 12), ("P.103", 13)],
|
|
"P.81": [("P.1", 4), ("P.103", 5)],
|
|
"P.82": [],
|
|
"P.84": [],
|
|
"P.87": [],
|
|
"P.77": [],
|
|
"P.85": [],
|
|
"P.76": [],
|
|
}
|
|
|
|
# Per-station usable-from dates: data before this cutoff is excluded from training
|
|
# because of known data-quality holes (see data-scout inventory).
|
|
TRAIN_START: Dict[str, str] = {"P.5": "2022-01-01"}
|
|
|
|
# Stations with data too sparse/broken to ever be a regression/classification
|
|
# target. They are still usable as upstream *input* features (HGB tolerates NaN).
|
|
NOT_TRAINABLE: Dict[str, str] = {"P.4A": "17% fill, dead 2019-2024"}
|
|
|
|
|
|
def get_thresholds(station_code: str) -> Tuple[float, float]:
|
|
"""Return (warning, danger) level thresholds for a station, falling back to the default."""
|
|
return THRESHOLDS.get(station_code, THRESHOLDS["*"])
|
|
|
|
|
|
# ---------------------------------------------------------------------------
|
|
# Hourly grid
|
|
# ---------------------------------------------------------------------------
|
|
|
|
|
|
@dataclass
|
|
class HourlyGrid:
|
|
"""A complete hourly time grid pivoted wide across stations.
|
|
|
|
observed: raw values, NaN where nothing was recorded that hour (pristine; used for labels).
|
|
filled: observed forward-filled per column with limit=FFILL_LIMIT_H (causal; used for features).
|
|
mask: boolean, True where `observed` has a real reading.
|
|
"""
|
|
|
|
observed: pd.DataFrame
|
|
filled: pd.DataFrame
|
|
mask: pd.DataFrame
|
|
|
|
|
|
def make_hourly_grid(df_long: pd.DataFrame) -> HourlyGrid:
|
|
"""Pivot a long station/timestamp measurement frame onto a complete hourly grid.
|
|
|
|
df_long columns: timestamp, station_code, water_level, discharge.
|
|
"""
|
|
if df_long.empty:
|
|
empty = pd.DataFrame(
|
|
index=pd.DatetimeIndex([], name="timestamp"),
|
|
columns=pd.MultiIndex.from_tuples([], names=["station_code", "field"]),
|
|
)
|
|
return HourlyGrid(observed=empty, filled=empty.copy(), mask=empty.copy())
|
|
|
|
df = df_long.copy()
|
|
df["timestamp"] = pd.to_datetime(df["timestamp"]).dt.floor("h")
|
|
df = df.drop_duplicates(subset=["station_code", "timestamp"], keep="last")
|
|
|
|
full_index = pd.date_range(
|
|
df["timestamp"].min(), df["timestamp"].max(), freq="h", name="timestamp"
|
|
)
|
|
|
|
wide = df.pivot(
|
|
index="timestamp", columns="station_code", values=["water_level", "discharge"]
|
|
)
|
|
wide = wide.reorder_levels([1, 0], axis=1).sort_index(axis=1)
|
|
wide = wide.reindex(full_index)
|
|
|
|
observed = wide
|
|
mask = observed.notna()
|
|
# Forward-fill only — never interpolate — so no row ever depends on a future value.
|
|
filled = observed.ffill(limit=FFILL_LIMIT_H)
|
|
|
|
return HourlyGrid(observed=observed, filled=filled, mask=mask)
|
|
|
|
|
|
def _series(
|
|
grid_frame: pd.DataFrame, station: str, field: str, index: pd.Index
|
|
) -> pd.Series:
|
|
"""Fetch a (station, field) column, or an all-NaN series if the station is absent."""
|
|
if (station, field) in grid_frame.columns:
|
|
return grid_frame[(station, field)]
|
|
return pd.Series(np.nan, index=index)
|
|
|
|
|
|
def _hours_since_observed(mask_col: pd.Series) -> pd.Series:
|
|
"""Hours since the last True in `mask_col` (0 at an observed hour; NaN if never observed yet)."""
|
|
idx = mask_col.index
|
|
obs_time = pd.Series(idx, index=idx).where(mask_col.to_numpy())
|
|
last_obs_time = obs_time.ffill()
|
|
age_hours = (idx.to_series() - last_obs_time).dt.total_seconds() / 3600.0
|
|
return age_hours
|
|
|
|
|
|
# ---------------------------------------------------------------------------
|
|
# Features
|
|
# ---------------------------------------------------------------------------
|
|
|
|
|
|
RAIN_FEATURES = ("rain_6h", "rain_24h", "rain_72h", "rain_fc24")
|
|
FORWARD_RAIN_H = 24 # rain_fc24 window
|
|
|
|
DAM_FEATURES = ("dam_storage_pct", "dam_storage_pct_d3", "dam_inflow", "dam_outflow")
|
|
# Stations hydrologically downstream of the Mae Ngat confluence (Ping mainstem
|
|
# at/below Mae Taeng) — the only ones where reservoir state is causal. West-
|
|
# tributary and upper-mainstem stations never receive dam columns.
|
|
DAM_STATIONS = frozenset({"P.1", "P.103", "P.67", "P.21", "P.5", "P.81"})
|
|
|
|
|
|
def build_features(
|
|
grid: HourlyGrid,
|
|
station: str,
|
|
rain: Optional[pd.Series] = None,
|
|
dam: Optional[pd.DataFrame] = None,
|
|
) -> pd.DataFrame:
|
|
"""Build the deterministic-order feature matrix for one target station.
|
|
|
|
``rain`` is the hourly catchment-average precipitation series (Open-Meteo,
|
|
src/ml/rain.py). Rain columns are added only when a series is passed:
|
|
HistGradientBoosting REJECTS all-NaN columns at fit time, so training
|
|
without rain must omit the columns entirely (bundles record their
|
|
feature_names, and serving subsets to them). At serving, pass an empty
|
|
series rather than None so the columns exist (as NaN) for rain-trained
|
|
bundles even when the live fetch fails. rain_fc24 is the forward 24 h
|
|
sum: the archived forecast series at training time, a real weather
|
|
forecast at serving time; it never contains river data.
|
|
|
|
``dam`` is the hourly Mae Ngat reservoir frame (src/ml/dam.py; columns
|
|
storage_pct/inflow_mcm/outflow_mcm, already leakage-shifted to 07:00
|
|
report time). Same contract as rain: None omits the columns, an empty
|
|
frame yields NaN columns; only DAM_STATIONS receive them.
|
|
"""
|
|
idx = grid.observed.index
|
|
cols: Dict[str, pd.Series] = {}
|
|
|
|
level = _series(grid.filled, station, "water_level", idx)
|
|
discharge = _series(grid.observed, station, "discharge", idx)
|
|
obs_mask = _series(grid.mask, station, "water_level", idx).fillna(False)
|
|
|
|
cols["level"] = level
|
|
for k in (1, 2, 3, 6, 12, 24, 48, 72):
|
|
cols[f"level_lag_{k}"] = level.shift(k)
|
|
for k in (1, 3, 6, 12, 24):
|
|
cols[f"rise_{k}"] = level - level.shift(k)
|
|
cols["roll_mean_6"] = level.rolling(6, min_periods=1).mean()
|
|
cols["roll_mean_24"] = level.rolling(24, min_periods=1).mean()
|
|
cols["roll_max_6"] = level.rolling(6, min_periods=1).max()
|
|
cols["roll_max_24"] = level.rolling(24, min_periods=1).max()
|
|
cols["roll_max_72"] = level.rolling(72, min_periods=1).max()
|
|
cols["roll_min_24"] = level.rolling(24, min_periods=1).min()
|
|
|
|
cols["discharge"] = discharge
|
|
cols["discharge_lag_6"] = discharge.shift(6)
|
|
cols["discharge_lag_24"] = discharge.shift(24)
|
|
cols["discharge_rise_6"] = discharge - discharge.shift(6)
|
|
|
|
obs_age_h = _hours_since_observed(obs_mask)
|
|
cols["obs_age_h"] = obs_age_h.where(obs_age_h <= FFILL_LIMIT_H)
|
|
cols["cov_24h"] = obs_mask.rolling(24, min_periods=1).mean()
|
|
|
|
for upstream_code, lead_h in UPSTREAM_LEADS.get(station, []):
|
|
u_level = _series(grid.filled, upstream_code, "water_level", idx)
|
|
u_rise_6 = u_level - u_level.shift(6)
|
|
u_rollmax_24 = u_level.rolling(24, min_periods=1).max()
|
|
near_lag = max(0, lead_h - 3)
|
|
cols[f"{upstream_code}_level_lag_{near_lag}"] = u_level.shift(near_lag)
|
|
cols[f"{upstream_code}_level_lag_{lead_h}"] = u_level.shift(lead_h)
|
|
cols[f"{upstream_code}_level_lag_{lead_h + 3}"] = u_level.shift(lead_h + 3)
|
|
cols[f"{upstream_code}_rise_6_lag_{lead_h}"] = u_rise_6.shift(lead_h)
|
|
cols[f"{upstream_code}_rollmax_24_lag_{near_lag}"] = u_rollmax_24.shift(
|
|
near_lag
|
|
)
|
|
|
|
if station != BASIN_ANCHOR:
|
|
p1_level = _series(grid.filled, BASIN_ANCHOR, "water_level", idx)
|
|
cols["P1_level"] = p1_level
|
|
cols["P1_rollmax_24"] = p1_level.rolling(24, min_periods=1).max()
|
|
cols["P1_rise_24"] = p1_level - p1_level.shift(24)
|
|
|
|
doy = idx.to_series().dt.dayofyear.astype(float)
|
|
cols["doy_sin"] = np.sin(2 * np.pi * doy / 365.25)
|
|
cols["doy_cos"] = np.cos(2 * np.pi * doy / 365.25)
|
|
cols["is_monsoon"] = idx.to_series().dt.month.isin(MONSOON_MONTHS).astype(float)
|
|
|
|
if rain is not None:
|
|
r = rain.reindex(idx)
|
|
cols["rain_6h"] = r.rolling(6, min_periods=1).sum()
|
|
cols["rain_24h"] = r.rolling(24, min_periods=1).sum()
|
|
cols["rain_72h"] = r.rolling(72, min_periods=1).sum()
|
|
# forward sum over (t, t+24]: shift(-1) starts the window at t+1.
|
|
# Computed on the grid extended by the 24 h after its last row: the
|
|
# last row is the one serving predicts from, and cutting the rain
|
|
# series to the grid made its forecast rain always NaN (unlike every
|
|
# training row since 2021). Earlier rows keep their values.
|
|
future = (
|
|
pd.date_range(idx[-1], periods=FORWARD_RAIN_H + 1, freq="h")[1:]
|
|
if len(idx)
|
|
else idx[:0]
|
|
)
|
|
r_ext = rain.reindex(idx.append(future))
|
|
cols["rain_fc24"] = (
|
|
r_ext.shift(-1)
|
|
.iloc[::-1]
|
|
.rolling(FORWARD_RAIN_H, min_periods=1)
|
|
.sum()
|
|
.iloc[::-1]
|
|
.iloc[: len(idx)]
|
|
.set_axis(idx)
|
|
)
|
|
|
|
if dam is not None and station in DAM_STATIONS:
|
|
d = dam.reindex(idx)
|
|
|
|
def _dam_col(name: str) -> pd.Series:
|
|
return d[name] if name in d.columns else pd.Series(np.nan, index=idx)
|
|
|
|
storage = _dam_col("storage_pct")
|
|
cols["dam_storage_pct"] = storage
|
|
cols["dam_storage_pct_d3"] = storage - storage.shift(72)
|
|
cols["dam_inflow"] = _dam_col("inflow_mcm")
|
|
cols["dam_outflow"] = _dam_col("outflow_mcm")
|
|
|
|
return pd.DataFrame(cols, index=idx)
|
|
|
|
|
|
# ---------------------------------------------------------------------------
|
|
# Labels
|
|
# ---------------------------------------------------------------------------
|
|
|
|
|
|
def _future_window_stats(col: pd.Series, horizon_h: int) -> Tuple[pd.Series, pd.Series]:
|
|
"""For every t, (max, count) of observed values in the OPEN window (t, t+horizon_h]."""
|
|
reversed_col = col.iloc[::-1]
|
|
shifted = reversed_col.shift(1) # excludes t itself
|
|
fut_max = shifted.rolling(horizon_h, min_periods=1).max().iloc[::-1]
|
|
fut_count = shifted.rolling(horizon_h, min_periods=1).count().iloc[::-1]
|
|
return fut_max, fut_count
|
|
|
|
|
|
def build_labels(
|
|
grid: HourlyGrid,
|
|
station: str,
|
|
horizons: Tuple[int, ...] = (6, 12, 24),
|
|
stats_end: Optional[str] = None,
|
|
) -> pd.DataFrame:
|
|
"""Build max-level and threshold-exceedance labels for one target station.
|
|
|
|
``stats_end`` bounds the data used for label-construction statistics (the
|
|
rescue quantile below): pass the training cutoff during evaluation so
|
|
test-period extremes cannot influence which training rows receive labels.
|
|
"""
|
|
idx = grid.observed.index
|
|
observed_level = _series(grid.observed, station, "water_level", idx)
|
|
warn_thr, danger_thr = get_thresholds(station)
|
|
# Low-coverage windows are still usable regression labels when they contain a
|
|
# rare high reading. Anchor "rare" to the station's own distribution (p97.5),
|
|
# NOT to warn_thr: coupling it to the configurable threshold made raising a
|
|
# station's threshold silently shrink its regression training set (P.5 lost
|
|
# 34% of rows and +46% MAE when its warning went 3.0 -> 4.55).
|
|
stats_level = (
|
|
observed_level.loc[: pd.Timestamp(stats_end)] if stats_end else observed_level
|
|
)
|
|
rescue_thr = (
|
|
float(stats_level.quantile(0.975)) if stats_level.notna().any() else np.inf
|
|
)
|
|
|
|
out: Dict[str, pd.Series] = {}
|
|
for horizon_h in horizons:
|
|
fut_max, fut_count = _future_window_stats(observed_level, horizon_h)
|
|
cov = fut_count / horizon_h
|
|
enough_cov = cov >= MIN_WINDOW_COVERAGE
|
|
|
|
exceed_warn = pd.Series(np.nan, index=idx)
|
|
exceed_warn[fut_max >= warn_thr] = 1.0
|
|
exceed_warn[enough_cov & exceed_warn.isna()] = 0.0
|
|
|
|
exceed_danger = pd.Series(np.nan, index=idx)
|
|
exceed_danger[fut_max >= danger_thr] = 1.0
|
|
exceed_danger[enough_cov & exceed_danger.isna()] = 0.0
|
|
|
|
max_level_valid = fut_max.where(enough_cov | (fut_max >= rescue_thr))
|
|
|
|
out[f"max_level_{horizon_h}"] = max_level_valid
|
|
out[f"exceed_warn_{horizon_h}"] = exceed_warn
|
|
out[f"exceed_danger_{horizon_h}"] = exceed_danger
|
|
|
|
return pd.DataFrame(out, index=idx)
|
|
|
|
|
|
# ---------------------------------------------------------------------------
|
|
# Glue
|
|
# ---------------------------------------------------------------------------
|
|
|
|
|
|
def build_matrix(
|
|
df_long: pd.DataFrame,
|
|
station: str,
|
|
horizons: Tuple[int, ...] = (6, 12, 24),
|
|
stats_end: Optional[str] = None,
|
|
rain: Optional[pd.Series] = None,
|
|
dam: Optional[pd.DataFrame] = None,
|
|
) -> Tuple[pd.DataFrame, pd.DataFrame, dict]:
|
|
"""Build (X, Y, meta) training/inference matrices for one station."""
|
|
grid = make_hourly_grid(df_long)
|
|
X = build_features(grid, station, rain=rain, dam=dam)
|
|
Y = build_labels(grid, station, horizons, stats_end=stats_end)
|
|
|
|
keep = X["obs_age_h"].notna()
|
|
train_start = TRAIN_START.get(station)
|
|
if train_start:
|
|
keep &= X.index >= pd.Timestamp(train_start)
|
|
|
|
X = X.loc[keep]
|
|
Y = Y.loc[keep]
|
|
|
|
positive_counts = {
|
|
col: int(Y[col].sum())
|
|
for col in Y.columns
|
|
if col.startswith("exceed_") and Y[col].notna().any()
|
|
}
|
|
meta = {
|
|
"station_code": station,
|
|
"n_rows": int(len(X)),
|
|
"span": (
|
|
(X.index.min().isoformat(), X.index.max().isoformat())
|
|
if len(X)
|
|
else (None, None)
|
|
),
|
|
"positive_counts": positive_counts,
|
|
}
|
|
return X, Y, meta
|