scripts/evaluation/identify_flood_events.py

"""
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()