feat: Open-Meteo catchment rainfall series + rain features + harness variant

src/ml/rain.py fetches hourly precipitation for five upper-Ping
catchment points (Chiang Dao, Mae Taeng, Mae Ngat, Mae Rim, city) from
the Open-Meteo forecast-model archive (2021-03 onward, no API key,
Bangkok-local timestamps, year-chunked local cache) plus the live
forecast endpoint for serving (trailing days + next 48h).

features.build_features/build_matrix accept the catchment-mean series
and add rain_6h/24h/72h trailing sums and rain_fc24 — the forward 24h
sum, a genuine forecast feature (archived forecasts at training time, a
real weather forecast at serving; never contains river data). Columns
exist only when a series is provided: HGB rejects all-NaN columns at
fit, so no-rain training omits them and serving passes an empty series
for alignment. evaluate.py gains the rise_rain variant (and --no-rain)
so the harness can judge whether rain beats the deployed rise baseline.

src/ml/evaluate.py

@@ -56,15 +56,24 @@ 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):
+                 quantile: bool = False, use_rain: bool = False):
         self.name = name
         self.target = target  # 'abs' or 'rise'
         self.weighted = weighted
         self.quantile = quantile
+        self.use_rain = use_rain
 
     def fit_predict(
         self, X_tr, y_abs_tr, X_te
     ) -> Tuple[np.ndarray, np.ndarray]:
+        if not self.use_rain:
+            drop = [c for c in features.RAIN_FEATURES if c in X_tr.columns]
+            X_tr = X_tr.drop(columns=drop)
+            X_te = X_te.drop(columns=drop)
+        elif "rain_24h" not in X_tr.columns:
+            raise ValueError(
+                f"{self.name} requires the rain series (run without --no-rain)"
+            )
         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
@@ -93,6 +102,7 @@ VARIANTS: Dict[str, Variant] = {
     "rise_weighted": Variant("rise_weighted", target="rise", weighted=True),
     "rise_quantile": Variant("rise_quantile", target="rise", weighted=True,
                              quantile=True),
+    "rise_rain": Variant("rise_rain", target="rise", use_rain=True),
 }
 
 
@@ -187,11 +197,12 @@ def evaluate_station(
     station: str,
     variants: Optional[List[str]] = None,
     seasons: Tuple[int, ...] = SEASONS,
+    rain: Optional[pd.Series] = 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)
+    X_all = features.build_features(grid, station, rain=rain)
     observed = grid.observed[(station, "water_level")]
 
     keep = X_all["obs_age_h"].notna()
@@ -339,6 +350,8 @@ def main(argv=None) -> int:
     parser.add_argument("--variants", default=None,
                         help="comma list; default all")
     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")
     args = parser.parse_args(argv)
 
     logging.basicConfig(
@@ -349,12 +362,25 @@ def main(argv=None) -> int:
         logger.error("no measurement data")
         return 1
 
+    rain_series = None
+    if not args.no_rain:
+        from . import rain as rain_mod
+
+        rain_series = rain_mod.catchment_mean(rain_mod.load_history())
+        if rain_series is None:
+            logger.warning("rain history unavailable; rain features will be NaN")
+        else:
+            logger.info(
+                f"rain series loaded: {rain_series.index.min()} .. "
+                f"{rain_series.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)
+        results = evaluate_station(df, station, variant_names, rain=rain_series)
         all_results.append(results)
         print(summarize(results))
 

src/ml/features.py

@@ -200,8 +200,24 @@ def _hours_since_observed(mask_col: pd.Series) -> pd.Series:
 # ---------------------------------------------------------------------------
 
 
-def build_features(grid: HourlyGrid, station: str) -> pd.DataFrame:
-    """Build the deterministic-order feature matrix for one target station."""
+RAIN_FEATURES = ("rain_6h", "rain_24h", "rain_72h", "rain_fc24")
+
+
+def build_features(
+    grid: HourlyGrid, station: str, rain: Optional[pd.Series] = 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.
+    """
     idx = grid.observed.index
     cols: Dict[str, pd.Series] = {}
 
@@ -254,6 +270,16 @@ def build_features(grid: HourlyGrid, station: str) -> pd.DataFrame:
     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
+        cols["rain_fc24"] = (
+            r.shift(-1).iloc[::-1].rolling(24, min_periods=1).sum().iloc[::-1]
+        )
+
     return pd.DataFrame(cols, index=idx)
 
 
@@ -331,10 +357,11 @@ def build_matrix(
     station: str,
     horizons: Tuple[int, ...] = (6, 12, 24),
     stats_end: Optional[str] = None,
+    rain: Optional[pd.Series] = 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)
+    X = build_features(grid, station, rain=rain)
     Y = build_labels(grid, station, horizons, stats_end=stats_end)
 
     keep = X["obs_age_h"].notna()

src/ml/rain.py

@@ -0,0 +1,162 @@
+"""Open-Meteo rainfall series for the upper Ping catchment.
+
+One consistent source for training AND serving: the Open-Meteo forecast-model
+archive (historical-forecast-api, 2021-03 onward) supplies hourly
+precipitation at five catchment points above P.1; the live forecast endpoint
+supplies the same series for recent days plus the next 48 h. Timestamps are
+Asia/Bangkok local, matching the measurement grid. Rows before 2021-03 simply
+have no rain data — HistGradientBoosting handles the NaNs natively.
+
+The forward-looking sum built from this series is a legitimate *forecast*
+feature, not label leakage: the series never contains river observations, and
+at serving time the future values come from an actual weather forecast.
+"""
+
+import datetime
+import logging
+from pathlib import Path
+from typing import Iterable, List, Optional, Tuple
+
+import pandas as pd
+import requests
+
+logger = logging.getLogger(__name__)
+
+# (name, lat, lon) — upper Ping catchment above P.1, headwaters to city
+CATCHMENT_POINTS: Tuple[Tuple[str, float, float], ...] = (
+    ("chiang_dao", 19.37, 98.97),
+    ("mae_taeng", 19.12, 98.94),
+    ("mae_ngat", 19.17, 99.05),
+    ("mae_rim", 18.92, 98.92),
+    ("chiang_mai", 18.79, 99.00),
+)
+
+HISTORY_URL = "https://historical-forecast-api.open-meteo.com/v1/forecast"
+FORECAST_URL = "https://api.open-meteo.com/v1/forecast"
+HISTORY_START = "2021-03-23"  # archive begins here
+CACHE_FILE = "rain_openmeteo.csv.gz"
+
+
+def _points_params() -> dict:
+    return {
+        "latitude": ",".join(str(lat) for _, lat, _ in CATCHMENT_POINTS),
+        "longitude": ",".join(str(lon) for _, _, lon in CATCHMENT_POINTS),
+        "hourly": "precipitation",
+        "timezone": "Asia/Bangkok",
+    }
+
+
+def _parse_multi(payload, columns: Iterable[str]) -> pd.DataFrame:
+    """Open-Meteo returns a list when multiple coordinates are requested."""
+    results = payload if isinstance(payload, list) else [payload]
+    frames = {}
+    for name, result in zip(columns, results):
+        hourly = result.get("hourly", {})
+        idx = pd.to_datetime(hourly.get("time", []))
+        frames[name] = pd.Series(hourly.get("precipitation", []), index=idx)
+    df = pd.DataFrame(frames)
+    df.index.name = "timestamp"
+    return df
+
+
+def fetch_history(
+    start: str, end: str, session: Optional[requests.Session] = None
+) -> pd.DataFrame:
+    """Hourly precipitation for all catchment points over [start, end]."""
+    session = session or requests.Session()
+    response = session.get(
+        HISTORY_URL,
+        params={**_points_params(), "start_date": start, "end_date": end},
+        timeout=120,
+    )
+    response.raise_for_status()
+    return _parse_multi(response.json(), [p[0] for p in CATCHMENT_POINTS])
+
+
+def fetch_forecast(
+    past_days: int = 5,
+    forecast_days: int = 2,
+    session: Optional[requests.Session] = None,
+) -> pd.DataFrame:
+    """Recent + next-48h precipitation from the live forecast endpoint."""
+    session = session or requests.Session()
+    response = session.get(
+        FORECAST_URL,
+        params={
+            **_points_params(),
+            "past_days": past_days,
+            "forecast_days": forecast_days,
+        },
+        timeout=60,
+    )
+    response.raise_for_status()
+    return _parse_multi(response.json(), [p[0] for p in CATCHMENT_POINTS])
+
+
+def load_history(
+    cache_dir: Path = Path("models/cache"),
+    end: Optional[datetime.date] = None,
+    refresh: bool = True,
+) -> Optional[pd.DataFrame]:
+    """Cached catchment rain history from 2021-03 to ~today.
+
+    Fetches year-sized chunks on first use (~6 requests), then only extends
+    the tail. Returns None when the API is unreachable and no cache exists.
+    """
+    cache_dir.mkdir(parents=True, exist_ok=True)
+    cache_path = cache_dir / CACHE_FILE
+    end = end or datetime.date.today()
+
+    cached: Optional[pd.DataFrame] = None
+    if cache_path.exists():
+        cached = pd.read_csv(cache_path, index_col=0, parse_dates=True)
+
+    fetch_from = pd.Timestamp(HISTORY_START)
+    if cached is not None and len(cached):
+        fetch_from = cached.index.max() - pd.Timedelta(days=2)  # re-fetch tail
+    if not refresh and cached is not None:
+        return cached
+
+    chunks: List[pd.DataFrame] = []
+    cursor = fetch_from.date()
+    try:
+        while cursor <= end:
+            chunk_end = min(
+                datetime.date(cursor.year, 12, 31), end
+            )
+            chunks.append(
+                fetch_history(cursor.isoformat(), chunk_end.isoformat())
+            )
+            cursor = datetime.date(cursor.year + 1, 1, 1)
+    except Exception as error:
+        logger.warning(f"Open-Meteo history fetch failed: {error}")
+        if not chunks and cached is None:
+            return None
+
+    if chunks:
+        fresh = pd.concat(chunks)
+        combined = (
+            pd.concat([cached[cached.index < fresh.index.min()], fresh])
+            if cached is not None
+            else fresh
+        )
+        combined = combined[~combined.index.duplicated(keep="last")].sort_index()
+        combined.to_csv(cache_path, compression="gzip")
+        return combined
+    return cached
+
+
+def catchment_mean(df: Optional[pd.DataFrame]) -> Optional[pd.Series]:
+    """Single catchment-average hourly rain series (mm)."""
+    if df is None or df.empty:
+        return None
+    return df.mean(axis=1)
+
+
+def serving_series() -> Optional[pd.Series]:
+    """Catchment rain for inference: trailing days + the next 48 h forecast."""
+    try:
+        return catchment_mean(fetch_forecast())
+    except Exception as error:
+        logger.warning(f"Open-Meteo forecast fetch failed: {error}")
+        return None