diff --git a/src/ml/evaluate.py b/src/ml/evaluate.py index ccc9d03..f13140a 100644 --- a/src/ml/evaluate.py +++ b/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)) diff --git a/src/ml/features.py b/src/ml/features.py index e1032b6..871171f 100644 --- a/src/ml/features.py +++ b/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() diff --git a/src/ml/rain.py b/src/ml/rain.py new file mode 100644 index 0000000..7182404 --- /dev/null +++ b/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