eval: quantile heads and fc48 on top of hgb-v3 (rejected/deferred); HII gauge-rain aggregate

Rolling-origin harness gains rise_rain_quantile, rise_rain_quantile_uw,
rise_rain_qsigma (L2 point + quantile sigma) and rise_rain_fc48, all
opt-in, plus --from-cache for reproducible offline reruns. Results in
models/eval_2026-09-12*.json, write-up in docs/FLOOD_FORECASTING.md:

- quantile point prediction: better MAE, worse first-alert lead at 5 of
  11 events (P.103 2022-08-14 +6h -> +1h) -> rejected
- quantile sigma only: Brier within noise (0.0031 -> 0.0029) -> not worth 3x heads
- rain_fc48: neutral everywhere except 2024-10-03 P.1 (+21h -> +72h), n=1
  -> deferred to after the 2026 season

src/ml/hii_rain.py: catchment-mean hourly rain from the ~130 HII gauges in
the upper-Ping box and a 24h-sum comparison against Open-Meteo. Not a
training feature (table exists only since 2026-08-11, no archive); exposed
at GET /api/hii/rainfall/catchment so the two sources' agreement is on
record by the time a fold can test it.

data._read_cache now skips non-station files in models/cache/ (the shared
dir also holds rain_openmeteo / dam_* caches, which crashed the reader).
scripts/summarize_eval.py prints per-variant lead/peak-error tables.

src/ml/data.py

@@ -282,6 +282,10 @@ def _read_cache(cache_dir: Path, stations: Optional[List[str]]) -> pd.DataFrame:
     frames = []
     for path in sorted(cache_dir.glob("*.csv.gz")):
         code = path.name[: -len(".csv.gz")]
+        # The dir is shared with rain.py / dam.py caches (rain_openmeteo,
+        # dam_<id>): only station files (P.<n>) are measurements.
+        if not code.startswith("P."):
+            continue
         if stations and code not in stations:
             continue
         with gzip.open(path, "rt", encoding="utf-8") as handle:

src/ml/evaluate.py

@@ -52,18 +52,31 @@ def _flood_weights(y_abs: pd.Series) -> np.ndarray:
     return 1.0 + 4.0 * np.clip((y_abs.to_numpy() - 2.5) / 1.2, 0.0, 1.0)
 
 
+# Experimental forward-48h rain sum, built in evaluate_station (not in
+# features.build_features) so the served feature set is untouched until the
+# harness says it helps. Serving could supply it: fetch_forecast() already
+# pulls forecast_days=2.
+EXTRA_RAIN_FEATURES = ("rain_fc48",)
+
+
 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,
-                 use_dam: bool = False):
+                 use_dam: bool = False, use_fc48: bool = False,
+                 qsigma: 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
+        self.use_fc48 = use_fc48
+        # Hybrid: L2 head for the point prediction (keeps the lead-time
+        # behaviour of the deployed model exactly, since p>=0.5 alerts are
+        # sigma-independent) and quantile heads ONLY for a per-row sigma.
+        self.qsigma = qsigma
 
     def fit_predict(
         self, X_tr, y_abs_tr, X_te
@@ -84,6 +97,12 @@ class Variant:
             raise ValueError(
                 f"{self.name} requires the dam series (rid_reservoir_daily backfilled)"
             )
+        if not self.use_fc48:
+            drop = [c for c in EXTRA_RAIN_FEATURES if c in X_tr.columns]
+            X_tr = X_tr.drop(columns=drop)
+            X_te = X_te.drop(columns=drop)
+        elif "rain_fc48" not in X_tr.columns:
+            raise ValueError(f"{self.name} requires the rain series")
         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
@@ -99,7 +118,13 @@ class Variant:
         else:
             reg = _make_regressor().fit(X_tr, y_tr, sample_weight=weights)
             pred = reg.predict(X_te)
-            sigma = np.full(len(X_te), FIXED_SIGMA)
+            if self.qsigma:
+                q50 = _quantile_regressor(0.5).fit(X_tr, y_tr, sample_weight=weights)
+                q90 = _quantile_regressor(0.9).fit(X_tr, y_tr, sample_weight=weights)
+                spread = np.maximum(q90.predict(X_te) - q50.predict(X_te), 0.0)
+                sigma = np.maximum(spread / 1.2816, 0.05)
+            else:
+                sigma = np.full(len(X_te), FIXED_SIGMA)
 
         pred_abs = pred + level_te if self.target == "rise" else pred
         pred_abs = np.maximum(pred_abs, level_te)  # peak >= current, as served
@@ -116,12 +141,29 @@ VARIANTS: Dict[str, Variant] = {
     "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),
+    # 2026-09-12 experiments on top of the deployed rise_rain configuration:
+    # per-row sigma from quantile heads (the served sigma sits on the 0.15
+    # floor at every P.1 horizon, so stage probabilities are constant-
+    # calibrated), and a longer forecast-rain window for the 24 h horizon.
+    "rise_rain_quantile": Variant("rise_rain_quantile", target="rise",
+                                  weighted=True, quantile=True, use_rain=True),
+    "rise_rain_quantile_uw": Variant("rise_rain_quantile_uw", target="rise",
+                                     quantile=True, use_rain=True),
+    "rise_rain_fc48": Variant("rise_rain_fc48", target="rise", use_rain=True,
+                              use_fc48=True),
+    "rise_rain_qsigma": Variant("rise_rain_qsigma", target="rise", use_rain=True,
+                                qsigma=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]
+# the 2026-08-13 ablation concluded them a negative result. The 2026-09-12
+# experiments are opt-in too (see their results in docs/FLOOD_FORECASTING.md).
+DEFAULT_VARIANTS = [
+    k for k, v in VARIANTS.items()
+    if not v.use_dam and not v.use_fc48 and not v.qsigma
+    and not (v.quantile and v.use_rain)
+]
 
 
 def _find_events(observed: pd.Series, thr: float) -> List[dict]:
@@ -222,6 +264,12 @@ def evaluate_station(
     warn_thr, _ = features.get_thresholds(station)
     grid = features.make_hourly_grid(df_long)
     X_all = features.build_features(grid, station, rain=rain, dam=dam)
+    if rain is not None:
+        # forward sum over (t, t+48]; same construction as rain_fc24
+        r = rain.reindex(X_all.index)
+        X_all["rain_fc48"] = (
+            r.shift(-1).iloc[::-1].rolling(48, min_periods=1).sum().iloc[::-1]
+        )
     observed = grid.observed[(station, "water_level")]
 
     keep = X_all["obs_age_h"].notna()
@@ -380,12 +428,19 @@ def main(argv=None) -> int:
                         help="skip loading the Open-Meteo rain series")
     parser.add_argument("--no-dam", action="store_true",
                         help="skip loading the Mae Ngat reservoir series")
+    parser.add_argument("--from-cache", action="store_true",
+                        help="offline: read models/cache/ only (no DB, no API, "
+                        "no Open-Meteo refresh) -- reproducible reruns")
     args = parser.parse_args(argv)
 
     logging.basicConfig(
         level=logging.INFO, format="%(asctime)s %(levelname)s %(message)s"
     )
-    df = data.load_measurements(db_url=args.db_url)
+    if args.from_cache:
+        df = data._read_cache(data.CACHE_DIR, None)
+        logger.info(f"measurements from cache: {len(df)} rows")
+    else:
+        df = data.load_measurements(db_url=args.db_url)
     if df.empty:
         logger.error("no measurement data")
         return 1
@@ -394,7 +449,9 @@ def main(argv=None) -> int:
     if not args.no_rain:
         from . import rain as rain_mod
 
-        rain_series = rain_mod.catchment_mean(rain_mod.load_history())
+        rain_series = rain_mod.catchment_mean(
+            rain_mod.load_history(refresh=not args.from_cache)
+        )
         if rain_series is None:
             logger.warning("rain history unavailable; rain features will be NaN")
         else:
@@ -404,7 +461,7 @@ def main(argv=None) -> int:
             )
 
     dam_frame = None
-    if not args.no_dam:
+    if not args.no_dam and not args.from_cache:
         from . import dam as dam_mod
 
         dam_frame = dam_mod.load_history(db_url=args.db_url)

src/ml/hii_rain.py

@@ -0,0 +1,119 @@
+"""Catchment-mean hourly rain from the HII/ThaiWater gauge network.
+
+Independent of Open-Meteo (src/ml/rain.py): those are model-analysis values,
+these are what the gauges measured. The `hii_rainfall` table has been filled
+by the hourly collector since 2026-08-11 and there is NO archive behind it
+(the api-v3 rain_24h_graph endpoint ignores its date range, see
+docs/DATA_SOURCES.md 2.1), so this series cannot yet be a training feature:
+every training row before 2026-08 would be NaN and HistGradientBoosting
+would learn nothing from the column. It becomes a candidate once a full
+monsoon season of gauge rows exists in the rolling-origin harness's test
+span -- the 2027 fold (train through 2027-04-30, test Jun-Nov 2027) is the
+first that could show anything.
+
+Until then it serves two purposes:
+  * a live cross-check of the Open-Meteo catchment mean (/api/hii/rainfall
+    already exposes the raw gauges; this gives the comparable aggregate);
+  * accumulating the comparison so the eventual feature evaluation has a
+    documented bias/variance relationship between the two sources.
+"""
+
+import logging
+from typing import Optional, Sequence, Tuple
+
+import pandas as pd
+
+from .data import resolve_db_url
+
+logger = logging.getLogger(__name__)
+
+# Same footprint as rain.CATCHMENT_POINTS: the upper Ping above P.1. Gauges
+# inside this box are averaged; there are ~130 with recent data (DWR, FOP,
+# HII, RID, TMD), far denser than the five Open-Meteo points.
+CATCHMENT_BOX: Tuple[float, float, float, float] = (18.75, 19.60, 98.60, 99.30)
+# A gauge that reports the same rain_24h for many hours is stuck; drop hours
+# where fewer than this many gauges reported at all.
+MIN_GAUGES_PER_HOUR = 5
+
+
+def load_gauge_mean(
+    db_url: Optional[str] = None,
+    start: Optional[pd.Timestamp] = None,
+    end: Optional[pd.Timestamp] = None,
+    box: Sequence[float] = CATCHMENT_BOX,
+    engine=None,
+) -> Optional[pd.Series]:
+    """Hourly catchment-mean rain_1h (mm) across HII gauges in `box`.
+
+    Pass `engine` (the API's HII store engine) to reuse a pool; otherwise a
+    connection is resolved from db_url / config. Returns None if the DB is
+    unavailable or the table is empty. Hours with fewer than
+    MIN_GAUGES_PER_HOUR reporting gauges are NaN.
+    """
+    if engine is None:
+        resolved = resolve_db_url(db_url)
+        if not resolved:
+            return None
+    lat_lo, lat_hi, lon_lo, lon_hi = box
+    try:
+        from sqlalchemy import create_engine, text
+
+        query = (
+            "SELECT m.timestamp, COUNT(m.rain_1h) AS n, AVG(m.rain_1h) AS rain_1h "
+            "FROM hii_rainfall m JOIN hii_rain_stations s ON s.id = m.station_id "
+            "WHERE s.latitude BETWEEN :lat_lo AND :lat_hi "
+            "AND s.longitude BETWEEN :lon_lo AND :lon_hi "
+            "AND m.rain_1h IS NOT NULL"
+        )
+        params = {"lat_lo": lat_lo, "lat_hi": lat_hi, "lon_lo": lon_lo, "lon_hi": lon_hi}
+        if start is not None:
+            query += " AND m.timestamp >= :start"
+            params["start"] = pd.Timestamp(start).to_pydatetime()
+        if end is not None:
+            query += " AND m.timestamp <= :end"
+            params["end"] = pd.Timestamp(end).to_pydatetime()
+        query += " GROUP BY m.timestamp ORDER BY m.timestamp"
+        if engine is None:
+            engine = create_engine(resolved, pool_pre_ping=True)
+        with engine.connect() as conn:
+            frame = pd.read_sql(text(query), conn, params=params)
+    except Exception as error:
+        logger.warning(f"HII gauge rain load failed: {error}")
+        return None
+    if frame.empty:
+        return None
+    frame["timestamp"] = pd.to_datetime(frame["timestamp"]).dt.floor("h")
+    frame = frame.groupby("timestamp").agg(n=("n", "sum"), rain_1h=("rain_1h", "mean"))
+    series = pd.to_numeric(frame["rain_1h"], errors="coerce")
+    series[frame["n"] < MIN_GAUGES_PER_HOUR] = float("nan")
+    series.name = "hii_gauge_mean"
+    return series
+
+
+def compare_with_openmeteo(
+    gauge: pd.Series, openmeteo: pd.Series, window_h: int = 24
+) -> dict:
+    """Bias/correlation of Open-Meteo against the gauges over the overlap.
+
+    Both are summed over trailing `window_h` so single-hour timing offsets
+    (gauges report at :00, the model's hour is an interval) do not dominate.
+    """
+    joined = pd.concat(
+        {"gauge": gauge, "openmeteo": openmeteo}, axis=1
+    ).dropna()
+    if joined.empty:
+        return {"overlap_hours": 0}
+    g = joined["gauge"].rolling(window_h, min_periods=window_h).sum()
+    o = joined["openmeteo"].rolling(window_h, min_periods=window_h).sum()
+    both = pd.concat({"g": g, "o": o}, axis=1).dropna()
+    if both.empty:
+        return {"overlap_hours": int(len(joined))}
+    return {
+        "overlap_hours": int(len(joined)),
+        "window_h": window_h,
+        "gauge_mean_mm": float(both["g"].mean()),
+        "openmeteo_mean_mm": float(both["o"].mean()),
+        "bias_mm": float((both["o"] - both["g"]).mean()),
+        "mae_mm": float((both["o"] - both["g"]).abs().mean()),
+        "corr": float(both["g"].corr(both["o"])),
+    }