""" Identifies real historical high-flow events from discharge and water-level records, station by station. This is the ground truth against which any flood-prediction claim eventually has to be checked -- "beats a naive mean-baseline" and "would have flagged the real floods that actually happened" are different questions, and only the second one matters for this project's actual purpose. Threshold: each station's own historical 95th/99th percentile, not a shared absolute value across stations. These rivers have very different scales (median discharge ranges from under 1,000 to over 9,000 L/s across real stations in this project), so a fixed absolute threshold would flag "high flow" at a small tributary while missing real extremes at a larger one. A station-relative threshold puts every station's real flood-scale events on the same footing. Events, not days: a real flood spans several consecutive days above threshold, and counting each day as a separate "event" would both inflate the count and obscure how long real high-flow periods actually last -- duration matters for anything downstream (evacuation windows, infrastructure impact), not just peak magnitude. Usage: python -m scripts.identify_flood_events --data-root datasets """ import argparse from pathlib import Path import numpy as np import pandas as pd try: from src.graph.dynamic_features import build_discharge_timeseries, build_waterlevel_timeseries except ImportError: import sys sys.path.insert(0, str(Path(__file__).resolve().parent.parent)) from src.graph.dynamic_features import build_discharge_timeseries, build_waterlevel_timeseries def find_events(series: pd.Series, threshold: float, min_duration_days: int = 1) -> pd.DataFrame: """ Groups consecutive above-threshold days into discrete events for a single station's real time series. A gap of even one day below threshold ends the current event -- real recession rarely means the flow dropped and immediately rose again, so treating a brief dip as "still the same event" would understate how variable real hydrographs actually are. Returns one row per event: start date, end date, duration, peak value, peak date. Events shorter than min_duration_days are dropped -- a single anomalous reading (a possible data error, or a genuinely brief spike) is a different thing from a sustained high- flow event, and conflating them would dilute what "event" means for later evaluation. """ above = series >= threshold if not above.any(): return pd.DataFrame(columns=["start_date", "end_date", "duration_days", "peak_value", "peak_date"]) # Assign a group id that increments each time the series drops # below threshold -- consecutive True values share a group. group_id = (~above).cumsum() events = [] for gid, group in series[above].groupby(group_id[above]): duration = len(group) if duration < min_duration_days: continue events.append({ "start_date": group.index.min(), "end_date": group.index.max(), "duration_days": duration, "peak_value": group.max(), "peak_date": group.idxmax(), }) return pd.DataFrame(events) def identify_station_events( wide: pd.DataFrame, percentile: float = 95.0, min_duration_days: int = 1, ) -> pd.DataFrame: """ Runs find_events for every station in a [date x station_code] wide table, using each station's own real historical percentile as its threshold. Stations with too little real data to compute a meaningful percentile (fewer than 30 real observations, an arbitrary but reasonable floor -- a percentile from a handful of readings isn't a real threshold) are skipped, not given a threshold computed from noise. """ all_events = [] for station_code in wide.columns: series = wide[station_code].dropna() if len(series) < 30: continue threshold = np.percentile(series, percentile) station_events = find_events(series, threshold, min_duration_days) if station_events.empty: continue station_events["station_code"] = station_code station_events["threshold"] = threshold station_events["threshold_percentile"] = percentile all_events.append(station_events) if not all_events: return pd.DataFrame(columns=["station_code", "start_date", "end_date", "duration_days", "peak_value", "peak_date", "threshold", "threshold_percentile"]) return pd.concat(all_events, ignore_index=True).sort_values(["station_code", "start_date"]) def main() -> None: parser = argparse.ArgumentParser(description="Identify real historical high-flow events") parser.add_argument("--data-root", type=Path, default=Path("datasets")) parser.add_argument("--stations-path", type=Path, default=None) parser.add_argument("--start-date", type=str, default="2013-01-01") parser.add_argument("--end-date", type=str, default="2026-12-31") parser.add_argument("--percentile", type=float, default=95.0, help="Station-relative percentile defining the threshold for a high-flow day.") parser.add_argument("--min-duration-days", type=int, default=1) parser.add_argument("--output-dir", type=Path, default=None) args = parser.parse_args() stations_path = args.stations_path or (args.data_root / "station_elevations.csv") nodes_df = pd.read_csv(stations_path) date_range = (args.start_date, args.end_date) output_dir = args.output_dir or args.data_root discharge_wide = build_discharge_timeseries(nodes_df, args.data_root / "hydrometric", date_range) discharge_events = identify_station_events(discharge_wide, args.percentile, args.min_duration_days) print(f"Discharge: {len(discharge_events)} real event(s) across " f"{discharge_events['station_code'].nunique() if not discharge_events.empty else 0} station(s)") if not discharge_events.empty: print(discharge_events.to_string(index=False)) discharge_out = output_dir / "flood_events_discharge.csv" discharge_events.to_csv(discharge_out, index=False) print(f"Saved to {discharge_out}\n") waterlevel_wide = build_waterlevel_timeseries(nodes_df, args.data_root / "hydrometric", date_range) waterlevel_events = identify_station_events(waterlevel_wide, args.percentile, args.min_duration_days) print(f"Water level: {len(waterlevel_events)} real event(s) across " f"{waterlevel_events['station_code'].nunique() if not waterlevel_events.empty else 0} station(s)") if not waterlevel_events.empty: print(waterlevel_events.to_string(index=False)) waterlevel_out = output_dir / "flood_events_waterlevel.csv" waterlevel_events.to_csv(waterlevel_out, index=False) print(f"Saved to {waterlevel_out}") if __name__ == "__main__": main()