diff --git a/docs/FLOOD_FORECASTING.md b/docs/FLOOD_FORECASTING.md index 6de304d..5c38013 100644 --- a/docs/FLOOD_FORECASTING.md +++ b/docs/FLOOD_FORECASTING.md @@ -584,6 +584,34 @@ span. `src/ml/hii_rain.py` builds the catchment mean and a 24 h-sum bias/MAE/correlation, so the two sources' relationship is on record by the time the 2027 fold (train ≤ 2027-04-30, test Jun–Nov 2027) can test it. +### 2026-10-18: forecast rain was missing at serving time — fixed + +`build_features` cut the rain series to the gauge grid before taking the +forward sum, and the grid ends at the newest reading. So the row serving +predicts from always had `rain_fc24 = NaN`, although training and every +evaluation above had the real value (the Open-Meteo forecast, already fetched +with `forecast_days=2`). The model was answering as if no weather forecast +existed. The forward sum now runs over the grid extended by 24 h; earlier rows +are unchanged, so training is unaffected. + +Cost of the bug, measured with the harness: the `rise_rain_fc24_missing` +variant is `rise_rain` with `rain_fc24` blanked on the test rows (results in +`models/eval_2026-10-18_fc24.json`, five monsoon folds, `--from-cache`): + +| Station | Variant | Events | Warned early | Median lead | False alarms | MAE | +|---|---|---|---|---|---|---| +| P.1 | fixed (`rise_rain`) | 4 | 4 | 7.5 h | 0 | 0.087 | +| P.1 | as served | 4 | 4 | 7.0 h | 0 | 0.087 | +| P.103 | fixed | 7 | 7 | 9.0 h | 2 | 0.152 | +| P.103 | as served | 7 | 7 | 9.0 h | 2 | 0.153 | +| P.21 | fixed | 10 | 8 | 2.5 h | 5 | 0.215 | +| P.21 | as served | 10 | 7 | 2.0 h | 10 | 0.232 | + +P.1 and P.103 lost little; P.21 doubled its false alarms. The two very long +"as served" leads (P.21 2025-08-26 +69 h, P.103 2025-10-03 +35 h) are alarm +episodes that happened to precede an event, not real early warnings. Since the +fix, the served features match what this section evaluates. + ## 6. Deployment ### API diff --git a/models/eval_2026-10-18_fc24.json b/models/eval_2026-10-18_fc24.json new file mode 100644 index 0000000..e6f487d --- /dev/null +++ b/models/eval_2026-10-18_fc24.json @@ -0,0 +1,737 @@ +[ + { + "station": "P.1", + "warn_thr": 3.7, + "folds": [ + { + "year": 2021, + "n_train": 20024, + "n_test": 4392, + "events": [], + "variants": { + "rise_rain": { + "mae": 0.07375926701460789, + "mae_above_2p5": null, + "brier_warn": 0.0, + "events": [], + "false_alarm_episodes": 0 + }, + "rise_rain_fc24_missing": { + "mae": 0.07349391292098781, + "mae_above_2p5": null, + "brier_warn": 0.0, + "events": [], + "false_alarm_episodes": 0 + } + } + }, + { + "year": 2022, + "n_train": 28784, + "n_test": 4392, + "events": [ + { + "crossing": "2022-10-02T19:00:00", + "peak_ts": "2022-10-03T15:00:00", + "peak_level": 4.65 + } + ], + "variants": { + "rise_rain": { + "mae": 0.08452847354970178, + "mae_above_2p5": 0.25082252888260664, + "brier_warn": 0.004300908725927739, + "events": [ + { + "crossing": "2022-10-02T19:00:00", + "lead_h": 5.0, + "peak_level": 4.65, + "peak_pred_24h_before": 3.8173954245046406 + } + ], + "false_alarm_episodes": 0 + }, + "rise_rain_fc24_missing": { + "mae": 0.08511964468019045, + "mae_above_2p5": 0.2527703432441232, + "brier_warn": 0.004315341268867166, + "events": [ + { + "crossing": "2022-10-02T19:00:00", + "lead_h": 5.0, + "peak_level": 4.65, + "peak_pred_24h_before": 3.8102898827310137 + } + ], + "false_alarm_episodes": 0 + } + } + }, + { + "year": 2023, + "n_train": 37539, + "n_test": 4392, + "events": [], + "variants": { + "rise_rain": { + "mae": 0.07696637416933236, + "mae_above_2p5": 0.10643655855133666, + "brier_warn": 1.919860722404625e-15, + "events": [], + "false_alarm_episodes": 0 + }, + "rise_rain_fc24_missing": { + "mae": 0.07281928830426289, + "mae_above_2p5": 0.12255431002563086, + "brier_warn": 1.859416740909973e-15, + "events": [], + "false_alarm_episodes": 0 + } + } + }, + { + "year": 2024, + "n_train": 46323, + "n_test": 4392, + "events": [ + { + "crossing": "2024-09-24T17:00:00", + "peak_ts": "2024-09-26T02:00:00", + "peak_level": 4.93 + }, + { + "crossing": "2024-10-03T09:00:00", + "peak_ts": "2024-10-05T12:00:00", + "peak_level": 5.3 + } + ], + "variants": { + "rise_rain": { + "mae": 0.0890019157543212, + "mae_above_2p5": 0.2093245927883516, + "brier_warn": 0.005826169840715695, + "events": [ + { + "crossing": "2024-09-24T17:00:00", + "lead_h": 10.0, + "peak_level": 4.93, + "peak_pred_24h_before": 4.817519939833057 + }, + { + "crossing": "2024-10-03T09:00:00", + "lead_h": 21.0, + "peak_level": 5.3, + "peak_pred_24h_before": 5.546588884631041 + } + ], + "false_alarm_episodes": 0 + }, + "rise_rain_fc24_missing": { + "mae": 0.0896559193295349, + "mae_above_2p5": 0.2184968282693981, + "brier_warn": 0.00611626529479383, + "events": [ + { + "crossing": "2024-09-24T17:00:00", + "lead_h": 9.0, + "peak_level": 4.93, + "peak_pred_24h_before": 4.7051072221287855 + }, + { + "crossing": "2024-10-03T09:00:00", + "lead_h": 21.0, + "peak_level": 5.3, + "peak_pred_24h_before": 5.563518687086929 + } + ], + "false_alarm_episodes": 0 + } + } + }, + { + "year": 2025, + "n_train": 55083, + "n_test": 4392, + "events": [ + { + "crossing": "2025-09-27T18:00:00", + "peak_ts": "2025-09-27T22:00:00", + "peak_level": 3.93 + } + ], + "variants": { + "rise_rain": { + "mae": 0.11202835461695825, + "mae_above_2p5": 0.2456782496767171, + "brier_warn": 0.005178356039745929, + "events": [ + { + "crossing": "2025-09-27T18:00:00", + "lead_h": 2.0, + "peak_level": 3.93, + "peak_pred_24h_before": 3.23 + } + ], + "false_alarm_episodes": 0 + }, + "rise_rain_fc24_missing": { + "mae": 0.11232147044143151, + "mae_above_2p5": 0.2457513740974484, + "brier_warn": 0.004988582804599792, + "events": [ + { + "crossing": "2025-09-27T18:00:00", + "lead_h": 4.0, + "peak_level": 3.93, + "peak_pred_24h_before": 3.3264027753671597 + } + ], + "false_alarm_episodes": 0 + } + } + } + ] + }, + { + "station": "P.103", + "warn_thr": 5.95, + "folds": [ + { + "year": 2021, + "n_train": 20009, + "n_test": 4392, + "events": [], + "variants": { + "rise_rain": { + "mae": 0.14438683486790602, + "mae_above_2p5": null, + "brier_warn": 0.0, + "events": [], + "false_alarm_episodes": 0 + }, + "rise_rain_fc24_missing": { + "mae": 0.14435605468178428, + "mae_above_2p5": null, + "brier_warn": 0.0, + "events": [], + "false_alarm_episodes": 0 + } + } + }, + { + "year": 2022, + "n_train": 28769, + "n_test": 4392, + "events": [ + { + "crossing": "2022-08-14T04:00:00", + "peak_ts": "2022-08-14T08:00:00", + "peak_level": 6.09 + }, + { + "crossing": "2022-10-02T16:00:00", + "peak_ts": "2022-10-03T16:00:00", + "peak_level": 7.54 + } + ], + "variants": { + "rise_rain": { + "mae": 0.14392334606606189, + "mae_above_2p5": 0.41075382386412596, + "brier_warn": 0.00764679988213082, + "events": [ + { + "crossing": "2022-08-14T04:00:00", + "lead_h": 6.0, + "peak_level": 6.09, + "peak_pred_24h_before": 4.825270553204425 + }, + { + "crossing": "2022-10-02T16:00:00", + "lead_h": 9.0, + "peak_level": 7.54, + "peak_pred_24h_before": 7.194064117976157 + } + ], + "false_alarm_episodes": 0 + }, + "rise_rain_fc24_missing": { + "mae": 0.1438805977972288, + "mae_above_2p5": 0.41655333933205874, + "brier_warn": 0.007670614570637225, + "events": [ + { + "crossing": "2022-08-14T04:00:00", + "lead_h": 5.0, + "peak_level": 6.09, + "peak_pred_24h_before": 4.822220911656332 + }, + { + "crossing": "2022-10-02T16:00:00", + "lead_h": 9.0, + "peak_level": 7.54, + "peak_pred_24h_before": 7.187888665288963 + } + ], + "false_alarm_episodes": 0 + } + } + }, + { + "year": 2023, + "n_train": 37524, + "n_test": 4392, + "events": [], + "variants": { + "rise_rain": { + "mae": 0.16195301512514512, + "mae_above_2p5": null, + "brier_warn": 3.0208441147560167e-07, + "events": [], + "false_alarm_episodes": 0 + }, + "rise_rain_fc24_missing": { + "mae": 0.15117839869986002, + "mae_above_2p5": null, + "brier_warn": 6.508469941612771e-08, + "events": [], + "false_alarm_episodes": 0 + } + } + }, + { + "year": 2024, + "n_train": 46308, + "n_test": 4058, + "events": [ + { + "crossing": "2024-09-24T10:00:00", + "peak_ts": "2024-09-26T00:00:00", + "peak_level": 8.27 + }, + { + "crossing": "2024-09-30T03:00:00", + "peak_ts": "2024-09-30T06:00:00", + "peak_level": 5.99 + }, + { + "crossing": "2024-10-03T06:00:00", + "peak_ts": "2024-10-05T07:00:00", + "peak_level": 9.93 + } + ], + "variants": { + "rise_rain": { + "mae": 0.13901842325128816, + "mae_above_2p5": 0.2887416310966684, + "brier_warn": 0.0072040857753137046, + "events": [ + { + "crossing": "2024-09-24T10:00:00", + "lead_h": 19.0, + "peak_level": 8.27, + "peak_pred_24h_before": 8.212734363860193 + }, + { + "crossing": "2024-09-30T03:00:00", + "lead_h": 9.0, + "peak_level": 5.99, + "peak_pred_24h_before": 5.47565489914091 + }, + { + "crossing": "2024-10-03T06:00:00", + "lead_h": 55.0, + "peak_level": 9.93, + "peak_pred_24h_before": 8.780278464915938 + } + ], + "false_alarm_episodes": 0 + }, + "rise_rain_fc24_missing": { + "mae": 0.15101681306989473, + "mae_above_2p5": 0.3221006670360109, + "brier_warn": 0.008559364941781224, + "events": [ + { + "crossing": "2024-09-24T10:00:00", + "lead_h": 18.0, + "peak_level": 8.27, + "peak_pred_24h_before": 7.997091483400909 + }, + { + "crossing": "2024-09-30T03:00:00", + "lead_h": 9.0, + "peak_level": 5.99, + "peak_pred_24h_before": 5.502287543201456 + }, + { + "crossing": "2024-10-03T06:00:00", + "lead_h": 62.0, + "peak_level": 9.93, + "peak_pred_24h_before": 9.103449390810866 + } + ], + "false_alarm_episodes": 0 + } + } + }, + { + "year": 2025, + "n_train": 54577, + "n_test": 4392, + "events": [ + { + "crossing": "2025-09-26T06:00:00", + "peak_ts": "2025-09-27T21:00:00", + "peak_level": 6.64 + }, + { + "crossing": "2025-10-03T06:00:00", + "peak_ts": "2025-10-03T12:00:00", + "peak_level": 6.14 + } + ], + "variants": { + "rise_rain": { + "mae": 0.17206685418009837, + "mae_above_2p5": 0.3840520085986099, + "brier_warn": 0.015824852891785323, + "events": [ + { + "crossing": "2025-09-26T06:00:00", + "lead_h": 6.0, + "peak_level": 6.64, + "peak_pred_24h_before": 5.73043665401637 + }, + { + "crossing": "2025-10-03T06:00:00", + "lead_h": 8.0, + "peak_level": 6.14, + "peak_pred_24h_before": 5.609919760117381 + } + ], + "false_alarm_episodes": 2 + }, + "rise_rain_fc24_missing": { + "mae": 0.17368174393168262, + "mae_above_2p5": 0.3756076666801506, + "brier_warn": 0.014794520234147523, + "events": [ + { + "crossing": "2025-09-26T06:00:00", + "lead_h": 5.0, + "peak_level": 6.64, + "peak_pred_24h_before": 5.8809952440025 + }, + { + "crossing": "2025-10-03T06:00:00", + "lead_h": 35.0, + "peak_level": 6.14, + "peak_pred_24h_before": 5.6215400644356235 + } + ], + "false_alarm_episodes": 2 + } + } + } + ] + }, + { + "station": "P.21", + "warn_thr": 3.2, + "folds": [ + { + "year": 2021, + "n_train": 19099, + "n_test": 4392, + "events": [], + "variants": { + "rise_rain": { + "mae": 0.31976019104736464, + "mae_above_2p5": 0.49280187683998045, + "brier_warn": 4.599376306869306e-06, + "events": [], + "false_alarm_episodes": 0 + }, + "rise_rain_fc24_missing": { + "mae": 0.3238133154815793, + "mae_above_2p5": 0.49187453544981, + "brier_warn": 4.609640062881583e-06, + "events": [], + "false_alarm_episodes": 0 + } + } + }, + { + "year": 2022, + "n_train": 27860, + "n_test": 4392, + "events": [ + { + "crossing": "2022-09-12T05:00:00", + "peak_ts": "2022-09-12T05:00:00", + "peak_level": 3.21 + }, + { + "crossing": "2022-10-02T04:00:00", + "peak_ts": "2022-10-03T21:00:00", + "peak_level": 3.92 + } + ], + "variants": { + "rise_rain": { + "mae": 0.188862265489962, + "mae_above_2p5": 0.22738565836625405, + "brier_warn": 0.009674093199871036, + "events": [ + { + "crossing": "2022-09-12T05:00:00", + "lead_h": 0.0, + "peak_level": 3.21, + "peak_pred_24h_before": 2.49615362708868 + }, + { + "crossing": "2022-10-02T04:00:00", + "lead_h": 1.0, + "peak_level": 3.92, + "peak_pred_24h_before": 3.907087039152045 + } + ], + "false_alarm_episodes": 0 + }, + "rise_rain_fc24_missing": { + "mae": 0.19791868293297726, + "mae_above_2p5": 0.2739869896442161, + "brier_warn": 0.0103072884296978, + "events": [ + { + "crossing": "2022-09-12T05:00:00", + "lead_h": 0.0, + "peak_level": 3.21, + "peak_pred_24h_before": 2.4722245571508683 + }, + { + "crossing": "2022-10-02T04:00:00", + "lead_h": 0.0, + "peak_level": 3.92, + "peak_pred_24h_before": 3.845235629130954 + } + ], + "false_alarm_episodes": 0 + } + } + }, + { + "year": 2023, + "n_train": 36614, + "n_test": 4392, + "events": [ + { + "crossing": "2023-10-10T06:00:00", + "peak_ts": "2023-10-10T06:00:00", + "peak_level": 3.2 + } + ], + "variants": { + "rise_rain": { + "mae": 0.17283257678933697, + "mae_above_2p5": 0.39237224972369483, + "brier_warn": 0.012908395601661503, + "events": [ + { + "crossing": "2023-10-10T06:00:00", + "lead_h": 15.0, + "peak_level": 3.2, + "peak_pred_24h_before": 1.8019095357375217 + } + ], + "false_alarm_episodes": 3 + }, + "rise_rain_fc24_missing": { + "mae": 0.1953027570197896, + "mae_above_2p5": 0.37345466270918226, + "brier_warn": 0.005336247430915456, + "events": [ + { + "crossing": "2023-10-10T06:00:00", + "lead_h": 14.0, + "peak_level": 3.2, + "peak_pred_24h_before": 1.8796882667592647 + } + ], + "false_alarm_episodes": 3 + } + } + }, + { + "year": 2024, + "n_train": 45398, + "n_test": 4385, + "events": [ + { + "crossing": "2024-09-01T03:00:00", + "peak_ts": "2024-09-01T08:00:00", + "peak_level": 4.3 + }, + { + "crossing": "2024-09-23T06:00:00", + "peak_ts": "2024-09-23T15:00:00", + "peak_level": 4.68 + }, + { + "crossing": "2024-10-03T05:00:00", + "peak_ts": "2024-10-04T20:00:00", + "peak_level": 5.21 + } + ], + "variants": { + "rise_rain": { + "mae": 0.20733600140331027, + "mae_above_2p5": 0.4173432645880738, + "brier_warn": 0.013730524899775963, + "events": [ + { + "crossing": "2024-09-01T03:00:00", + "lead_h": 0.0, + "peak_level": 4.3, + "peak_pred_24h_before": 1.7973215452596953 + }, + { + "crossing": "2024-09-23T06:00:00", + "lead_h": 2.0, + "peak_level": 4.68, + "peak_pred_24h_before": 2.670881094688805 + }, + { + "crossing": "2024-10-03T05:00:00", + "lead_h": 33.0, + "peak_level": 5.21, + "peak_pred_24h_before": 4.999142187693036 + } + ], + "false_alarm_episodes": 0 + }, + "rise_rain_fc24_missing": { + "mae": 0.2222046566812092, + "mae_above_2p5": 0.43423994130380844, + "brier_warn": 0.01611358822165686, + "events": [ + { + "crossing": "2024-09-01T03:00:00", + "lead_h": 0.0, + "peak_level": 4.3, + "peak_pred_24h_before": 1.7847018175539242 + }, + { + "crossing": "2024-09-23T06:00:00", + "lead_h": 1.0, + "peak_level": 4.68, + "peak_pred_24h_before": 2.640820102168781 + }, + { + "crossing": "2024-10-03T05:00:00", + "lead_h": 32.0, + "peak_level": 5.21, + "peak_pred_24h_before": 4.749354278189956 + } + ], + "false_alarm_episodes": 0 + } + } + }, + { + "year": 2025, + "n_train": 54150, + "n_test": 4392, + "events": [ + { + "crossing": "2025-07-23T15:00:00", + "peak_ts": "2025-07-24T04:00:00", + "peak_level": 3.63 + }, + { + "crossing": "2025-08-26T19:00:00", + "peak_ts": "2025-08-26T23:00:00", + "peak_level": 3.48 + }, + { + "crossing": "2025-09-12T21:00:00", + "peak_ts": "2025-09-12T22:00:00", + "peak_level": 3.21 + }, + { + "crossing": "2025-09-22T20:00:00", + "peak_ts": "2025-09-22T21:00:00", + "peak_level": 3.27 + } + ], + "variants": { + "rise_rain": { + "mae": 0.1851599197920176, + "mae_above_2p5": 0.29411810271343447, + "brier_warn": 0.021849546400589632, + "events": [ + { + "crossing": "2025-07-23T15:00:00", + "lead_h": 3.0, + "peak_level": 3.63, + "peak_pred_24h_before": 2.4024371626847 + }, + { + "crossing": "2025-08-26T19:00:00", + "lead_h": 4.0, + "peak_level": 3.48, + "peak_pred_24h_before": 2.492383184457359 + }, + { + "crossing": "2025-09-12T21:00:00", + "lead_h": 18.0, + "peak_level": 3.21, + "peak_pred_24h_before": 2.5690976144688706 + }, + { + "crossing": "2025-09-22T20:00:00", + "lead_h": 1.0, + "peak_level": 3.27, + "peak_pred_24h_before": 2.382265050508778 + } + ], + "false_alarm_episodes": 2 + }, + "rise_rain_fc24_missing": { + "mae": 0.21936157981899435, + "mae_above_2p5": 0.3121582241464527, + "brier_warn": 0.025980409031233827, + "events": [ + { + "crossing": "2025-07-23T15:00:00", + "lead_h": 2.0, + "peak_level": 3.63, + "peak_pred_24h_before": 2.096100137296138 + }, + { + "crossing": "2025-08-26T19:00:00", + "lead_h": 69.0, + "peak_level": 3.48, + "peak_pred_24h_before": 2.3678295230657067 + }, + { + "crossing": "2025-09-12T21:00:00", + "lead_h": 19.0, + "peak_level": 3.21, + "peak_pred_24h_before": 2.783241642058761 + }, + { + "crossing": "2025-09-22T20:00:00", + "lead_h": 2.0, + "peak_level": 3.27, + "peak_pred_24h_before": 2.279238629032059 + } + ], + "false_alarm_episodes": 7 + } + } + } + ] + } +] \ No newline at end of file diff --git a/src/ml/data.py b/src/ml/data.py index 906ba66..1ce3b77 100644 --- a/src/ml/data.py +++ b/src/ml/data.py @@ -146,15 +146,22 @@ def _fetch_hii_levels( with engine.connect() as connection: df = pd.read_sql(text(query), connection, params=params) + return _prepare_hii_levels(df) + + +def _prepare_hii_levels(df: pd.DataFrame) -> pd.DataFrame: + """Hourly mirror rows. Several readings in one hour: the latest wins + (stable sort on the raw time, so the choice never depends on row order).""" df = df.dropna(subset=["wl_msl"]) if df.empty: return df - df["timestamp"] = pd.to_datetime(df["timestamp"]).dt.floor("h") + df = df.copy() + df["timestamp"] = pd.to_datetime(df["timestamp"]) + df = df.sort_values("timestamp", kind="stable") + df["timestamp"] = df["timestamp"].dt.floor("h") df["wl_msl"] = pd.to_numeric(df["wl_msl"], errors="coerce") df["discharge"] = pd.to_numeric(df["discharge"], errors="coerce") - df = df.sort_values("timestamp").drop_duplicates( - subset=["station_code", "timestamp"], keep="last" - ) + df = df.drop_duplicates(subset=["station_code", "timestamp"], keep="last") return df diff --git a/src/ml/evaluate.py b/src/ml/evaluate.py index 0ae8e0c..bc5dfdb 100644 --- a/src/ml/evaluate.py +++ b/src/ml/evaluate.py @@ -72,6 +72,7 @@ class Variant: use_dam: bool = False, use_fc48: bool = False, qsigma: bool = False, + fc24_missing_at_serve: bool = False, ): self.name = name self.target = target # 'abs' or 'rise' @@ -84,6 +85,10 @@ class Variant: # 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 + # Test rows get rain_fc24 = NaN: what production served until the + # 2026-10 fix (the rain series was cut to the gauge grid before the + # forward sum, so the newest row never had forecast rain). + self.fc24_missing_at_serve = fc24_missing_at_serve def fit_predict(self, X_tr, y_abs_tr, X_te) -> Tuple[np.ndarray, np.ndarray]: if not self.use_rain: @@ -108,6 +113,9 @@ class Variant: X_te = X_te.drop(columns=drop) elif "rain_fc48" not in X_tr.columns: raise ValueError(f"{self.name} requires the rain series") + if self.fc24_missing_at_serve and "rain_fc24" in X_te.columns: + X_te = X_te.copy() + X_te["rain_fc24"] = np.nan 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 @@ -164,6 +172,14 @@ VARIANTS: Dict[str, Variant] = { "rise_rain_qsigma": Variant( "rise_rain_qsigma", target="rise", use_rain=True, qsigma=True ), + # 2026-10-18: the deployed rise_rain model as it was actually served + # before the rain_fc24 serving fix (forecast rain missing at test time). + "rise_rain_fc24_missing": Variant( + "rise_rain_fc24_missing", + target="rise", + use_rain=True, + fc24_missing_at_serve=True, + ), } # Dam variants are opt-in by name: they require dam columns that only exist @@ -176,6 +192,7 @@ DEFAULT_VARIANTS = [ if not v.use_dam and not v.use_fc48 and not v.qsigma + and not v.fc24_missing_at_serve and not (v.quantile and v.use_rain) ] diff --git a/src/ml/features.py b/src/ml/features.py index 3b9dd3b..603a300 100644 --- a/src/ml/features.py +++ b/src/ml/features.py @@ -209,6 +209,7 @@ def _hours_since_observed(mask_col: pd.Series) -> pd.Series: RAIN_FEATURES = ("rain_6h", "rain_24h", "rain_72h", "rain_fc24") +FORWARD_RAIN_H = 24 # rain_fc24 window DAM_FEATURES = ("dam_storage_pct", "dam_storage_pct_d3", "dam_inflow", "dam_outflow") # Stations hydrologically downstream of the Mae Ngat confluence (Ping mainstem @@ -297,9 +298,25 @@ def build_features( 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 + # forward sum over (t, t+24]: shift(-1) starts the window at t+1. + # Computed on the grid extended by the 24 h after its last row: the + # last row is the one serving predicts from, and cutting the rain + # series to the grid made its forecast rain always NaN (unlike every + # training row since 2021). Earlier rows keep their values. + future = ( + pd.date_range(idx[-1], periods=FORWARD_RAIN_H + 1, freq="h")[1:] + if len(idx) + else idx[:0] + ) + r_ext = rain.reindex(idx.append(future)) cols["rain_fc24"] = ( - r.shift(-1).iloc[::-1].rolling(24, min_periods=1).sum().iloc[::-1] + r_ext.shift(-1) + .iloc[::-1] + .rolling(FORWARD_RAIN_H, min_periods=1) + .sum() + .iloc[::-1] + .iloc[: len(idx)] + .set_axis(idx) ) if dam is not None and station in DAM_STATIONS: diff --git a/tests/test_rain_serving.py b/tests/test_rain_serving.py new file mode 100644 index 0000000..55c0a36 --- /dev/null +++ b/tests/test_rain_serving.py @@ -0,0 +1,64 @@ +"""rain_fc24 at the newest grid row (the one serving predicts from).""" + +import numpy as np +import pandas as pd + +from src.ml import features +from src.ml.data import _prepare_hii_levels + + +def _grid(hours=48): + idx = pd.date_range("2026-09-01", periods=hours, freq="h") + df = pd.DataFrame( + { + "timestamp": idx, + "station_code": "P.1", + "water_level": np.linspace(1.0, 2.0, hours), + "discharge": np.nan, + } + ) + return features.make_hourly_grid(df) + + +def test_forward_rain_reaches_past_the_last_reading(): + grid = _grid() + idx = grid.observed.index + rain = pd.Series( + np.arange(1, len(idx) + 49, dtype=float), + index=pd.date_range(idx[0], periods=len(idx) + 48, freq="h"), + ) + X = features.build_features(grid, "P.1", rain=rain) + last = idx[-1] + expected = rain.loc[last + pd.Timedelta(hours=1) : last + pd.Timedelta(hours=24)] + assert X.loc[last, "rain_fc24"] == expected.sum() + # rows whose window lies inside the grid are what they always were + old = rain.reindex(idx).shift(-1).iloc[::-1].rolling(24, min_periods=1).sum() + inner = idx[:-24] + np.testing.assert_allclose(X.loc[inner, "rain_fc24"], old.iloc[::-1].loc[inner]) + + +def test_forward_rain_nan_without_forecast_hours(): + grid = _grid() + idx = grid.observed.index + rain = pd.Series(1.0, index=idx) # nothing after the last reading + X = features.build_features(grid, "P.1", rain=rain) + assert np.isnan(X.loc[idx[-1], "rain_fc24"]) + assert X.loc[idx[-2], "rain_fc24"] == 1.0 + + +def test_hii_rows_latest_reading_in_an_hour_wins(): + raw = pd.DataFrame( + { + "timestamp": pd.to_datetime( + ["2026-10-01 10:40", "2026-10-01 10:10", "2026-10-01 11:00"] + ), + "station_code": "P.1", + "wl_msl": [301.5, 301.2, 301.7], + "discharge": [10.0, 9.0, 11.0], + } + ) + out = _prepare_hii_levels(raw) + assert list(out["wl_msl"]) == [301.5, 301.7] + assert list(out["timestamp"]) == list( + pd.to_datetime(["2026-10-01 10:00", "2026-10-01 11:00"]) + )