Spaces:
Running on Zero
Running on Zero
| """ | |
| 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() |