B4L/Northern-Thailand-Ping-River-Monitor
5
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 15 Wiki Activity

feat: ML flood-event forecasting from 8 years of gauge history
Security & Dependency Updates / Dependency Security Scan (push) Successful in 1m8s
Details
Security & Dependency Updates / License Compliance (push) Successful in 25s
Details
Security & Dependency Updates / Check for Dependency Updates (push) Successful in 20s
Details
Security & Dependency Updates / Code Quality Metrics (push) Successful in 17s
Details
Security & Dependency Updates / Security Summary (push) Successful in 9s
Details

Browse Source 
Add src/ml/ package predicting, per station and per 6/12/24 h horizon,
the probability of exceeding warning (3.0 m) and danger (4.5 m) levels
plus expected peak level, trained on the 592k-row PostgreSQL history:

- features.py: hourly grid with coverage gating and no future leakage;
  upstream stations enter at empirically measured travel-time lags
  (P.20 +17h ... P.103 +1h vs P.1); hour-of-day deliberately excluded
  (it encodes the scrape schedule, not hydrology)
- train.py: HistGradientBoosting regression + warn/danger classifier
  heads per station x horizon, >=30-positives gate with calibrated
  sigmoid-on-regression fallback, strict temporal splits, per-event
  lead-time evaluation; guards against sklearn 1.9.0 crash on
  degenerate feature columns
- predict.py: bundle loading with feature-name checks, heuristic
  fallback tier, get_latest_forecasts() for the API; raises when no
  models are trained so the endpoint 503s instead of serving
  persistence output as forecasts
- data.py: Postgres-first loader (FLOOD_ML_DB_URL override), HTTP API
  fallback (flagged: that path backfills synthetic discharge), csv.gz
  cache
- /forecast endpoint (15-min TTL cache) + dashboard flood-risk panel
  (hidden until models exist)
- docs/FLOOD_FORECASTING.md: full system doc with measured deployment
  numbers (~335 MB RSS, CPU negligible, ~6 min full retrain) and
  retraining policy

Validation: out-of-sample backtest of the record 2024 flood season
(train <= Aug 2024) alerted 24-48 h ahead of the Oct 5 peak; 2025-26
test split: P.1 6h PR-AUC 0.974, recall 98.3% at 1% false-alarm rate.

Also: fix P.81 station coordinates (was Ban Pong/Ratchaburi, 493 km
out of basin; now 18.6936 N 99.0819 E per RID station page), pin
scikit-learn==1.9.0 and numpy<2, gitignore model artifacts (~100 MB,
train on the server via scripts/train_flood_model.py).
This commit is contained in:
grabowski
2026-08-10 12:49:47 +07:00
parent 49a3de0087
commit 4358d52d55
15 changed files with 2188 additions and 2 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/ml/__init__.py
+1
View File
@@ -0,0 +1 @@
"""Flood forecasting ML package: data loading, feature/label engineering, training, and prediction."""
src/ml/data.py
+223
View File
@@ -0,0 +1,223 @@
"""Loaders for flood-model training/inference data.
Primary path reads raw measurements straight from PostgreSQL (keeping NULL
discharge as NULL). HTTP fallback goes through the public API's history
endpoint, which backfills missing discharge with a synthetic rating-curve
estimate -- callers are told about that via the `discharge_maybe_synthetic`
cache metadata flag.
"""
import datetime
import gzip
import json
import logging
import os
from pathlib import Path
from typing import Dict, List, Optional
import pandas as pd
from sqlalchemy import create_engine, text
from ..config import Config
from .features import UPSTREAM_LEADS
logger = logging.getLogger(__name__)
DEFAULT_API_URL = "http://100.81.167.42:8000"
CACHE_DIR = Path("models/cache")
_MEASUREMENT_COLUMNS = ["timestamp", "station_code", "water_level", "discharge"]
def resolve_db_url(db_url: Optional[str] = None) -> Optional[str]:
"""Resolve a Postgres connection string: explicit param > FLOOD_ML_DB_URL env >
Config's postgresql connection string > None (caller should fall back to HTTP)."""
if db_url:
return db_url
env_url = os.getenv("FLOOD_ML_DB_URL")
if env_url:
return env_url
try:
db_config = Config.get_database_config()
except Exception as error:
logger.warning(f"Could not resolve database config: {error}")
return None
if db_config.get("type") == "postgresql":
return db_config.get("connection_string")
return None
def _default_stations() -> List[str]:
return list(UPSTREAM_LEADS.keys())
def _normalize_long(df: pd.DataFrame) -> pd.DataFrame:
if df.empty:
return pd.DataFrame(columns=_MEASUREMENT_COLUMNS)
df = df.copy()
df["timestamp"] = pd.to_datetime(df["timestamp"]).dt.floor("h")
df["water_level"] = pd.to_numeric(df["water_level"], errors="coerce")
df["discharge"] = pd.to_numeric(df["discharge"], errors="coerce")
df = df.drop_duplicates(subset=["station_code", "timestamp"], keep="last")
df = df.sort_values("timestamp").reset_index(drop=True)
return df[_MEASUREMENT_COLUMNS]
def _fetch_from_db(
db_url: str,
stations: Optional[List[str]],
start: Optional[datetime.datetime],
end: Optional[datetime.datetime],
) -> pd.DataFrame:
engine = create_engine(db_url, pool_pre_ping=True)
query = (
"SELECT m.timestamp, s.station_code, m.water_level, m.discharge "
"FROM water_measurements m JOIN stations s ON m.station_id = s.id WHERE 1=1"
)
params: Dict = {}
if start is not None:
query += " AND m.timestamp >= :start_time"
params["start_time"] = start
if end is not None:
query += " AND m.timestamp <= :end_time"
params["end_time"] = end
if stations:
placeholders = ", ".join(f":station_{i}" for i in range(len(stations)))
query += f" AND s.station_code IN ({placeholders})"
for i, code in enumerate(stations):
params[f"station_{i}"] = code
query += " ORDER BY m.timestamp"
with engine.connect() as connection:
df = pd.read_sql(text(query), connection, params=params)
return _normalize_long(df)
def _fetch_station_from_api(api_url: str, station_code: str, hours: int, limit: int = 100000) -> pd.DataFrame:
import requests
response = requests.get(
f"{api_url}/measurements/history/{station_code}",
params={"hours": hours, "limit": limit},
timeout=30,
)
response.raise_for_status()
rows = response.json()
for row in rows:
row["station_code"] = station_code
return pd.DataFrame(rows, columns=_MEASUREMENT_COLUMNS + ["discharge_percent"])
def _fetch_from_api(
api_url: str,
stations: List[str],
start: Optional[datetime.datetime],
end: Optional[datetime.datetime],
) -> pd.DataFrame:
now = datetime.datetime.now()
reference_end = end or now
reference_start = start or (reference_end - datetime.timedelta(days=365 * 8))
hours = max(1, int((reference_end - reference_start).total_seconds() // 3600) + 1)
frames = []
for code in stations:
try:
frames.append(_fetch_station_from_api(api_url, code, hours))
except Exception as error:
logger.warning(f"HTTP fallback failed for station {code}: {error}")
if not frames:
return pd.DataFrame(columns=_MEASUREMENT_COLUMNS)
df = pd.concat(frames, ignore_index=True)
return _normalize_long(df)
def _write_cache(df: pd.DataFrame, cache_dir: Path, source: str, discharge_maybe_synthetic: bool) -> None:
cache_dir.mkdir(parents=True, exist_ok=True)
for code, group in df.groupby("station_code"):
path = cache_dir / f"{code}.csv.gz"
with gzip.open(path, "wt", encoding="utf-8", newline="") as handle:
group.to_csv(handle, index=False)
meta = {
"fetched_at": datetime.datetime.now().isoformat(),
"source": source,
"discharge_maybe_synthetic": discharge_maybe_synthetic,
}
with open(cache_dir / "meta.json", "w", encoding="utf-8") as handle:
json.dump(meta, handle)
def _read_cache(cache_dir: Path, stations: Optional[List[str]]) -> pd.DataFrame:
if not cache_dir.exists():
return pd.DataFrame(columns=_MEASUREMENT_COLUMNS)
frames = []
for path in sorted(cache_dir.glob("*.csv.gz")):
code = path.name[: -len(".csv.gz")]
if stations and code not in stations:
continue
with gzip.open(path, "rt", encoding="utf-8") as handle:
frames.append(pd.read_csv(handle, parse_dates=["timestamp"]))
if not frames:
return pd.DataFrame(columns=_MEASUREMENT_COLUMNS)
return _normalize_long(pd.concat(frames, ignore_index=True))
def load_measurements(
db_url: Optional[str] = None,
stations: Optional[List[str]] = None,
start: Optional[datetime.datetime] = None,
end: Optional[datetime.datetime] = None,
use_cache: bool = True,
cache_dir: Path = CACHE_DIR,
api_url: str = DEFAULT_API_URL,
) -> pd.DataFrame:
"""Load the long-format [timestamp, station_code, water_level, discharge] history.
Tries PostgreSQL first, then the HTTP API, then the on-disk cache as a last
resort. A successful DB/API fetch refreshes the cache; the cache itself is
never treated as a source of fresh data.
"""
resolved_db_url = resolve_db_url(db_url)
if resolved_db_url:
try:
df = _fetch_from_db(resolved_db_url, stations, start, end)
if use_cache:
_write_cache(df, cache_dir, source="postgres", discharge_maybe_synthetic=False)
return df
except Exception as error:
logger.warning(f"PostgreSQL fetch failed, falling back to HTTP API: {error}")
try:
api_stations = stations or _default_stations()
df = _fetch_from_api(api_url, api_stations, start, end)
if not df.empty:
if use_cache:
_write_cache(df, cache_dir, source="api", discharge_maybe_synthetic=True)
return df
except Exception as error:
logger.warning(f"HTTP API fetch failed: {error}")
if use_cache:
logger.warning("Falling back to on-disk cache for measurement history")
return _read_cache(cache_dir, stations)
return pd.DataFrame(columns=_MEASUREMENT_COLUMNS)
def load_latest(
db_url: Optional[str] = None,
hours: int = 336,
stations: Optional[List[str]] = None,
) -> pd.DataFrame:
"""Load the last `hours` of history for all (or given) stations. Never cached to disk."""
end = datetime.datetime.now()
start = end - datetime.timedelta(hours=hours)
return load_measurements(
db_url=db_url,
stations=stations,
start=start,
end=end,
use_cache=False,
)
src/ml/features.py
+272
View File
@@ -0,0 +1,272 @@
"""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) level thresholds in meters. "*" is the default
# applied to any station without an explicit override.
THRESHOLDS: Dict[str, Tuple[float, float]] = {
"*": (3.0, 4.5),
}
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
# ---------------------------------------------------------------------------
def build_features(grid: HourlyGrid, station: str) -> pd.DataFrame:
"""Build the deterministic-order feature matrix for one target station."""
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)
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)) -> pd.DataFrame:
"""Build max-level and threshold-exceedance labels for one target station."""
idx = grid.observed.index
observed_level = _series(grid.observed, station, "water_level", idx)
warn_thr, danger_thr = get_thresholds(station)
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 >= warn_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),
) -> 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)
Y = build_labels(grid, station, horizons)
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
src/ml/predict.py
+284
View File
@@ -0,0 +1,284 @@
"""Flood forecast inference.
Integration contract (see get_forecasts / get_latest_forecasts): callers pass
raw station readings, get back one forecast dict per station x horizon. A
station with a stale, missing, or version-mismatched model transparently
falls back to a simple persistence heuristic instead of raising -- this
module must never crash the caller (e.g. the web API).
"""
import datetime
import logging
from pathlib import Path
from typing import Dict, List, Optional, Tuple, Union
import joblib
import numpy as np
import pandas as pd
from . import features
logger = logging.getLogger(__name__)
DEFAULT_HORIZONS: Tuple[int, ...] = (6, 12, 24)
STALE_AFTER_H = 6.0
HEURISTIC_SIGMA = 0.3
HEURISTIC_VERSION = "heuristic-v1"
# Keyed by (path, mtime) so a retrained model (new mtime) invalidates the old entry.
_MODEL_CACHE: Dict[Tuple[str, float], dict] = {}
def _load_bundle(path: Path) -> dict:
# joblib.load runs arbitrary pickle code; safe here because `path` is always
# models/flood_{station}.joblib, an artifact this pipeline's own train.py wrote --
# never a user- or network-supplied file.
key = (str(path), path.stat().st_mtime)
cached = _MODEL_CACHE.get(key)
if cached is not None:
return cached
bundle = joblib.load(path)
for stale_key in [k for k in _MODEL_CACHE if k[0] == str(path)]:
del _MODEL_CACHE[stale_key]
_MODEL_CACHE[key] = bundle
return bundle
def _readings_to_long_df(readings_by_station: Dict[str, List[dict]]) -> pd.DataFrame:
rows = []
for station_code, readings in readings_by_station.items():
for reading in readings:
timestamp = reading.get("timestamp")
if isinstance(timestamp, str):
timestamp = pd.to_datetime(timestamp)
rows.append(
{
"timestamp": timestamp,
"station_code": station_code,
"water_level": reading.get("water_level"),
"discharge": reading.get("discharge"),
}
)
if not rows:
return pd.DataFrame(columns=["timestamp", "station_code", "water_level", "discharge"])
df = pd.DataFrame(rows)
return df.dropna(subset=["timestamp"])
def _clip_probability(value: float) -> float:
return float(min(max(value, 0.0), 1.0))
def _sigmoid_probability(predicted_max: float, threshold: float, sigma: float) -> float:
return 1.0 / (1.0 + np.exp(-(predicted_max - threshold) / sigma))
def _heuristic_forecast(
station_code: str,
as_of: pd.Timestamp,
current_level: float,
level_t_minus_3: Optional[float],
warn_thr: float,
danger_thr: float,
horizons: Tuple[int, ...],
) -> List[dict]:
if level_t_minus_3 is None:
rate = 0.0
else:
rate = max(0.0, (current_level - level_t_minus_3) / 3.0)
results = []
for horizon_h in horizons:
predicted_max = max(current_level + rate * horizon_h * 0.7, current_level)
p_warning = _clip_probability(_sigmoid_probability(predicted_max, warn_thr, HEURISTIC_SIGMA))
p_danger = _clip_probability(_sigmoid_probability(predicted_max, danger_thr, HEURISTIC_SIGMA))
p_danger = min(p_danger, p_warning)
results.append(
{
"station_code": station_code,
"horizon_hours": horizon_h,
"p_warning": p_warning,
"p_danger": p_danger,
"predicted_max_level": predicted_max,
"current_level": current_level,
"as_of": as_of.isoformat(),
"model_version": HEURISTIC_VERSION,
"trained_at": None,
"source": "heuristic",
"threshold_warning": warn_thr,
"threshold_danger": danger_thr,
}
)
return results
def _model_forecast(
station_code: str,
grid: features.HourlyGrid,
bundle: dict,
as_of: pd.Timestamp,
current_level: float,
) -> List[dict]:
warn_thr = bundle["thresholds"]["warning"]
danger_thr = bundle["thresholds"]["danger"]
feature_row = features.build_features(grid, station_code).loc[[as_of]]
expected_columns = bundle["feature_names"]
missing = [c for c in expected_columns if c not in feature_row.columns]
if missing:
logger.error(f"Feature mismatch for {station_code} (missing {missing}); falling back to heuristic")
return None
feature_row = feature_row[expected_columns]
results = []
for horizon_h in bundle["horizons"]:
reg = bundle["heads"].get(f"max_{horizon_h}")
if reg is None:
results.append(None)
continue
predicted_max = max(float(reg.predict(feature_row)[0]), current_level)
sigma_h = bundle["sigma"].get(horizon_h, HEURISTIC_SIGMA)
warn_head = bundle["heads"].get(f"warn_{horizon_h}")
if warn_head is not None:
p_warning = float(warn_head.predict_proba(feature_row)[0][1])
else:
p_warning = _sigmoid_probability(predicted_max, warn_thr, sigma_h)
danger_head = bundle["heads"].get(f"danger_{horizon_h}")
if danger_head is not None:
p_danger = float(danger_head.predict_proba(feature_row)[0][1])
else:
p_danger = _sigmoid_probability(predicted_max, danger_thr, sigma_h)
p_warning = _clip_probability(p_warning)
p_danger = min(_clip_probability(p_danger), p_warning)
results.append(
{
"station_code": station_code,
"horizon_hours": horizon_h,
"p_warning": p_warning,
"p_danger": p_danger,
"predicted_max_level": predicted_max,
"current_level": current_level,
"as_of": as_of.isoformat(),
"model_version": bundle["model_version"],
"trained_at": bundle["trained_at"],
"source": "model",
"threshold_warning": warn_thr,
"threshold_danger": danger_thr,
}
)
return results
def _forecast_station(
station_code: str,
grid: features.HourlyGrid,
models_dir: Path,
now: pd.Timestamp,
horizons: Tuple[int, ...],
) -> List[dict]:
level_col = (station_code, "water_level")
if level_col not in grid.observed.columns:
logger.warning(f"No data for station {station_code}; omitting")
return []
observed_level = grid.observed[level_col].dropna()
if observed_level.empty:
logger.warning(f"No observed readings for station {station_code}; omitting")
return []
as_of = observed_level.index.max()
current_level = float(observed_level.loc[as_of])
staleness_h = (pd.Timestamp(now) - as_of).total_seconds() / 3600.0
warn_thr, danger_thr = features.get_thresholds(station_code)
t_minus_3 = as_of - pd.Timedelta(hours=3)
level_t_minus_3 = float(observed_level.loc[t_minus_3]) if t_minus_3 in observed_level.index else None
bundle_path = models_dir / f"flood_{station_code}.joblib"
if not bundle_path.exists() or staleness_h > STALE_AFTER_H:
return _heuristic_forecast(
station_code, as_of, current_level, level_t_minus_3, warn_thr, danger_thr, horizons
)
bundle = _load_bundle(bundle_path)
model_results = _model_forecast(station_code, grid, bundle, as_of, current_level)
if model_results is None:
return _heuristic_forecast(
station_code, as_of, current_level, level_t_minus_3, warn_thr, danger_thr, horizons
)
# Per-horizon heads that were skipped at train time (e.g. too few positives) still
# need a forecast row -- fall back to the single-horizon heuristic for just that row.
filled = []
for horizon_h, row in zip(bundle["horizons"], model_results):
if row is not None:
filled.append(row)
else:
filled.extend(
_heuristic_forecast(
station_code, as_of, current_level, level_t_minus_3, warn_thr, danger_thr, (horizon_h,)
)
)
return filled
def get_forecasts(
readings_by_station: Dict[str, List[dict]],
models_dir: Union[str, Path] = "models",
now: Optional[Union[datetime.datetime, str]] = None,
) -> List[dict]:
"""Produce flood forecasts for every station present in `readings_by_station`.
Each reading dict needs at least {timestamp, water_level, discharge}; extra
keys are ignored so raw API/DB rows can be passed straight through. At
least 96 hours of span is required to populate every feature; 336 hours
(14 days) is recommended.
"""
models_dir = Path(models_dir)
if now is None:
now = datetime.datetime.now()
now = pd.Timestamp(now)
df_long = _readings_to_long_df(readings_by_station)
if df_long.empty:
return []
grid = features.make_hourly_grid(df_long)
results: List[dict] = []
for station_code in readings_by_station.keys():
try:
results.extend(_forecast_station(station_code, grid, models_dir, now, DEFAULT_HORIZONS))
except Exception as error:
logger.error(f"Forecast failed for station {station_code}: {error}")
return results
def get_latest_forecasts(
db_url: Optional[str] = None,
models_dir: Union[str, Path] = "models",
hours: int = 336,
) -> List[dict]:
"""Convenience wrapper for web_api: load the latest window from the DB/API and forecast.
Raises FileNotFoundError when no trained model bundle exists at all, so the
API can 503 instead of serving purely heuristic output as if it were a forecast.
"""
from .data import load_latest
if not sorted(Path(models_dir).glob("flood_*.joblib")):
raise FileNotFoundError(f"no trained model bundles in {models_dir}")
df_long = load_latest(db_url=db_url, hours=hours)
readings_by_station: Dict[str, List[dict]] = {}
if not df_long.empty:
for station_code, group in df_long.groupby("station_code"):
readings_by_station[station_code] = group[["timestamp", "water_level", "discharge"]].to_dict("records")
expected_stations = set(features.UPSTREAM_LEADS.keys())
for missing_station in expected_stations - set(readings_by_station.keys()):
logger.warning(f"No recent data for station {missing_station}; omitting from forecasts")
return get_forecasts(readings_by_station, models_dir=models_dir)
src/ml/train.py
+414
View File
@@ -0,0 +1,414 @@
"""Training CLI for the Ping River flood forecast models.
Per station: build the feature/label matrix once, evaluate with a strict
temporal holdout (Split B), then refit each head on the full record for the
deployed artifact. Hyperparameters are fixed (chosen via an earlier Split A
sweep, not repeated here) -- no random search, no shuffling, no sklearn
early_stopping (its internal validation split is random and would leak
across time).
"""
import argparse
import datetime
import json
import logging
import subprocess
from pathlib import Path
from typing import Dict, List, Optional, Tuple
import joblib
import numpy as np
import pandas as pd
import sklearn
from sklearn.ensemble import HistGradientBoostingClassifier, HistGradientBoostingRegressor
from sklearn.metrics import average_precision_score, brier_score_loss, mean_absolute_error, mean_squared_error
from . import features
from .data import DEFAULT_API_URL, load_measurements, resolve_db_url
logger = logging.getLogger(__name__)
HORIZONS: Tuple[int, ...] = (6, 12, 24)
SPLIT_B_TRAIN_END = "2024-12-31"
SPLIT_B_TEST_START = "2025-01-01"
SPLIT_B_TEST_END = "2026-08-10"
MIN_POSITIVES_FOR_CLASSIFIER = 30
MIN_SIGMA = 0.15
MIN_ROWS_TO_TRAIN = 200
MIN_ROWS_FOR_HEAD = 50
HGB_PARAMS = {
"max_iter": 300,
"learning_rate": 0.06,
"max_leaf_nodes": 31,
"min_samples_leaf": 50,
"l2_regularization": 1.0,
"early_stopping": False,
"random_state": 42,
}
def _git_short_sha() -> str:
try:
result = subprocess.run(
["git", "rev-parse", "--short", "HEAD"], capture_output=True, text=True, timeout=5, check=True
)
sha = result.stdout.strip()
return sha or "nogit"
except Exception:
return "nogit"
def _make_regressor(overrides: Optional[dict] = None) -> HistGradientBoostingRegressor:
params = {**HGB_PARAMS, **(overrides or {})}
return HistGradientBoostingRegressor(loss="squared_error", **params)
def _make_classifier(overrides: Optional[dict] = None) -> HistGradientBoostingClassifier:
params = {**HGB_PARAMS, **(overrides or {})}
return HistGradientBoostingClassifier(**params)
def _safe_fit(estimator, X: pd.DataFrame, y: pd.Series, head_key: str, skipped_heads: Dict[str, str]):
"""Fit an estimator, converting any failure (e.g. HistGradientBoosting's binning
step rejecting an all-NaN/constant feature column) into a recorded skip rather
than a station-killing exception."""
try:
estimator.fit(X, y)
return estimator
except Exception as error:
skipped_heads[head_key] = f"fit failed: {error}"
logger.warning(f"{head_key}: fit failed, skipping ({error})")
return None
def _recall_at_far(y_true: np.ndarray, y_score: np.ndarray, target_far: float) -> Optional[float]:
"""Recall at the score threshold whose false-positive rate over true negatives is <= target_far."""
y_true = np.asarray(y_true)
y_score = np.asarray(y_score)
neg_scores = np.sort(y_score[y_true == 0])[::-1]
n_pos = int((y_true == 1).sum())
n_neg = len(neg_scores)
if n_pos == 0 or n_neg == 0:
return None
k = int(np.floor(target_far * n_neg))
threshold = neg_scores[k - 1] if k > 0 else neg_scores[0] + 1e-9
predicted_positive = y_score >= threshold
tp = int(np.sum(predicted_positive & (y_true == 1)))
return tp / n_pos
def _p_warning_series(head, reg, X: pd.DataFrame, threshold: float, sigma: float) -> pd.Series:
"""Model score if a classifier head exists, else the sigmoid-derived fallback probability."""
if head is not None:
return pd.Series(head.predict_proba(X)[:, 1], index=X.index)
predicted_max = pd.Series(reg.predict(X), index=X.index)
return 1.0 / (1.0 + np.exp(-(predicted_max - threshold) / sigma))
def _find_events(observed_level: pd.Series, warn_thr: float) -> List[dict]:
"""Group contiguous observed hours >= warn_thr into flood events."""
above = observed_level >= warn_thr
events: List[dict] = []
start = None
prev_t = None
for t, is_above in above.items():
if is_above and start is None:
start = t
elif not is_above and start is not None:
window = observed_level.loc[start:prev_t]
events.append({"crossed_warn_at": start, "peak_time": window.idxmax(), "peak_level": float(window.max())})
start = None
prev_t = t
if start is not None:
window = observed_level.loc[start:]
events.append({"crossed_warn_at": start, "peak_time": window.idxmax(), "peak_level": float(window.max())})
return events
def _first_alert_at(p_series: pd.Series, crossed_at, lookback_h: int = 48):
"""Earliest time p_warning was sustained (>=0.5 for 2 consecutive hours) within the prior lookback_h."""
window = p_series.loc[crossed_at - pd.Timedelta(hours=lookback_h) : crossed_at]
sustained = (window >= 0.5) & (window.shift(1) >= 0.5)
hits = sustained[sustained].index
if len(hits) == 0:
return None
return hits.min() - pd.Timedelta(hours=1)
def _events_with_lead_time(
observed_level_test: pd.Series, warn_thr: float, p_warning_test: pd.Series
) -> List[dict]:
events = _find_events(observed_level_test, warn_thr)
for event in events:
first_alert_at = _first_alert_at(p_warning_test, event["crossed_warn_at"])
event["first_alert_at"] = first_alert_at.isoformat() if first_alert_at is not None else None
if first_alert_at is not None:
lead_hours = (event["crossed_warn_at"] - first_alert_at).total_seconds() / 3600.0
else:
lead_hours = None
event["lead_hours"] = lead_hours
event["crossed_warn_at"] = event["crossed_warn_at"].isoformat()
event["peak_time"] = event["peak_time"].isoformat()
return events
def train_station(
df_long: pd.DataFrame,
station: str,
horizons: Tuple[int, ...] = HORIZONS,
skip_eval: bool = False,
hgb_overrides: Optional[dict] = None,
split_train_end: str = SPLIT_B_TRAIN_END,
split_test_start: str = SPLIT_B_TEST_START,
split_test_end: str = SPLIT_B_TEST_END,
) -> Tuple[Optional[dict], dict]:
"""Train every head for one station. Returns (bundle_or_None, station_metrics)."""
X, Y, meta = features.build_matrix(df_long, station, horizons)
if meta["n_rows"] < MIN_ROWS_TO_TRAIN:
return None, {"status": "failed", "reason": f"only {meta['n_rows']} usable rows (< {MIN_ROWS_TO_TRAIN})"}
warn_thr, danger_thr = features.get_thresholds(station)
feature_names = list(X.columns)
if skip_eval:
train_mask = pd.Series(True, index=X.index)
test_mask = pd.Series(False, index=X.index)
else:
train_mask = X.index <= pd.Timestamp(split_train_end)
test_mask = (X.index >= pd.Timestamp(split_test_start)) & (X.index <= pd.Timestamp(split_test_end))
X_train, Y_train = X.loc[train_mask], Y.loc[train_mask]
X_test, Y_test = X.loc[test_mask], Y.loc[test_mask]
eval_X, eval_Y = (X, Y) if skip_eval else (X_train, Y_train)
heads: Dict[str, object] = {}
sigma: Dict[int, float] = {}
skipped_heads: Dict[str, str] = {}
per_horizon: Dict[int, dict] = {}
observed_grid = features.make_hourly_grid(df_long).observed
for h in horizons:
max_col, warn_col, danger_col = f"max_level_{h}", f"exceed_warn_{h}", f"exceed_danger_{h}"
horizon_metrics: dict = {}
# --- regression head (max level) ---
reg_labeled = eval_Y[max_col].notna()
reg = None
if reg_labeled.sum() >= MIN_ROWS_FOR_HEAD:
reg = _safe_fit(
_make_regressor(hgb_overrides),
eval_X.loc[reg_labeled],
eval_Y.loc[reg_labeled, max_col],
f"max_{h}",
skipped_heads,
)
else:
skipped_heads[f"max_{h}"] = f"only {int(reg_labeled.sum())} labeled rows"
sigma_h = MIN_SIGMA
if reg is not None and not skip_eval:
test_labeled = Y_test[max_col].notna()
if test_labeled.sum() > 0:
y_true = Y_test.loc[test_labeled, max_col]
y_pred = reg.predict(X_test.loc[test_labeled])
residuals = y_true.to_numpy() - y_pred
sigma_h = max(float(np.std(residuals)), MIN_SIGMA)
horizon_metrics["n_test"] = int(test_labeled.sum())
horizon_metrics["mae"] = float(mean_absolute_error(y_true, y_pred))
horizon_metrics["rmse"] = float(np.sqrt(mean_squared_error(y_true, y_pred)))
above_2m = y_true >= 2.0
horizon_metrics["mae_above_2m"] = (
float(mean_absolute_error(y_true[above_2m], y_pred[above_2m])) if above_2m.any() else None
)
sigma[h] = sigma_h
horizon_metrics["sigma"] = sigma_h
# --- classification heads (warn / danger) ---
p_warning_test = None
for label_name, col, thr in (("warn", warn_col, warn_thr), ("danger", danger_col, danger_thr)):
train_labeled = eval_Y[col].notna()
n_pos = int(eval_Y.loc[train_labeled, col].sum()) if train_labeled.any() else 0
head_key = f"{label_name}_{h}"
clf = None
if n_pos >= MIN_POSITIVES_FOR_CLASSIFIER:
clf = _safe_fit(
_make_classifier(hgb_overrides),
eval_X.loc[train_labeled],
eval_Y.loc[train_labeled, col],
head_key,
skipped_heads,
)
else:
skipped_heads[head_key] = f"only {n_pos} positives in train span (< {MIN_POSITIVES_FOR_CLASSIFIER})"
heads[head_key] = clf
if not skip_eval:
test_labeled = Y_test[col].notna()
horizon_metrics[f"base_rate_{label_name}"] = (
float(Y_test.loc[test_labeled, col].mean()) if test_labeled.any() else None
)
if clf is not None and test_labeled.sum() > 0 and Y_test.loc[test_labeled, col].nunique() > 1:
y_true = Y_test.loc[test_labeled, col]
y_score = clf.predict_proba(X_test.loc[test_labeled])[:, 1]
horizon_metrics[f"pr_auc_{label_name}"] = float(average_precision_score(y_true, y_score))
horizon_metrics[f"brier_{label_name}"] = float(brier_score_loss(y_true, y_score))
horizon_metrics[f"recall_{label_name}_at_far1pct"] = _recall_at_far(y_true, y_score, 0.01)
horizon_metrics[f"recall_{label_name}_at_far5pct"] = _recall_at_far(y_true, y_score, 0.05)
else:
horizon_metrics[f"pr_auc_{label_name}"] = None
horizon_metrics[f"brier_{label_name}"] = None
horizon_metrics[f"recall_{label_name}_at_far1pct"] = None
horizon_metrics[f"recall_{label_name}_at_far5pct"] = None
if label_name == "warn" and not skip_eval and reg is not None:
p_warning_test = _p_warning_series(clf, reg, X_test, thr, sigma_h)
per_horizon[h] = horizon_metrics
heads[f"max_{h}"] = reg
if not skip_eval and reg is not None and p_warning_test is not None:
observed_test_level = observed_grid.get((station, "water_level"))
if observed_test_level is not None:
observed_test_level = observed_test_level.loc[observed_test_level.index.isin(X_test.index)]
per_horizon[h]["events"] = _events_with_lead_time(observed_test_level, warn_thr, p_warning_test)
# --- full refit on the ENTIRE record for the deployed artifact ---
# This may include/exclude different heads than the eval-phase gate above (the
# full record has more labeled rows), so skip reasons are re-derived here --
# skipped_heads must reflect what actually ends up in the saved bundle.
final_heads: Dict[str, object] = {}
for h in horizons:
max_col, warn_col, danger_col = f"max_level_{h}", f"exceed_warn_{h}", f"exceed_danger_{h}"
head_key = f"max_{h}"
labeled = Y[max_col].notna()
if labeled.sum() >= MIN_ROWS_FOR_HEAD:
reg = _safe_fit(_make_regressor(hgb_overrides), X.loc[labeled], Y.loc[labeled, max_col], head_key, skipped_heads)
final_heads[head_key] = reg
if reg is not None:
skipped_heads.pop(head_key, None)
else:
skipped_heads[head_key] = f"only {int(labeled.sum())} labeled rows"
final_heads[head_key] = None
for label_name, col in (("warn", warn_col), ("danger", danger_col)):
head_key = f"{label_name}_{h}"
train_labeled = Y[col].notna()
n_pos = int(Y.loc[train_labeled, col].sum()) if train_labeled.any() else 0
if n_pos >= MIN_POSITIVES_FOR_CLASSIFIER:
clf = _safe_fit(
_make_classifier(hgb_overrides), X.loc[train_labeled], Y.loc[train_labeled, col], head_key, skipped_heads
)
final_heads[head_key] = clf
if clf is not None:
skipped_heads.pop(head_key, None)
else:
skipped_heads[head_key] = f"only {n_pos} positives in train span (< {MIN_POSITIVES_FOR_CLASSIFIER})"
final_heads[head_key] = None
bundle = {
"station_code": station,
"model_version": f"hgb-v1+{_git_short_sha()}",
"trained_at": datetime.datetime.now().isoformat(),
"sklearn_version": sklearn.__version__,
"feature_names": feature_names,
"horizons": list(horizons),
"thresholds": {"warning": warn_thr, "danger": danger_thr},
"heads": final_heads,
"sigma": sigma,
"skipped_heads": skipped_heads,
"train_span": meta["span"],
"n_train_rows": meta["n_rows"],
}
station_metrics = {"status": "trained", "per_horizon": per_horizon}
return bundle, station_metrics
def train_all(
df_long: pd.DataFrame,
stations: List[str],
horizons: Tuple[int, ...] = HORIZONS,
models_dir: Path = Path("models"),
skip_eval: bool = False,
hgb_overrides: Optional[dict] = None,
) -> dict:
"""Train and save every requested station's models. Returns the metrics.json payload."""
models_dir = Path(models_dir)
models_dir.mkdir(parents=True, exist_ok=True)
model_version = f"hgb-v1+{_git_short_sha()}"
station_results: Dict[str, dict] = {}
for station in stations:
if station in features.NOT_TRAINABLE:
reason = features.NOT_TRAINABLE[station]
logger.info(f"{station}: heuristic ({reason})")
station_results[station] = {"status": "heuristic", "reason": reason}
continue
try:
bundle, station_metrics = train_station(
df_long, station, horizons, skip_eval=skip_eval, hgb_overrides=hgb_overrides
)
if bundle is None:
logger.warning(f"{station}: failed ({station_metrics.get('reason')})")
station_results[station] = station_metrics
continue
joblib.dump(bundle, models_dir / f"flood_{station}.joblib")
logger.info(
f"{station}: trained, {bundle['n_train_rows']} rows, "
f"{len(bundle['skipped_heads'])} heads skipped"
)
station_results[station] = station_metrics
except Exception as error:
logger.error(f"{station}: failed with exception: {error}")
station_results[station] = {"status": "failed", "reason": str(error)}
metrics_payload = {
"generated_at": datetime.datetime.now().isoformat(),
"model_version": model_version,
"split": {
"train_end": SPLIT_B_TRAIN_END,
"test_start": SPLIT_B_TEST_START,
"test_end": SPLIT_B_TEST_END,
},
"stations": station_results,
}
with open(models_dir / "metrics.json", "w", encoding="utf-8") as handle:
json.dump(metrics_payload, handle, indent=2, default=str)
return metrics_payload
def main(argv: Optional[List[str]] = None) -> None:
logging.basicConfig(level=logging.INFO, format="%(asctime)s %(levelname)s %(name)s: %(message)s")
parser = argparse.ArgumentParser(description="Train Ping River flood forecast models")
parser.add_argument("--stations", default="all", help="'all' or a comma-separated list of station codes")
parser.add_argument("--models-dir", default="models")
parser.add_argument("--db-url", default=None)
parser.add_argument("--api-url", default=DEFAULT_API_URL)
parser.add_argument("--skip-eval", action="store_true", help="Refit-only fast path; skip Split B evaluation")
parser.add_argument("--start", default=None, help="ISO date; earliest measurement to load")
parser.add_argument("--end", default=None, help="ISO date; latest measurement to load")
args = parser.parse_args(argv)
if args.stations == "all":
stations = list(features.UPSTREAM_LEADS.keys())
else:
stations = [s.strip() for s in args.stations.split(",") if s.strip()]
start = datetime.datetime.fromisoformat(args.start) if args.start else None
end = datetime.datetime.fromisoformat(args.end) if args.end else None
logger.info(f"Loading measurements for {len(stations)} stations...")
df_long = load_measurements(
db_url=resolve_db_url(args.db_url), stations=None, start=start, end=end, api_url=args.api_url
)
logger.info(f"Loaded {len(df_long)} rows spanning {df_long['timestamp'].min()} .. {df_long['timestamp'].max()}")
metrics_payload = train_all(
df_long, stations, models_dir=Path(args.models_dir), skip_eval=args.skip_eval
)
trained = sum(1 for s in metrics_payload["stations"].values() if s["status"] == "trained")
logger.info(f"Done: {trained}/{len(stations)} stations trained. metrics.json written to {args.models_dir}")
if __name__ == "__main__":
main()