feat: Mae Ngat dam features — built, evaluated, defaulted OFF

src/ml/dam.py loads rid_reservoir_daily into a leakage-safe hourly frame
(daily row visible from 07:00 its own date, ffill capped at 48 h) and is
plumbed through features/train/predict/evaluate exactly like rain, gated
to the six mainstem stations below the Mae Ngat confluence.

The experiment concludes as a documented NEGATIVE result: on the 2024
record-flood backtest every dam-feature subset costs 1-3 h of first-alert
lead (13h -> 10-12h) for <=3 cm of peak-error gain, because the daily RID
report lags up to 31 h and describes yesterday's benign absorbing
reservoir during fast onset. Features therefore default OFF (--dam
opt-in on the training and backtest CLIs; rise_rain_dam/rise_dam harness
variants, excluded from the default variant set). The ablation also
isolated the HII gap-fill as lead-neutral: the acceptance gate holds at
13 h with fill enabled, and docs/img charts are regenerated with the
shipping configuration. Full table in docs/FLOOD_FORECASTING.md §5.

Review-swarm fixes: evaluate.py skips variants whose feature family is
absent instead of crashing the run; --dam forwards --db-url and warns
loudly when no dam history loads; an empty DB result can no longer wipe
a good dam cache; run-level metrics version claims v4 only when a dam
station is actually in the set.

src/ml/dam.py

@@ -0,0 +1,111 @@
+"""Mae Ngat reservoir series for the flood models.
+
+rid_reservoir_daily (collected hourly by src/rid_reservoir.py, backfilled to
+2018) holds daily storage/inflow/outflow for every RID large dam. Mae Ngat
+Somboon Chon (DAM_ID 200103) is the only large dam upstream of Chiang Mai:
+in Oct 2024 its inflow hit 19-22 MCM/day and storage 114% of usable capacity
+days around the P.1 crossing — upstream state no river gauge carries.
+
+Leakage rule: RID publishes the daily report for date D on the morning of D,
+so the row becomes visible to features at D 07:00 local time, never earlier.
+Forward-fill is capped at FFILL_LIMIT_H so a stalled collector degrades to
+NaN (HGB-native) instead of silently serving stale reservoir state.
+
+Known residual optimism: the collector upserts keep-last (and re-fetches
+yesterday), so the stored row for date D is RID's FINAL revision, which
+training then back-dates to D 07:00 — values live serving may not have had
+that morning. This bias works IN FAVOR of dam features, so the 2026-08-13
+negative result (they cost 1-3 h of alert lead) holds a fortiori; but any
+future POSITIVE result must first validate intraday row stability or shift
+the flow columns to D+1 07:00.
+"""
+
+import datetime
+import logging
+from pathlib import Path
+from typing import Optional
+
+import pandas as pd
+
+from ..rid_reservoir import MAE_NGAT_DAM_ID
+from .data import CACHE_DIR, resolve_db_url
+
+logger = logging.getLogger(__name__)
+
+REPORT_HOUR = 7  # daily value valid from 07:00 local on its own date
+FFILL_LIMIT_H = 48  # two missed daily reports -> NaN, not stale state
+DAM_COLUMNS = ("storage_pct", "inflow_mcm", "outflow_mcm")
+CACHE_FILE = f"dam_{MAE_NGAT_DAM_ID}.csv.gz"
+
+
+def load_daily(
+    db_url: Optional[str] = None,
+    dam_id: str = MAE_NGAT_DAM_ID,
+    start: Optional[datetime.date] = None,
+    cache_dir: Path = CACHE_DIR,
+) -> Optional[pd.DataFrame]:
+    """Daily dam rows indexed by date. DB first, on-disk cache as fallback."""
+    cache_path = Path(cache_dir) / CACHE_FILE
+    resolved = resolve_db_url(db_url)
+    if resolved:
+        try:
+            from sqlalchemy import create_engine, text
+
+            query = (
+                "SELECT date, storage_pct, inflow_mcm, outflow_mcm "
+                "FROM rid_reservoir_daily WHERE dam_id = :dam_id"
+            )
+            params = {"dam_id": dam_id}
+            if start is not None:
+                query += " AND date >= :start"
+                params["start"] = start
+            engine = create_engine(resolved, pool_pre_ping=True)
+            with engine.connect() as conn:
+                daily = pd.read_sql(
+                    text(query + " ORDER BY date"), conn, params=params
+                )
+            daily["date"] = pd.to_datetime(daily["date"])
+            daily = daily.set_index("date")
+            for col in DAM_COLUMNS:
+                daily[col] = pd.to_numeric(daily[col], errors="coerce")
+            # Only full, NON-EMPTY loads refresh the cache: a truncated or
+            # freshly-recreated table must not wipe a good fallback archive.
+            if start is None and not daily.empty:
+                cache_path.parent.mkdir(parents=True, exist_ok=True)
+                daily.to_csv(cache_path, compression="gzip")
+            return daily
+        except Exception as error:
+            logger.warning(f"dam series DB load failed: {error}")
+    if cache_path.exists():
+        logger.warning("falling back to on-disk cache for the dam series")
+        return pd.read_csv(cache_path, index_col=0, parse_dates=True)
+    return None
+
+
+def hourly_frame(daily: Optional[pd.DataFrame]) -> Optional[pd.DataFrame]:
+    """Step the daily rows onto an hourly grid, each valid from D 07:00."""
+    if daily is None or daily.empty:
+        return None
+    frame = daily.copy()
+    frame.index = pd.to_datetime(frame.index) + pd.Timedelta(hours=REPORT_HOUR)
+    frame = frame[~frame.index.duplicated(keep="last")].sort_index()
+    hourly_index = pd.date_range(
+        frame.index.min(),
+        frame.index.max() + pd.Timedelta(hours=FFILL_LIMIT_H),
+        freq="h",
+    )
+    return frame.reindex(hourly_index).ffill(limit=FFILL_LIMIT_H)
+
+
+def load_history(db_url: Optional[str] = None) -> Optional[pd.DataFrame]:
+    """Full hourly Mae Ngat history for training; None when unavailable."""
+    return hourly_frame(load_daily(db_url))
+
+
+def serving_frame(
+    db_url: Optional[str] = None, days: int = 21
+) -> Optional[pd.DataFrame]:
+    """Recent hourly dam state for inference (covers the 336 h feature window
+    plus the 72 h storage-delta lag)."""
+    start = datetime.date.today() - datetime.timedelta(days=days)
+    return hourly_frame(load_daily(db_url, start=start))

src/ml/evaluate.py

@@ -56,12 +56,14 @@ class Variant:
     """A trainable candidate producing (pred_abs, sigma_per_row) on test rows."""
 
     def __init__(self, name: str, target: str, weighted: bool = False,
-                 quantile: bool = False, use_rain: bool = False):
+                 quantile: bool = False, use_rain: bool = False,
+                 use_dam: bool = False):
         self.name = name
         self.target = target  # 'abs' or 'rise'
         self.weighted = weighted
         self.quantile = quantile
         self.use_rain = use_rain
+        self.use_dam = use_dam
 
     def fit_predict(
         self, X_tr, y_abs_tr, X_te
@@ -74,6 +76,14 @@ class Variant:
             raise ValueError(
                 f"{self.name} requires the rain series (run without --no-rain)"
             )
+        if not self.use_dam:
+            drop = [c for c in features.DAM_FEATURES if c in X_tr.columns]
+            X_tr = X_tr.drop(columns=drop)
+            X_te = X_te.drop(columns=drop)
+        elif "dam_storage_pct" not in X_tr.columns:
+            raise ValueError(
+                f"{self.name} requires the dam series (rid_reservoir_daily backfilled)"
+            )
         level_tr = X_tr["level"]
         level_te = X_te["level"].to_numpy()
         y_tr = (y_abs_tr - level_tr) if self.target == "rise" else y_abs_tr
@@ -103,8 +113,16 @@ VARIANTS: Dict[str, Variant] = {
     "rise_quantile": Variant("rise_quantile", target="rise", weighted=True,
                              quantile=True),
     "rise_rain": Variant("rise_rain", target="rise", use_rain=True),
+    "rise_rain_dam": Variant("rise_rain_dam", target="rise", use_rain=True,
+                             use_dam=True),
+    "rise_dam": Variant("rise_dam", target="rise", use_dam=True),
 }
 
+# Dam variants are opt-in by name: they require dam columns that only exist
+# for features.DAM_STATIONS and only when the reservoir series loaded, and
+# the 2026-08-13 ablation concluded them a negative result.
+DEFAULT_VARIANTS = [k for k, v in VARIANTS.items() if not v.use_dam]
+
 
 def _find_events(observed: pd.Series, thr: float) -> List[dict]:
     """Contiguous >=thr episodes (gaps under EVENT_GAP_H merged)."""
@@ -198,11 +216,12 @@ def evaluate_station(
     variants: Optional[List[str]] = None,
     seasons: Tuple[int, ...] = SEASONS,
     rain: Optional[pd.Series] = None,
+    dam: Optional[pd.DataFrame] = None,
 ) -> Dict:
     """Run every fold x variant for one station; returns the results tree."""
     warn_thr, _ = features.get_thresholds(station)
     grid = features.make_hourly_grid(df_long)
-    X_all = features.build_features(grid, station, rain=rain)
+    X_all = features.build_features(grid, station, rain=rain, dam=dam)
     observed = grid.observed[(station, "water_level")]
 
     keep = X_all["obs_age_h"].notna()
@@ -211,7 +230,7 @@ def evaluate_station(
         keep &= X_all.index >= pd.Timestamp(train_start)
     X_all = X_all.loc[keep]
 
-    chosen = {k: VARIANTS[k] for k in (variants or VARIANTS)}
+    chosen = {k: VARIANTS[k] for k in (variants or DEFAULT_VARIANTS)}
     results: Dict = {"station": station, "warn_thr": warn_thr, "folds": []}
 
     for year in seasons:
@@ -253,7 +272,14 @@ def evaluate_station(
         }
 
         for name, variant in chosen.items():
-            pred_abs, sigma = variant.fit_predict(X_tr, y_tr, X_te)
+            try:
+                pred_abs, sigma = variant.fit_predict(X_tr, y_tr, X_te)
+            except ValueError as error:
+                # A variant whose required feature family is absent (e.g. a
+                # dam variant on a non-DAM_STATIONS target) skips this fold
+                # instead of killing the whole run and its finished results.
+                logger.warning(f"{station} {year} {name}: skipped ({error})")
+                continue
             pred_series = pd.Series(pred_abs, index=X_te.index)
             p_warn = pd.Series(
                 1.0 - _phi((warn_thr - pred_abs) / sigma), index=X_te.index
@@ -352,6 +378,8 @@ def main(argv=None) -> int:
     parser.add_argument("--out", default="models/eval_variants.json")
     parser.add_argument("--no-rain", action="store_true",
                         help="skip loading the Open-Meteo rain series")
+    parser.add_argument("--no-dam", action="store_true",
+                        help="skip loading the Mae Ngat reservoir series")
     args = parser.parse_args(argv)
 
     logging.basicConfig(
@@ -375,12 +403,27 @@ def main(argv=None) -> int:
                 f"{rain_series.index.max()}"
             )
 
+    dam_frame = None
+    if not args.no_dam:
+        from . import dam as dam_mod
+
+        dam_frame = dam_mod.load_history(db_url=args.db_url)
+        if dam_frame is None:
+            logger.warning("dam history unavailable; dam features will be absent")
+        else:
+            logger.info(
+                f"dam series loaded: {dam_frame.index.min()} .. "
+                f"{dam_frame.index.max()}"
+            )
+
     variant_names = args.variants.split(",") if args.variants else None
     all_results = []
     for station in args.stations.split(","):
         station = station.strip()
         logger.info(f"Evaluating {station}...")
-        results = evaluate_station(df, station, variant_names, rain=rain_series)
+        results = evaluate_station(
+            df, station, variant_names, rain=rain_series, dam=dam_frame
+        )
         all_results.append(results)
         print(summarize(results))
 

src/ml/features.py

@@ -202,9 +202,18 @@ def _hours_since_observed(mask_col: pd.Series) -> pd.Series:
 
 RAIN_FEATURES = ("rain_6h", "rain_24h", "rain_72h", "rain_fc24")
 
+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
+    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.
 
@@ -217,6 +226,11 @@ def build_features(
     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] = {}
@@ -280,6 +294,18 @@ def build_features(
             r.shift(-1).iloc[::-1].rolling(24, min_periods=1).sum().iloc[::-1]
         )
 
+    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)
 
 
@@ -358,10 +384,11 @@ def build_matrix(
     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)
+    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()

src/ml/predict.py

@@ -128,6 +128,7 @@ def _model_forecast(
     as_of: pd.Timestamp,
     current_level: float,
     rain: Optional[pd.Series] = None,
+    dam: Optional[pd.DataFrame] = None,
 ) -> List[dict]:
     warn_thr = bundle["thresholds"]["warning"]
     danger_thr = bundle["thresholds"]["danger"]
@@ -146,7 +147,9 @@ def _model_forecast(
         )
         warn_thr, danger_thr = cfg_warn, cfg_danger
 
-    feature_row = features.build_features(grid, station_code, rain=rain).loc[[as_of]]
+    feature_row = features.build_features(grid, station_code, rain=rain, dam=dam).loc[
+        [as_of]
+    ]
     expected_columns = bundle["feature_names"]
     missing = [c for c in expected_columns if c not in feature_row.columns]
     if missing:
@@ -231,6 +234,7 @@ def _forecast_station(
     now: pd.Timestamp,
     horizons: Tuple[int, ...],
     rain: Optional[pd.Series] = None,
+    dam: Optional[pd.DataFrame] = None,
 ) -> List[dict]:
     level_col = (station_code, "water_level")
     if level_col not in grid.observed.columns:
@@ -267,7 +271,7 @@ def _forecast_station(
 
     bundle = _load_bundle(bundle_path)
     model_results = _model_forecast(
-        station_code, grid, bundle, as_of, current_level, rain=rain
+        station_code, grid, bundle, as_of, current_level, rain=rain, dam=dam
     )
     if model_results is None:
         return _heuristic_forecast(
@@ -306,6 +310,7 @@ def get_forecasts(
     models_dir: Union[str, Path] = DEFAULT_MODELS_DIR,
     now: Optional[Union[datetime.datetime, str]] = None,
     rain: Optional[pd.Series] = None,
+    dam: Optional[pd.DataFrame] = None,
 ) -> List[dict]:
     """Produce flood forecasts for every station present in `readings_by_station`.
 
@@ -329,7 +334,13 @@ def get_forecasts(
         try:
             results.extend(
                 _forecast_station(
-                    station_code, grid, models_dir, now, DEFAULT_HORIZONS, rain=rain
+                    station_code,
+                    grid,
+                    models_dir,
+                    now,
+                    DEFAULT_HORIZONS,
+                    rain=rain,
+                    dam=dam,
                 )
             )
         except Exception as error:
@@ -376,4 +387,19 @@ def get_latest_forecasts(
         logger.warning("live rain unavailable; rain features will be NaN")
         rain = pd.Series(dtype=float)
 
-    return get_forecasts(readings_by_station, models_dir=models_dir, rain=rain)
+    # Recent Mae Ngat reservoir state; same empty-not-None contract so
+    # dam-trained bundles keep their columns (NaN) when the DB read fails.
+    from . import dam as dam_mod
+
+    try:
+        dam = dam_mod.serving_frame(db_url=db_url)
+    except Exception as error:
+        logger.warning(f"dam serving frame failed: {error}")
+        dam = None
+    if dam is None:
+        logger.warning("dam state unavailable; dam features will be NaN")
+        dam = pd.DataFrame()
+
+    return get_forecasts(
+        readings_by_station, models_dir=models_dir, rain=rain, dam=dam
+    )

src/ml/train.py

@@ -204,9 +204,10 @@ def train_station(
     split_test_start: str = SPLIT_B_TEST_START,
     split_test_end: str = SPLIT_B_TEST_END,
     rain: Optional[pd.Series] = None,
+    dam: Optional[pd.DataFrame] = None,
 ) -> 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, rain=rain)
+    X, Y, meta = features.build_matrix(df_long, station, horizons, rain=rain, dam=dam)
     if meta["n_rows"] < MIN_ROWS_TO_TRAIN:
         return None, {
             "status": "failed",
@@ -412,8 +413,13 @@ def train_station(
                 ] = f"only {n_pos} positives in train span (< {MIN_POSITIVES_FOR_CLASSIFIER})"
                 final_heads[head_key] = None
 
-    # v3 = rise target + Open-Meteo rain features; v2 = rise target only
-    version_prefix = "hgb-v3" if "rain_24h" in feature_names else "hgb-v2"
+    # v4 = + Mae Ngat dam features; v3 = rise + rain; v2 = rise target only
+    if "dam_storage_pct" in feature_names:
+        version_prefix = "hgb-v4"
+    elif "rain_24h" in feature_names:
+        version_prefix = "hgb-v3"
+    else:
+        version_prefix = "hgb-v2"
     bundle = {
         "station_code": station,
         "model_version": f"{version_prefix}+{_git_short_sha()}",
@@ -443,6 +449,8 @@ def train_all(
     skip_eval: bool = False,
     hgb_overrides: Optional[dict] = None,
     use_rain: bool = True,
+    use_dam: bool = False,
+    db_url: Optional[str] = None,
 ) -> dict:
     """Train and save every requested station's models. Returns the metrics.json payload."""
     models_dir = Path(models_dir)
@@ -462,7 +470,42 @@ def train_all(
         logger.info(
             f"rain series: {rain_series.index.min()} .. {rain_series.index.max()}"
         )
-    version_prefix = "hgb-v3" if rain_series is not None else "hgb-v2"
+
+    # Mae Ngat reservoir state (rid_reservoir_daily, 2018+). OFF by default:
+    # the 2026-08-13 backtest ablation showed every dam-feature subset COSTS
+    # 1-3 h of first-alert lead on the 2024 record flood (the daily report
+    # lags up to 31 h, so during fast onset the columns describe yesterday's
+    # benign reservoir and damp the alarm). Kept as an opt-in for post-monsoon
+    # re-evaluation once the 2026 season adds dam-era flood events.
+    dam_frame = None
+    if use_dam:
+        try:
+            from . import dam as dam_mod
+
+            dam_frame = dam_mod.load_history(db_url=db_url)
+        except Exception as error:
+            logger.warning(f"dam history unavailable, training without it: {error}")
+        if dam_frame is None:
+            # load_history returns None (no raise) when both DB and cache
+            # miss — an explicitly requested experiment must say so loudly.
+            logger.warning(
+                "--dam requested but no dam history available; "
+                "training v3-style bundles WITHOUT dam features"
+            )
+    if dam_frame is not None:
+        logger.info(
+            f"dam series: {dam_frame.index.min()} .. {dam_frame.index.max()}"
+        )
+
+    # Run-level version: v4 only if some requested station actually receives
+    # dam columns (they are gated to DAM_STATIONS; per-bundle versions are
+    # derived from each station's own feature_names and remain authoritative).
+    if dam_frame is not None and any(s in features.DAM_STATIONS for s in stations):
+        version_prefix = "hgb-v4"
+    elif rain_series is not None:
+        version_prefix = "hgb-v3"
+    else:
+        version_prefix = "hgb-v2"
     model_version = f"{version_prefix}+{_git_short_sha()}"
 
     station_results: Dict[str, dict] = {}
@@ -480,6 +523,7 @@ def train_all(
                 skip_eval=skip_eval,
                 hgb_overrides=hgb_overrides,
                 rain=rain_series,
+                dam=dam_frame,
             )
             if bundle is None:
                 logger.warning(f"{station}: failed ({station_metrics.get('reason')})")
@@ -542,6 +586,12 @@ def main(argv: Optional[List[str]] = None) -> None:
         action="store_true",
         help="train without the Open-Meteo rain features (v2-style bundles)",
     )
+    parser.add_argument(
+        "--dam",
+        action="store_true",
+        help="EXPERIMENTAL: include Mae Ngat reservoir features (v4 bundles); "
+        "the 2026-08 ablation showed they cost 1-3 h of alert lead",
+    )
     args = parser.parse_args(argv)
 
     if args.stations == "all":
@@ -570,6 +620,8 @@ def main(argv: Optional[List[str]] = None) -> None:
         models_dir=Path(args.models_dir),
         skip_eval=args.skip_eval,
         use_rain=not args.no_rain,
+        use_dam=args.dam,
+        db_url=resolve_db_url(args.db_url),
     )
     trained = sum(
         1 for s in metrics_payload["stations"].values() if s["status"] == "trained"