"""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"])), }