README.md

metadata

title: La Eure / La Risle Hydrometric GNN
emoji: 📈
colorFrom: indigo
colorTo: indigo
sdk: gradio
sdk_version: 5.47.0
app_file: src/gradio_app.py
pinned: true

La Eure / La Risle Hydrometric GNN

A physics-informed graph neural network that predicts streamflow (discharge, water level) at gauged and ungauged points along two Normandy rivers, La Eure and La Risle — each modeled as a real reach-based network (confluences, braided splits/rejoins, ~4,500 nodes per basin including virtual infill points), not a single chain of gauges, with covariates pulled from nine independent data sources.

The model itself is a spatiotemporal GNN trained with quantile regression (not a single point estimate) and four graph-tied physics losses, evaluated directly against real historical flood events rather than only against a naive baseline. See §4 for the architecture and real, current evaluation results.

flowchart LR
    hubeau["Hub'Eau<br/>discharge · water level<br/>catchment area"]
    ades["ADES<br/>groundwater levels"]
    era5["Copernicus ERA5<br/>climate reanalysis"]
    otd["Open Topo Data<br/>station elevation"]
    brgm["BRGM<br/>IDPR · BD Charm-50 geology"]
    bdtopo["IGN BD TOPO<br/>real reach topology + catchment polygons"]
    bdcav["Géorisques<br/>BDCavités (sinkholes)"]
    wc["ESA WorldCover<br/>landcover · NDVI"]
    mf["Météo-France<br/>real historical precipitation"]
    ecmwf["ECMWF (Open-Meteo)<br/>real forecast archive"]

    bdtopo --> brg["build_reach_graph.py<br/>real confluences, splits/rejoins,<br/>gauge snapping"]
    brg --> brgs["build_reach_graphs.py<br/>~4,500 nodes/basin"]

    hubeau --> nf
    ades --> nf
    era5 --> nf
    otd --> nf
    brgm --> nf
    bdcav --> nf
    wc --> nf
    brgs --> nf["node_features.py /<br/>enrich_reach_graph.py<br/>date-filtered 2013-2026"]

    bdtopo --> cc["compute_cumulative_catchment.py<br/>graph-wide catchment area"]
    cc --> nf

    nf --> pyg["build_pyg_graph<br/>x_static / x_dynamic split"]
    pyg --> phys["physics_losses.py<br/>confluence · split-rejoin ·<br/>routing · manning · water balance"]

    mf --> train
    ecmwf --> train
    phys --> train["train_spatiotemporal_gnn.py<br/>quantile regression + physics losses,<br/>100-node subgraph"]
    train --> model["model.pt<br/>trained checkpoint"]
    model --> eval["evaluate_flood_detection.py<br/>vs. real historical flood events"]

    pyg --> app["src/app.py<br/>Streamlit explorer +<br/>network validation view"]
    pyg --> testsuite["test_build_graph.py<br/>validation"]

---

1. Repository layout

PoC_v1/
├── scripts/                     # organized by purpose, not a flat list — see the
│   │                             #   category note below the tree
│   ├── data_acquisition/            # everything that talks to an external API/WFS/S3
│   │   ├── download_hubeau.py           # discharge + water level, Hub'Eau API v2
│   │   ├── download_elevation.py        # point elevations, Open Topo Data
│   │   ├── download_era5_sample.py      # ERA5 sanity-check pull (Jan 2020 only)
│   │   ├── download_era5_full.py        # ERA5 1960–2026, split instant/accum vars
│   │   ├── extract_era5.py              # unzips CDS API's zipped NetCDF output
│   │   ├── download_catchment.py        # Hub'Eau referentiel/sites -> surface_bv
│   │   ├── download_bdtopo_hydro.py     # IGN WFS -> tronçons, surfaces, catchments
│   │   ├── download_bdcavites.py        # Géorisques BDCavités (sinkhole/cavity inventory)
│   │   ├── download_bdcharm.py          # BRGM BD Charm-50 harmonized geology, per department
│   │   ├── fetch_landcover.py           # ESA WorldCover landcover class, real gauges
│   │   ├── fetch_worldcover_ndvi.py     # ESA WorldCover NDVI percentile composite
│   │   ├── fetch_ecmwf_forecast.py      # live ECMWF forecast, Open-Meteo (current conditions only)
│   │   ├── fetch_previous_runs.py       # historical ECMWF forecast archive, Open-Meteo Previous Runs API
│   │   └── fetch_meteofrance_precipitation.py  # real historical precipitation, Météo-France public API
│   │
│   ├── graph_building/              # everything that builds or enriches the real reach graph
│   │   ├── build_reach_graphs.py        # real reach-based topology, both basins
│   │   ├── enrich_reach_graph.py        # runs node_features.py against the reach graph
│   │   ├── compute_cumulative_catchment.py  # graph-wide catchment area from BD TOPO polygons
│   │   ├── compute_edge_width.py        # real channel width from BD TOPO polygons, per edge
│   │   ├── diagnose_confluences.py      # verify real vs. artifact confluences
│   │   └── build_dynamic_tensors.py     # genuine [n_nodes, T] tensors, wired into physics_losses.py
│   │
│   ├── analysis/                    # one-off diagnostics, not part of the main pipeline
│   │   ├── analyze_bdtopo_hydro.py      # centerline export + karst check
│   │   ├── run_bdtopo_checks.py         # karst + catchment cross-check, one shot
│   │   └── cross_check_catchments.py    # spatial join: station -> containing polygon
│   │
│   ├── training/
│   │   └── train_spatiotemporal_gnn.py  # main model: quantile regression + physics losses (§4)
│   │
│   └── evaluation/
│       ├── identify_flood_events.py     # real historical high-flow events — ground truth
│       └── evaluate_flood_detection.py  # precision/recall of the trained model vs. real historical floods
│
├── src/
│   ├── app.py                       # Streamlit river explorer + network validation view
│   ├── generate_plots.py            # batch plot generation across all loaders
│   ├── test_build_graph.py          # graph-construction test/validation suite
│   ├── extract_river_centerline.py  # digitizes a traced map image into a centerline
│   │
│   ├── data/
│   │   ├── loaders/
│   │   │   ├── base.py                  # BaseDataLoader — shared load()/get_metadata()
│   │   │   ├── hydrometric.py           # discharge & water level (Hub'Eau)
│   │   │   ├── ades.py                  # groundwater levels (ADES)
│   │   │   ├── safran.py                # ERA5 reanalysis, vectorized station interpolation
│   │   │   ├── idpr.py                  # infiltration/runoff tendency (BRGM)
│   │   │   ├── catchment.py             # per-station catchment area (Hub'Eau)
│   │   │   ├── bdtopo_hydro.py          # IGN BD TOPO hydrography (GeoJSON)
│   │   │   ├── shapefile.py             # watershed boundary polygon
│   │   │   └── station_elevations.py    # station coordinates + elevation
│   │   │
│   │   ├── river_graph.py           # basin assignment, elevation ordering, edges
│   │   ├── river_line.py            # straight-line interpolation between gauges
│   │   └── river_centerline.py      # real-centerline interpolation + gauge snapping
│   │
│   └── graph/
│       ├── build_graph.py           # PyG conversion: x_static/x_dynamic split, structural columns
│       ├── build_reach_graph.py     # real reach topology: confluences, splits/rejoins, MultiDiGraph
│       ├── node_features.py         # pulls every loader into one feature table (static, one row/node)
│       ├── dynamic_features.py      # genuine [n_nodes, T] series: discharge, groundwater, climate,
│       │                            #   precipitation, forecast (§4.4, §4.5)
│       ├── physics_losses.py        # confluence/split-rejoin/routing/manning/water-balance loss terms,
│       │                            #   NaN-masked (§2.5, §4.3)
│       ├── spatiotemporal_gnn.py    # the model — TypeAwareInputEncoder, RiverMessagePassing (spatial),
│       │                            #   per-node GRU (temporal), quantile decoder (§4.1)
│       ├── pinball_loss.py          # quantile regression loss, asymmetric and NaN-masked (§4.2)
│       └── subgraph_selection.py    # real-gauge + confluence/braid subgraph reduction for
│                                     #   tractable training (§4.6)
│
├── datasets/                    # not checked in; populated by the scripts above
│   ├── station_list.csv             # raw station roster (X, Y, names, INSEE, etc.)
│   ├── station_elevations.csv       # station_code, lat, lon, elevation_m
│   ├── idpr.csv
│   ├── catchment_area.csv
│   ├── edge_widths.csv              # real channel width per edge, from compute_edge_width.py
│   ├── ades/
│   ├── hydrometric/
│   ├── safran/
│   ├── bdtopo_hydro/
│   ├── bdcavites/
│   ├── bdcharm50/
│   ├── centerlines/
│   ├── meteofrance_precipitation/    # real historical precipitation, Météo-France (§3.12)
│   │   └── daily_precipitation.csv
│   ├── previous_runs_forecast.csv    # real ECMWF forecast archive, partial coverage (§3.13)
│   ├── flood_events_discharge.csv    # real historical high-flow events — ground truth (§4.7)
│   ├── flood_events_waterlevel.csv
│   ├── flood_detection_evaluation.csv  # confusion matrix by lead time, real held-out test data
│   ├── reach_graph/                 # {eure,risle}_{nodes,edges}.csv, _nodes_enriched.csv
│   └── spatiotemporal_gnn/
│       └── model.pt                 # trained model checkpoint

scripts/ talks to the outside world (APIs, WFS, S3) or trains/evaluates a model against real data; src/ doesn't — nothing under src/ makes a network call, and a script under src/ that wants one is a bug. Most of src/data/loaders/ predates the graph work — general-purpose readers/plotters for each dataset, with node_features.py stitching them together afterward rather than the other way around. scripts/ itself is organized into subfolders by purpose (data_acquisition/, graph_building/, analysis/, training/, evaluation/) — every category listed in the tree above needs its own __init__.py for python -m scripts.category.script_name to resolve correctly.

---

2. The graph

This is the part everything else in the repo exists to feed. Two graphs, one per river — H4xx… stations feed the La Eure graph, H6xx… feed La Risle — built with no edge between them, because there's no surface connection between the two basins to model.

The graph is now built from real reach topology, not a single ordered chain of gauges. build_reach_graph.py constructs it directly from BD TOPO's own tronçon-to-node linkage (lien_vers_noeud_hydrographique_ini/fin) — the NEXT_DOWN-equivalent approach — rather than inferring station order from position along a digitized line. That means real branching, real confluences, and real braided-channel structure fall directly out of the data instead of needing to be modeled separately.

2.1 Node types

Four kinds of node, not one:

Type	What it is	Column
Real gauge	one of the 27 hydrometric stations	is_gauged
Real confluence	a genuinely different, independently-sourced river joins	is_confluence
Split / rejoin	a channel divides and later recombines (braiding, an anabranch) — same water, no new mass	is_split_point / is_rejoin_point, paired via braid_id
Virtual (infill)	inserted along long confluence-free stretches so "predict at any point" has real spatial resolution	none of the above

A confluence requires more than a shared node with in-degree ≥ 2 — BD TOPO's fine tronçon segmentation produces plenty of same-river multi-inflow points with no real branching involved (confirmed against real data: incoming-edge distances as short as 4.6 m at some falsely-flagged "confluences"). The real test (find_real_confluences in build_reach_graph.py) requires (a) more than one distinct normalized river name among the incoming edges — river-name normalization strips articles, parenthetical qualifiers, and "bras de/du/d'" (arm-of) prefixes, since a named secondary channel of the same river ("Bras de la Charentonne") isn't a different river — and (b) that those branches don't trace back to a common upstream split within 15 km, which would mean it's a rejoin, not a confluence. Splits themselves need no such disambiguation: out-degree ≥ 2 is an unambiguous physical definition on its own, since a split by construction has exactly one thing flowing in.

Real branching topology also meant the underlying graph had to move from a plain DiGraph to a MultiDiGraph — two distinct tronçons directly connecting the same two hydrographic nodes (exactly the shape a short braid takes) is real data, not a collision, and a plain DiGraph was silently overwriting the second such edge's data on add_edge rather than keeping both. Confirmed as a real bug with real impact, not just a synthetic-test concern: fixing it recovered dozens of previously-invisible parallel edges per basin on the actual data.

2.2 Node and edge features

The feature set now spans several independent sources, each merged onto the node table by node_features.py's add_*_features functions. Every column lands in exactly one of four places once build_pyg_graph processes it:

flowchart TD
    raw["Enriched node table<br/>(node_features.py)"]

    raw --> struct{"structural /<br/>graph-role column?"}
    struct -->|"is_gauged, is_confluence,<br/>is_split_point, is_rejoin_point,<br/>braid_id, snap_distance_km"| structout["data.is_gauged, data.is_confluence, ...<br/>own Data attribute — never in x"]

    raw --> tgt{"target_* column?"}
    tgt -->|"target_discharge_m3s_*<br/>target_waterlevel_mm_*"| y["data.y<br/>never in x — label leakage otherwise"]

    raw --> feat{"real model input"}
    feat -->|"static: elevation_m, idpr_*,<br/>catchment_area_km2, landcover_*,<br/>geology_*, cavites distance/count"| xstatic["data.x_static"]
    feat -->|"dynamic: climate_*,<br/>avg_groundwater_*, ndvi_*<br/>(period-aggregate, not a real series yet)"| xdynamic["data.x_dynamic"]

    xstatic --> x["data.x — full combined tensor,<br/>z-scored"]
    xdynamic --> x

    edges["Edge table<br/>(build_reach_graph_tables)"] --> eattr{"numeric edge<br/>attribute?"}
    eattr -->|"distance_km,<br/>elevation_drop_m,<br/>verified_continuous,<br/>width_m"| edgeattr["data.edge_attr<br/>[n_edges, 4]"]
    eattr -->|"toponym, cleabs<br/>(diagnostic metadata)"| meta["not used by build_pyg_graph —<br/>stays in edges_df only"]

Node features:

Feature	Source	Coverage
latitude, longitude, elevation_m	station coords / real BD TOPO tronçon Z	every node
idpr_value, idpr_nearest_point_distance	BRGM IDPR	every node (spatial fallback for non-gauge codes)
catchment_area_km2	Hub'Eau, cumulative, real gauges only	27 stations
cumulative_catchment_area_km2	BD TOPO incremental polygons, summed upstream via real graph topology	graph-wide (~98% of nodes)
landcover_* (one-hot)	ESA WorldCover 10 m classification	real gauges only, for now
ndvi_p10, ndvi_p50, ndvi_p90	ESA WorldCover NDVI percentile composite	real gauges only, for now
geology_* (one-hot)	BRGM BD Charm-50, point-in-polygon	real gauges only, for now
distance_to_nearest_cavity_km, n_cavities_within_20km	Géorisques BDCavités, KD-tree + haversine	real gauges only, for now
avg_groundwater_level_m, avg_groundwater_depth_m, n_nearby_wells	ADES, radius-averaged, KD-tree + exact haversine	every node
climate_* (temp/wind/solar/precip/evap/snow/runoff)	ERA5, vectorized station interpolation	every node (needs safran_path)
{col}__was_missing	auto-generated	any feature column with real gaps

Edge features — four numeric attributes per edge, from build_reach_graph.py's build_reach_graph_tables (three original plus real channel width, added by compute_edge_width.py):

Feature	Meaning
distance_km	along-river distance between the two endpoint nodes
elevation_drop_m	elevation difference, upstream minus downstream — negated on the reverse edge when bidirectional=True
verified_continuous	False for any edge deliberately flagged via known_losing_reaches (the bétoire stretch — §3.7)
width_m (+ width_missing_flag)	real channel width, derived from BD TOPO polygon area ÷ tronçon length; NaN (not zero) where no real polygon match exists, with an explicit missingness flag rather than a silent fill — real coverage is partial, not universal

toponym and cleabs also live on the real edges table (the tronçon's river name and unique BD TOPO ID) but are diagnostic metadata, not model input — build_pyg_graph selects edge_attr columns by explicit name, so extra columns like these pass through harmlessly rather than needing to be stripped out first.

Structural columns never enter x. is_gauged, is_confluence, is_split_point, is_rejoin_point, snap_distance_km, braid_id describe node role, not a physical covariate — build_pyg_graph's auto-detection excludes them explicitly (confirmed as a real, not hypothetical, bug once: pandas treats bool as a numeric dtype, so without this exclusion these columns were being silently z-scored and fed to the model as if they were elevation or precipitation). They're still attached to the returned Data object as their own typed attributes, for masking supervised loss to gauged nodes and for the physics-loss index builders.

Landcover and geology are one-hot, not a raw class code. Both are nominal categories (10 = Tree cover, 50 = Built-up; a geological formation code), not an ordered quantity — leaving either as a raw integer would let auto-detection z-score it as if one category were numerically "more" than another, the same class of error as the structural-column bug, just subtler since these are meant to be real model input.

Targets are not features. target_discharge_m3s_mean/std/count and target_waterlevel_mm_mean/std/count exist on the enriched table but never enter x — they're pulled out into data.y separately, and attach only to real gauge rows (verified: gauge codes, BD TOPO hydrographic node IDs, and virtual-node marker strings occupy structurally distinct namespaces, so a left-merge on station_code can never mislabel a confluence or virtual node).

2.3 Static vs. dynamic features — and a real temporal pipeline

build_pyg_graph splits every feature by physical temporal nature:

• data.x_static / data.static_feature_names — genuinely time-invariant: elevation, IDPR, catchment area, landcover, geology, cavité proximity, coordinates.
• data.x_dynamic / data.dynamic_feature_names — physically time-varying quantities, still as a single period-aggregated number here (mean/sum over the whole date range, or a latest well reading) — this tensor is a static snapshot of dynamic-natured quantities, not a real series. data.x remains the full combined tensor unchanged; the split is additional, not a replacement.

A genuine [n_nodes, T] series exists separately, in src/graph/dynamic_features.py — build_discharge_timeseries, build_waterlevel_timeseries, build_groundwater_timeseries, build_climate_timeseries, and (added for model training, §4) build_precipitation_timeseries and build_forecast_tensor. Same loaders as everywhere else, no re-fetching; the difference is that these functions pivot to wide [date x station_code] form (or, for the forecast archive, a real per-lead-time lookup — see §3.13) instead of collapsing to one aggregate the way node_features.py's add_*_features do.

Two real challenges, not incidental engineering:

• Groundwater reports on wildly irregular schedules (confirmed: 13 different "latest dates" among 18 real wells within 20 km of one station). Each well is resampled to a common daily grid via forward-fill (a water table changes slowly — carrying the last known reading forward is standard practice, not an invented shortcut) before spatial averaging, not after — averaging raw irregular readings per exact calendar date is exactly what made the static version undercount real coverage by 5–10x before that was fixed (§3.3). The spatial neighbor-set per node is computed once, reused across every date — verified fast at real reach-graph scale (21.6s for 2,900 nodes × 14 years daily, real 272k-row ADES data).
• Discharge deliberately does not get forward-filled the way groundwater does — a missing daily reading stays missing, since discharge genuinely changes day to day and papering over a gap with yesterday's value would misrepresent it.

scripts/build_dynamic_tensors.py is the actual wiring: builds these tensors for a basin, saves them, and feeds discharge directly into routing_consistency_loss alongside build_routing_index — the real integration point, not just parallel unconnected pieces. train_spatiotemporal_gnn.py (§4) is the fuller integration point, wiring discharge, water level, groundwater, precipitation, and (partially) forecast data together into one training loop.

Climate is untested against real data — no real ERA5/safran_path files were available to validate build_climate_timeseries against in this project's development environment; the logic mirrors the already-tested discharge pivot directly, but verify the real output before trusting it.

2.4 Date-range filtering

build_node_features/enrich_reach_graph.py accept a date_range applied to every time-varying source (groundwater, climate, hydrometric targets) together, so all three describe the same period rather than each silently aggregating over its own full, differently-shaped history (ADES wells reporting from the 1970s to 2026 on wildly different schedules; ERA5 spanning 1960–2026; hydrometric records with their own per-station ranges entirely).

Default: 2013-01-01 to 2026-12-31 — computed, not guessed, via a brute-force interval-overlap check across all 8 discharge-gauged stations' real date ranges. This is the window that maximizes simultaneous station coverage: 6 of 8 stations, 8,923 real, quality-filtered observations (code_qualification >= 16, the same threshold HydrometricLoader itself applies — a naive raw count that skips this filter gives 13,084, which is what an earlier pass at this analysis originally reported before the discrepancy was traced and corrected). Two stations (H403301101: 1969–1985, H605022010: 1970–1980) are permanently excluded by any reasonable window — a ~35–40 year dead gap separates them from every other station's record, so including them would mean spanning six mostly-empty decades, not a genuine improvement.

2.5 Physics-informed loss terms (physics_losses.py)

Five constraints, each tied to real graph structure, not generic — four original, plus Manning's equation added alongside the model (§4.3):

Term	Constraint	Applies to
confluence_mass_balance_loss	Q_confluence ≈ sum(Q_upstream_branches) — new mass genuinely enters	is_confluence nodes
split_rejoin_conservation_loss	Q_split ≈ Q_rejoin — same water, no new mass	paired braid_id nodes
routing_consistency_loss	Q_downstream[t] ≈ Q_upstream[t - lag], lag from real distance_km/slope	every edge, real [n_nodes, T] via dynamic_features.py (§2.3)
manning_consistency_loss	Q ≈ (1/n)·A·R^(2/3)·S^(1/2), learned roughness n, with a bank-height relaxation gate so overbank flow isn't penalized against in-channel physics	every edge with real width data
water_balance_loss	P - ET - Q - ΔS ≈ 0 in volume terms, normalized to a relative (not absolute) imbalance	nodes with cumulative_catchment_area_km2 and real precipitation

Confluence and split/rejoin are deliberately different constraints, not one generic "conserve mass everywhere" rule — a model that only learned "sum the inflows" would get a split/rejoin wrong, since a rejoin's two branches together should equal the split's value, not add something new on top. All four original terms apply graph-wide, not just at the 27 labeled gauges — that's the actual mechanism by which sparse supervision generalizes to the ungauged nodes, not an incidental detail. ΔS (storage change) defaults to zero, a named steady-state approximation — this project has no direct basin-wide storage measurement, only sparse well levels, which aren't the same thing; likewise ET defaults to 0.0, since real evapotranspiration data exists (from the ERA5 climate work) but lives on a different grid than this graph's nodes, and mapping between them is unstarted work.

All five are NaN-masked, not just tolerant of complete data. Real ground- truth Q is ~93.5% NaN by construction (only real gauges with real observations ever have a value — confirmed against the real discharge tensor) — that's the normal shape of the data, not a rare edge case. The shared _mse helper every loss function uses previously computed a plain mean, so a single NaN anywhere in a residual silently poisoned the entire loss to NaN — confirmed as a real, not hypothetical, failure. _mse now masks NaN out before averaging (returning NaN only if truly nothing usable exists at all, which is a real "no data" signal worth keeping, not silently averaging to a misleading 0).

Masking the value isn't sufficient on its own, either — a second, separate real bug found during model training (§4.3): a masked-out NaN upstream input can still poison a gradient, since 0 × NaN = NaN in IEEE754 arithmetic even when the outer loss correctly excludes that position. The real fix, applied in water_balance_loss, is to replace NaN inputs with a safe finite placeholder before any arithmetic touches them, and apply the real valid mask only at the very end, for loss selection — not to rely on masking alone to keep NaN out of the graph.

This also means these functions are directly usable as a diagnostic against real historical data alone, independent of any trained model — e.g. "does real observed discharge at two connected gauges actually satisfy the routing physics" — a genuine, model-free sanity check on both the physics math and the graph topology, not just a training-time loss term.

2.6 Two graphs, not one

build_pyg_graphs_per_basin() returns {0: eure_graph, 1: risle_graph}, each with its own local 0..n-1 node indexing, rather than one merged Data object with two disconnected components. La Eure and La Risle are distinct hydrographic systems with nothing connecting them at the surface, and PyTorch Geometric's own batching (Batch.from_data_list) expects a list of separate small graphs — building two graphs from the start matches that convention directly.

2.7 What it looks like

The Streamlit explorer (src/app.py) has two views. "Explore" renders the original click-to-read interface over real course geometry. "Network validation" renders the full reach graph — confluences as diamonds, gauges as elevation-colored circles, every edge as one line trace regardless of edge count (verified fast at real scale: 0.29s to build a figure for ~2,900 edges) — specifically for visually confirming the topology looks like a real river network before trusting it as model input.

---

3. Datasets

Every dataset here has its own quirks, and in a couple of cases the quirks materially affect what the data means.

3.1 Station roster (station_list.csv, station_elevations.csv)

27 stations across the two basins, spanning three French departments — verified directly against the real roster: 12 in Eure (27), 10 in Eure-et-Loir (28, the Eure's southern tributaries near Chartres/Dreux — Voise, Drouette, and others), 2 in Orne (61). Split roughly by Hub'Eau code prefix (H4xx… for La Eure, H6xx… for La Risle — a heuristic based on observed codes, not a documented rule). Elevation comes from Open Topo Data's eudem25m endpoint (scripts/download_elevation.py), queried per station coordinate — a point lookup, not a raster, so there's no slope or catchment information hiding in it.

The row order in station_list.csv does not follow the river's course — verified directly, it jumps around in both latitude and elevation. Anything that needs upstream/downstream ordering has to derive it from elevation, latitude, or real centerline/graph position; never from file order.

Not every station in this list is actively gauged. Cross-referencing station names against the hydrometric data turned up three categories worth knowing about:

• Manual "observateur" stations — read by a person, not telemetered, so there's no digital time series to have. Four of these in the current roster.
• Partner-network ("SEBV") stations — operated outside the standard Hub'Eau telemetry network, likely need a different data source entirely if you want their readings. H431021010 is one of these.
• Everything else with no data is unexplained from the name alone and worth a direct check on Hub'Eau's site before assuming it's just a gap.

14 of the 27 stations have neither real discharge nor real water-level data at all. These aren't excluded from the graph — they still get real static attributes (landcover, geology, catchment area) and, since §4.4, a real precipitation-based discharge estimate as an input feature (never a target), plus whatever signal reaches them through message passing from nearby gauged nodes.

3.2 Hydrometric data (hydrometric/, via scripts/download_hubeau.py)

Discharge and water level from Hub'Eau's obs_elab endpoint. One thing that trips up a naive read: both files contain a mix of grandeur_hydro_elab codes, not just the variable implied by the filename. discharge_observations.csv has HIXnJ/HIXM (water-level codes) sitting right alongside QmnJ (discharge) rows for the same stations — HydrometricLoader filters each file down to its intended grandeur code explicitly (QmnJ from the discharge file, HIXnJ from the water-level file) rather than trusting the filename.

HydrometricLoader also filters on Hub'Eau's own code_qualification field, keeping only >= 16 (their "acceptable"/"good" threshold) and dropping lower- quality/provisional readings. This is real and meaningful, not a rounding detail — traced directly against the raw discharge file for the 2013–2026 window (§2.4): 13,084 raw QmnJ rows in range, of which 4,022 have code_qualification == 12 (below the threshold) and get correctly excluded, leaving 8,923. Any manual read of the raw CSVs that skips this filter will overcount real usable observations by close to a third.

Only 8 of the 27 stations have any QmnJ (daily mean discharge) rows at all. Several others report water level only. This isn't evenly distributed and matters a lot for anything downstream that assumes "gauged" means "has both variables" — it usually doesn't.

3.3 Groundwater (ades/)

ADES piezometer data: 113 wells in the watershed extract, 102 of them with actual level readings, going back as far as 1967 for some wells and to 2026-07-05 for the most recent reading at time of writing. ADESLoader.load() merges the levels file against the stations file on code_bss and renames x/y to lon/lat — those columns are already in degrees in this dataset, not a projected CRS, so no reprojection happens or is needed.

node_features.py's add_groundwater_features does not use ADESLoader.aggregate_to_stations — that method loops per station and does a full haversine .apply() over the entire groundwater dataframe for each one. At 27 stations against ~272k readings that's slow but tolerable; at the reach graph's ~4,500 nodes it's over a billion row-wise Python calls, confirmed as a genuine, not hypothetical, multi-hour hang. The fix (reduce to each well's latest reading first, then a KD-tree coarse prefilter + exact haversine on the small candidate set) turned out to also fix a real accuracy bug: the old method required wells to share the exact same reporting date before averaging, but real wells report on wildly different schedules (18 real wells within 20 km of one station spanned 13 different "latest dates," one from 1972) — silently discarding most real coverage every time.

Groundwater is used as a station-level input covariate, not as a graph edge. Well proximity alone isn't sufficient grounds for a subsurface/karst connectivity edge — that would need either correlated well hydrographs over time or a shared BDLISA aquifer-unit code (groundwater_stations.csv has a codes_bdlisa column available for exactly this kind of check; still unused).

Well coverage is not uniform across the two basins. The Eure's southern reach (south of roughly 48.68°N, toward Chartres) has essentially zero wells within range in this extract.

3.4 Climate (safran/, via download_era5_sample.py / download_era5_full.py)

Despite the safran naming throughout this codebase (a holdover from an earlier plan to use Météo-France's SAFRAN reanalysis), the actual data is ERA5 from Copernicus's Climate Data Store, pulled via cdsapi. ERA5 splits instantaneous variables (temperature, wind) from accumulated ones (precipitation, evaporation, radiation, snowfall, runoff) at the API level — download_era5_full.py downloads each set separately per year and merges them, because the CDS API rejects mixed requests. The full pull spans 1960–2026 and is genuinely slow.

SAFRANLoader interpolates the ERA5 grid to every station in one vectorized xarray call per file, not one .sel() + .to_dataframe() call per station — the per-station loop version does real per-call work (an index lookup, then a full DataFrame conversion) that's tolerable at 27 stations (1,800 calls across ~67 year-files) but was confirmed to actually hang at the reach graph's ~4,500 nodes (193,000 calls). Vectorized indexing with DataArray indexers sharing a station dimension does every station in one call per file instead.

Note: this is real historical reanalysis — a model-based reconstruction of what already happened, not a forecast. §3.13 covers the separate, genuine forecast archive used for real advance-warning signal in the model (§4.5).

3.5 IDPR (idpr.csv)

BRGM's Indice de Développement et de Persistance des Réseaux — an infiltration-vs-runoff tendency index, and the closest thing this project has to a real soil/drainage covariate. It's an integrated hydrological behavior indicator (infiltration tendency), not raw soil texture data, but arguably more directly useful for a streamflow model than a texture map would be on its own — paired with BD Charm-50 geology (§3.9) for the broader hydrological context soil data would otherwise provide. The file used here is already one row per station (station_id matching station_code exactly, verified 1:1 against all 27 stations), so node_features.py does a direct ID join when possible rather than nearest-neighbor search, falling back to spatial nearest-neighbor for any station code that isn't an exact match (every non-gauge reach-graph node, and — a real, minor precision trade-off worth knowing — every gauge too, once the table also contains non-gauge codes, since the exact-match path requires the entire table to match IDPR's station list).

3.6 Catchment area — two independent sources

Hub'Eau (catchment_area.csv, via scripts/download_catchment.py): published on the site referentiel, not the station referentiel — surface_bv on hydrometrie/referentiel/sites, in km². Since one site can have several stations, the download script does two passes: station → code_site, then code_site → surface_bv. 16 of 27 stations have a value. This number is cumulative — the total catchment area draining to that point, all the way to the source.

BD TOPO, graph-wide (cumulative_catchment_area_km2, via scripts/compute_cumulative_catchment.py): sums BD TOPO's incremental catchment polygons upstream of any node, via the real graph topology — distinct polygons counted once even when many nodes/edges share the same coarse polygon (verified with a hand-computed test case specifically checking this). Covers ~98% of nodes graph-wide, not just the 27 gauges — the actual fix for the "confluences and virtual nodes have no catchment area at all" gap.

Cross-checked against Hub'Eau's real values on real gauges — and there's a real, identified bias, not a clean match. Ratio (BD-TOPO-summed ÷ Hub'Eau) runs from about 0.75 to 1.25 for smaller catchments (< ~800 km², plausibly normal polygon-boundary/digitization precision) but drops to 0.75–0.89 for the largest catchments (> ~3,500 km²) — a clean, monotonic pattern, not noise. Most likely cause: bounding-box truncation — the original BD TOPO pull bbox had only a 9.6 km margin on its southern edge (the tightest of all four directions, and south is exactly where the Eure's longest upstream tributaries run, toward Chartres/Dreux), not a safe margin for real watershed extent. The bbox in download_bdtopo_hydro.py was widened afterward (from (0.3, 48.3, 1.7, 49.5) to (-0.1, 47.7, 2.1, 49.9), ~2.9x the area) — the full download_bdtopo_hydro.py → build_reach_graphs.py → enrich_reach_graph.py → compute_cumulative_catchment.py chain needs re-running against the wider box to actually resolve this, which had not yet happened as of the last verified run in this project.

In practice, the model (§4) currently uses a per-node fallback between the two sources, since cumulative_catchment_area_km2 isn't present at all on the real enriched node tables used for training as of the last verified run — catchment_area_km2 (Hub'Eau, real but gauge-only) is used wherever the cumulative BD TOPO value is missing. This recovers real catchment area for real gauges but not for confluences or braid points without a Hub'Eau delineation; those still have no real catchment area at all until the BD TOPO re-pull above is actually run.

3.7 BD TOPO hydrography (bdtopo_hydro/, via scripts/download_bdtopo_hydro.py)

IGN's BD TOPO / BD TOPAGE hydrographic network, pulled from the Geoplateforme WFS (https://data.geopf.fr/wfs) rather than downloaded as a national bulk file — the download script queries a bounding box around the two basins instead (see §3.6 for why that box was widened). Three layers, all scoped to that bbox:

• troncon_hydrographique.geojson — river centerline reaches, now the primary source for graph topology too (§2), via lien_vers_noeud_ hydrographique_ini/fin and sens_de_l_ecoulement. Real per-vertex altitude data doubles as a fine-grained elevation profile, denser than anything derivable from the 27 gauge points alone.
• surface_hydrographique.geojson — hydrographic surfaces, including a Nature attribute that's supposed to flag karst-influenced reaches (IGN documents this attribute as provisional and incomplete), and the real source for compute_edge_width.py's channel-width derivation (§2.2): polygon area ÷ tronçon length, joined by cleabs.
• bassin_versant_topographique.geojson — catchment polygons, incremental (see §3.6).

WFS axis order: when a BBOX parameter's CRS is given via the URN form, the OGC spec requires latitude, longitude axis order — the opposite of the lon,lat order most GIS tools use by default. Getting this backwards doesn't raise an error; it silently matches zero real features. download_bdtopo_hydro.py and scripts/download_bdcavites.py both try lon,lat first and automatically retry with the axes swapped if that comes back empty.

Real branching topology fixed a naive assumption. Filtering 30,045 tronçons down to a single named river and building a graph from their endpoints does not give one connected line — for "Risle" alone, 1,195 name-matched tronçons split into 132 disconnected components. Broadening the name filter to include known tributaries (§2.1) initially made this worse (487/214 components), traced to short/generic tributary names ("Bec", "Avre") matching unrelated streams elsewhere within the ~100×130 km bbox — fixed by requiring every name-matched tronçon to also fall within a real distance of a known gauge (load_troncons_for_basin's anchor_radius_km), and by selecting the connected component actually containing the most real gauges rather than the component with the most raw tronçons (best_component_for_stations) — proven to matter, not just theoretically: a synthetic adversarial test showed the naive "biggest component" approach picking a larger but entirely unrelated decoy network over the real one.

The bétoire finding: two stations in the roster are explicitly named "[amont bétoire]" and "[aval bétoire]" in Hub'Eau's own site names — bétoire being the Normandy dialect term for a karst swallow-hole. Three edges spanning that stretch on La Risle (H605641101 → H605022010 → H605641401 → H605641201) are flagged verified_continuous=False. BD TOPO's own karst attribute doesn't currently confirm it (see the provisional- attribute note above) — scripts/download_bdcavites.py (§3.8) exists specifically to get an independent, purpose-built second check on this, rather than relying only on naming inference.

3.8 BDCavités (bdcavites/, via scripts/download_bdcavites.py)

BRGM's national underground cavity inventory (sinkholes, quarries, natural cavities), via Géorisques' WFS (georisques.gouv.fr/services, typeName CAVITE_LOCALISEE, confirmed live and GeoJSON-capable directly against the real service). Built specifically as an independent check on the bétoire finding (§3.7) — a purpose-built cavity dataset, not inference from station naming or a provisional BD TOPO attribute. One real caveat: departments 75/78/91/92/93/94/95 (Paris region, unrelated to this project) are excluded from BDCavités entirely, and the Eure department's own inventory was among the later batches of the national 2001–2013 completion program — worth checking coverage density before treating a sparse result as a negative finding rather than incomplete data.

3.9 Geology (bdcharm50/, via scripts/download_bdcharm.py)

BRGM's BD Charm-50, harmonized 1:50,000 geological maps — free, open (Licence Ouverte), no authentication, direct per-department ZIP download from InfoTerre (a genuinely different access pattern than the WFS sources elsewhere in this project: fixed URL per department, no bbox query, no axis- order ambiguity). Departments 27 (Eure), 28 (Eure-et-Loir), 61 (Orne) — verified directly against the real, complete station roster (§3.1), not guessed. A separate, CIGAL-membership-gated distribution of similar data exists for at least one other French region; this project only uses the free InfoTerre path.

3.10 Landcover and NDVI (scripts/fetch_landcover.py, scripts/fetch_worldcover_ndvi.py)

ESA WorldCover, sampled at real gauge points from the public AWS S3 Cloud- Optimized GeoTIFFs.

Landcover classification uses the product's 3°×3° tile grid; every real station coordinate falls inside exactly one tile (N48E000), verified directly against all 27 real coordinates. NDVI uses the annual composites' 1°×1° tile grid instead — genuinely different from the classification grid, looked up per-station via VITO's own authoritative tile-index grid file (esa_worldcover_grid_composites.fgb) rather than a hand-guessed S3 key pattern. Both need AWS_NO_SIGN_REQUEST=YES for s3://-scheme tile URLs specifically — a plain HTTPS URL to the same public bucket needs no signing at all.

Both currently cover only the 27 real gauges (exact station_code match), same limitation as Hub'Eau's catchment_area_km2 before the cumulative-BD- TOPO fix (§3.6) — extending either script to the full reach graph is unstarted work, not a design decision.

3.11 Centerline generation

Only relevant to the older single-chain pipeline (build_surface_edges, still available for direct comparison/debugging) — the reach graph (§2) derives its topology directly from BD TOPO's own node linkage and doesn't use these centerline files at all.

centerlines/eure_centerline.csv and centerlines/risle_centerline.csv — the geometry build_surface_edges orders stations against — are generated by scripts/analyze_bdtopo_hydro.py --export-centerline. It filters troncon_hydrographique.geojson (§3.7) down to the named river, builds a graph from the tronçon endpoints, and walks the longest path through it via double-BFS shortest-path to get one continuous, correctly-ordered sequence of real coordinates. This is the accurate method: real BD TOPO vector geometry snaps stations to within 0.05 km on average.

scripts/extract_river_centerline.py is a separate, standalone technique for deriving a centerline directly from a traced map image, for a river or region without BD TOPO coverage: color-threshold the image to isolate a traced route, skeletonize it, walk the end-to-end path the same double-BFS way, then georeference by fitting a least-squares affine transform from a handful of manually-read reference-point pixel positions to their known coordinates. This produces a reliable shape, but the absolute position is only as good as the georeferencing step — residuals at the reference points run to a few kilometers with a handful of manually-read points, giving roughly a 0.88 km average snap distance rather than 0.05 km. It's the fallback when a real vector source isn't available, not the method used for the current centerlines/ files.

3.12 Real historical precipitation (meteofrance_precipitation/)

Real, directly-measured rain-gauge observations from Météo-France's public "Données climatologiques de base — horaires" dataset — deliberately not ERA5 reanalysis (§3.4) for this specific variable: precipitation varies sharply over short distances in ways a ~9 km reanalysis grid cell can miss, so a real nearby gauge is more locally authoritative where one exists. Official source, Open Licence 2.0, no authentication — fetch_meteofrance_precipitation.py fetches only the departments this project needs (27, 28, 61) and only the decades overlapping the real study period, aggregating real hourly RR1 readings to daily sums per station.

Each graph node reads from its nearest real Météo-France gauge by haversine distance (map_stations_to_precipitation_gauges) — several nodes can share one real gauge if they're close together, which is expected, not a bug, the same way several climate stations already share one ERA5 grid cell elsewhere in this project.

Also used to estimate discharge for stations with no real gauge data at all (§3.1): fit_specific_discharge_coefficient fits a simple, standard hydrological scaling relationship (Q ≈ k · precipitation · catchment_area, a real regionalization technique) using only real (discharge, precipitation, catchment area) triples from gauged stations during the training period, then applies it to every node with real precipitation and catchment area, including the 14 with no real discharge history. This is wired in as a model input feature, never a target — it's a real, physically-grounded estimate, not a measurement, and reaches the model exactly the way every other imperfect-but-real input does.

3.13 Real ECMWF forecast archive (previous_runs_forecast.csv)

Real historical forecasts — not the live current forecast, and a genuinely different thing from ERA5 reanalysis (§3.4): reanalysis is a best estimate of what already happened, with zero forecast skill baked in; this archive is what was actually predicted in advance, which is what a model needs to learn genuine multi-day-ahead skill from. Fetched via fetch_previous_runs.py (Open-Meteo's Previous Runs API).

Real, confirmed schema: one row per real observation date, with several precipitation_previous_dayN columns — each holding what was forecasted N days before that date, for that date. For an anchor date t and a horizon h, the real forecast value lives in the row for date = t + h, column precipitation_previous_day{h} — not the row for t itself, since that forecast was issued on (t+h) - h = t, the anchor. station_code in this file already matches this project's real graph node codes directly — no nearest-neighbor mapping needed here, unlike §3.12's precipitation gauges.

Real, significant coverage limits, not implementation gaps: only lead times [1, 2, 3, 5, 7] days have any real data at all — none of this project's other horizons (4, 10, 30, 60, 90, 180) have real forecast data, ever. Coverage also only starts 2024-03-01 — the large majority of this project's 2013+ training window has no real forecast data regardless of lead time. Both gaps are handled with an explicit real/missing flag per (example, node, horizon), the same value-plus-missingness pattern used throughout this project, not a silent fill.

A real, confirmed gotcha worth flagging directly: with this project's default --train-end/--val-end split, the entire real forecast archive falls after the validation cutoff — meaning the model can see 0% real forecast coverage during training and validation, then 100% at evaluation time. That mismatch produced a genuine, measured collapse in evaluation performance at exactly the horizons where real forecast coverage is densest, since the model's forecast-conditioning pathway was never trained on any real (non-missing) input. The fix is a --train-end/--val-end shift so all three splits genuinely overlap the real forecast-covered period — see §4.5 and §7.

---

4. The model

4.1 Architecture (spatiotemporal_gnn.py)

SpatiotemporalGNN predicts discharge and water level at any node — gauged or ungauged — for a requested lead time, trained end-to-end against the real graph and physics substrate above.

• Input encoding (TypeAwareInputEncoder): gauged nodes get a dedicated encoder path over real static features plus dynamic history; non-gauged nodes get a learned placeholder, rather than forcing every node through the same encoder as if the ~14 uninstrumented stations (§3.1) had real time series to offer.
• Spatial (RiverMessagePassing, two stacked layers): directed message passing over the real reach graph, respecting real flow direction, with sum aggregation, not mean — deliberate, not a default left unconsidered: a confluence's downstream discharge is the sum of its upstream branches (§2.5), so the aggregation itself reflects that physical reality rather than averaging it away.
• Temporal (a per-node GRU): runs over the real 30-day lookback window, fed by the spatial layers' output at each timestep — spatial happens first, then the resulting sequence of embeddings goes through the GRU, so the recurrent step learns genuine spatiotemporal patterns rather than each node's isolated history.
• Decoder: combines the GRU's final hidden state with a basin embedding, a horizon embedding (letting one trained model answer at ten different lead times without retraining separately for each), and — when available — a real forecast-precipitation embedding (§4.5), then outputs quantile predictions for discharge and water level.

4.2 Quantile regression, not a point estimate (pinball_loss.py)

The model predicts several quantiles ([0.5, 0.9, 0.95, 0.99] by default) per target, not a single value — a real, evidence-driven architecture decision, not a default choice. Real discharge is heavily right-skewed (confirmed against real data: the 99th percentile runs 4–12x the median), and a model trained on plain MSE structurally minimizes average error, which means predicting close to the mean is a genuinely good strategy for MSE and a genuinely bad one for flood risk — confirmed directly, not just argued: an earlier MSE-trained version of this model caught 0 of 203 real historical flood events at every tested lead time (§4.7), despite showing real average- case skill.

Pinball loss is asymmetric by construction: for a high quantile (e.g. 0.95), under-predicting is penalized far more than over-predicting by the same amount, which is the actual mechanism that pushes the model to cover the real tail rather than hide near the mean.

Quantile crossing is architecturally impossible, not just discouraged. A well-known failure mode of naive quantile regression is the predicted 90th percentile coming out lower than the predicted median — physically nonsensical. The lowest quantile is predicted directly; every subsequent one is the previous quantile plus a softplus-transformed (always non-negative) offset, so quantile[k] ≥ quantile[k-1] holds for every real input, including deliberately adversarial ones.

The naive baseline this model is compared against is also quantile-aware — each quantile's naive prediction is that quantile's own real empirical value from the training-period distribution (the direct quantile-regression equivalent of climatology as a baseline in the earlier climate work), not zero or the mean for every quantile alike.

4.3 Physics losses, wired into training

All five terms from §2.5 are combined per training example, horizon-weighted rather than with one fixed weight applied uniformly — routing and water balance specifically only apply to the short contiguous-horizon block (the first 5 daily steps), since both need a real, physically comparable period, not a mismatched one.

Two real numerical-stability bugs were found and fixed during training, worth naming since both generalize beyond this specific project:

• Manning's equation (1/n, hydraulic_radius^(2/3)) has genuinely steep local derivatives near its natural boundaries — an untrained model's raw output is close to unconstrained noise, and a single bad step could push n toward zero, producing a real, not hypothetical, inf loss (confirmed: two real training runs failed exactly this way before the fix). n, bank height, hydraulic radius, and depth are all explicitly clamped to physically real ranges before use.
• A gradient-masking gotcha in water_balance_loss (§2.5): a value correctly excluded from the loss can still corrupt its gradient if NaN entered the computation before masking, since 0 × NaN = NaN. Fixed by replacing NaN with a safe finite placeholder before any arithmetic, not by relying on the final mask alone.

A real, current limitation worth being direct about: Manning's learned roughness n has landed on implausible values across real training runs (0.38–0.92, against a real natural-stream range of roughly 0.025–0.1), and varies meaningfully run to run rather than converging to a stable estimate. The likely mechanism is clamp saturation — the value used inside the loss computation is clamped to [0.01, 0.3] for numerical safety, and once the underlying parameter drifts past that ceiling, the clamp's own gradient is zero there, removing further corrective signal. This means Manning's equation has likely been running at its artificial ceiling for meaningful stretches of training rather than converging to a physically informed value — worth investigating (a wider clamp range, or a gradient that doesn't vanish at the boundary) before trusting this specific learned parameter.

4.4 Real precipitation and a specific-discharge estimate

Covered fully in §3.12. In short: real Météo-France precipitation feeds the temporal encoder directly as a dynamic input channel, and — for the 14 stations with no real discharge history — a real, physically-grounded specific-discharge estimate (fit only on real gauged data, applied uniformly to every node) gives the model a genuine signal to work with instead of nothing, without ever fabricating a value the supervised loss is held accountable to match.

4.5 Real forecast integration

Covered fully in §3.13. The model's decoder accepts real forecasted precipitation per (node, horizon), defaulting to fully-missing when no real forecast dataset is wired in — backward compatible with every training run before this was added. Given the real, confirmed train/test coverage mismatch documented in §3.13, using this correctly requires shifting --train-end/--val-end so training and validation genuinely overlap the real 2024-03-01+ forecast-covered period (see §7 for the actual command) — using the project's older default split dates with forecast data wired in will reproduce the coverage-mismatch failure, not avoid it.

4.6 Subgraph reduction (subgraph_selection.py)

Training runs against a real, reduced 100-node subgraph (50 nodes per basin: real gauges plus real confluences/braid points only — virtual infill nodes are deliberately excluded, not because they're uninteresting, but because training on the full ~4,500-node graph per basin was not tractable within this project's real compute constraints. Chains through dropped nodes are collapsed into single edges (summing distance/elevation, distance- weighting average width) rather than simply deleted, so real connectivity between kept nodes is preserved exactly. A full-graph retrain, once the 100-node architecture is validated, remains real, unstarted future work.

4.7 Real evaluation

identify_flood_events.py establishes real ground truth independent of any model: real historical high-flow events, defined per station by that station's own real historical percentile (95th by default — a station-relative, not a shared absolute, threshold, since these rivers have very different real scales), with consecutive above-threshold days grouped into discrete events rather than counted individually.

evaluate_flood_detection.py checks the trained model against this real ground truth on the held-out real test split — flagging a real event when the model's predicted upper quantile (0.95 by default, matching the same percentile identify_flood_events.py uses) exceeds that station's real threshold, not the median, since flood risk is a tail question. Precision and recall are reported per lead time, since catching a flood a day out and catching it a week out are different questions that would be hidden by one combined number.

Real, current results (subject to change as the model continues to be refined — treat these as a snapshot, not a final claim):

• The naive-baseline comparison (§4.2) shows real, meaningful skill: the quantile model beats the real empirical-quantile baseline by roughly 30–50% on both discharge and water level across real training runs.
• Real flood-detection recall at short lead times (1–5 days) has reached roughly 65–70% in real evaluation runs — a genuine, validated improvement over the earlier MSE-trained model's 0/203.
• Precision remains low in absolute terms (under 10% at every tested horizon in the most recent full evaluation) — real, meaningful skill above the base rate at short horizons, but still many false alarms; not yet something to treat as deployment-ready without further tuning (e.g. a higher --flag-quantile, trading recall for precision).
• Long-horizon "skill" needs a base-rate check, not just a recall number: at the longest tested horizon (180 days), precision has landed below the real base rate in at least one evaluation run — a sign the model is likely hedging with a wide quantile spread rather than genuinely predicting that far out, not a sign of real long-range foresight. evaluate_flood_ detection.py's output should always be read against the real base rate ((tp+fn)/total), not the recall number alone.

---

5. Applications

src/app.py (Streamlit) has two views, selected by a radio at the top:

Explore — the original click-to-read UI over real course geometry: pick a river, click (or slide) along its course, see interpolated elevation, estimated groundwater level, and — for whichever real gauge is nearest that point — water level, discharge, and rating-curve plots pulled directly from HydrometricLoader's own plotting methods rather than reimplemented.

Click support uses Streamlit's native chart-selection (st.plotly_chart(..., on_select="rerun")), not a third-party click-handling package. The click handler and the position slider share a single source of truth by design: Streamlit only honors a slider's value= argument the first time that widget is created, and on every later rerun returns whatever's stored under that widget's own session-state key — so the click handler writes directly into the slider's own key before it's instantiated, rather than a separate key. It also de-duplicates incoming click events, since Streamlit's chart-selection state persists across reruns caused by other widgets and would otherwise re-fire on every unrelated interaction.

Network validation — renders the full reach graph (§2): every edge as one Plotly line trace regardless of edge count (a trace-per-edge approach doesn't hold up at ~5,000+ edges; verified fast at real scale — 0.29s to build a figure for ~2,900 edges), real confluences as diamond markers, real gauges as elevation-colored circles. Metrics card reports node/edge/confluence/gauge counts and, when available, IDPR and cumulative-catchment coverage. Virtual infill nodes are deliberately not drawn individually — at ~2,400 per basin, markers for each would bury the actual validation signal (do confluences sit where a tributary visibly joins the line? do gauges sit on the network, not offset from it?) rather than help it. Reads directly from reach_graph/ {basin}_nodes_enriched.csv, keyed on file modification time so a re-run of build_reach_graphs.py/enrich_reach_graph.py is picked up automatically — st.cache_data otherwise keys purely on function arguments, not file contents, and this was confirmed to actually cause stale numbers once during development, not just a theoretical risk.

---

6. Testing (src/test_build_graph.py)

Not a unit test suite in the pytest sense — a script with two independent sections, both run from main().

run_checks — the original single-chain pipeline: runs node_features → build_surface_edges → build_pyg_graph(s) against real data and checks the result is sane — no NaN/Inf in the feature tensor, no accidental cross-basin edges, targets genuinely excluded from the model input, edge indices within bounds, bidirectional edge count exactly double the directed count, per-basin node counts summing to the combined total, standardized features actually landing near zero mean / unit variance, the known_losing_reaches flag actually taking effect, and mean/max snap_distance_km per basin against whatever centerline is currently in centerlines/.

run_reach_graph_checks — the reach graph pipeline, gracefully skipped (not a failure) if reach_graph/ doesn't exist yet. Mostly regression tests for three bugs found and fixed during development, kept here specifically so they can't silently reintroduce themselves:

• structural columns (is_gauged/is_confluence/etc.) never leak into feature_names, but remain accessible as their own Data attributes
• target values never attach to a non-gauge node, and target coverage never exceeds the real gauge count
• edge_attr stays exactly the real numeric edge columns despite extra edge metadata (toponym, cleabs) sitting on the real edges table
• physics_losses.py's build_confluence_index/build_braid_index produce counts matching is_confluence/is_rejoin_point sums, with every index within node bounds and every confluence having ≥ 2 upstream branches
• IDPR and cumulative_catchment_area_km2 presence/coverage are reported explicitly (the latter compared against the Hub'Eau-only baseline it's meant to exceed)

Exits 0 on a clean pass across both sections, 1 otherwise — usable as a pre-commit or CI gate if that's ever set up.

identify_flood_events.py and physics_losses.py's own functions (§2.5, §4.7) double as an informal, model-independent validation layer too — real historical data checked directly against real physics and real event definitions, without needing a trained model in the loop at all.

---

7. Known limitations and open questions

• cumulative_catchment_area_km2 underestimates the largest catchments by up to ~25%, traced to the BD TOPO pull's original bounding box having an insufficient southern margin. The bbox has been widened in download_bdtopo_hydro.py; the full re-pull-and-rebuild chain needs re-running for this to actually resolve — and, in the meantime, this column is entirely absent from the real enriched node tables the model currently trains against (§3.6), which is why the model falls back to the gauge-only Hub'Eau value.
• Climate's genuine time series (build_climate_timeseries) is untested against real data — discharge and groundwater's equivalents are; verify climate's real output before relying on it.
• Landcover and NDVI only cover the 27 real gauges, not the full reach graph — same scope catchment_area_km2 had before its cumulative-BD-TOPO extension.
• The model trains and evaluates against a 100-node subgraph, not the full ~4,500-node-per-basin graph — a real, deliberate scope limit for tractable compute, not a claim the full graph wouldn't matter. A full-graph retrain is real, unstarted future work.
• Manning's learned roughness n is currently not trustworthy as a physical finding — real evidence (§4.3) points to clamp saturation, not genuine convergence. Needs either a wider clamp range or a non-vanishing-at-the-boundary gradient before the learned value means anything physically.
• Real forecast precipitation coverage is genuinely partial — five lead times, roughly 2.5 years — and using it correctly requires a real, documented train/val/test split shift (§3.13, §4.5); the project's original default split dates predate the entire forecast archive.
• Flood-detection precision is low in absolute terms even where recall is real and meaningful (§4.7) — worth tuning the flag quantile and continuing to validate before treating this as deployment-ready.
• Real evapotranspiration is not wired into water_balance_loss — passed as an explicit 0.0 placeholder; real ERA5 evaporation data exists but lives on a different grid than this graph's nodes, and building that mapping is real, unstarted work.
• The karst losing-reach flag rests on naming evidence and a BDCavités cross-check (§3.7–3.8), not a fully confirmed BD TOPO classification.
• The groundwater-well BDLISA aquifer-unit field is unused. First place to look if a subsurface connectivity edge is ever justified with real evidence rather than proximity.

---

8. Running things

Data acquisition (from repo root, in roughly dependency order):

python -m scripts.data_acquisition.download_hubeau
python -m scripts.data_acquisition.download_elevation
python -m scripts.data_acquisition.download_era5_full        # slow; download_era5_sample.py first if just testing
python -m scripts.data_acquisition.extract_era5
python -m scripts.data_acquisition.download_catchment
python -m scripts.data_acquisition.download_bdtopo_hydro --check   # verify typeNames before the real pull
python -m scripts.data_acquisition.download_bdtopo_hydro
python -m scripts.data_acquisition.download_bdcavites --check
python -m scripts.data_acquisition.download_bdcavites
python -m scripts.data_acquisition.download_bdcharm
python -m scripts.data_acquisition.fetch_meteofrance_precipitation --check --department 27
python -m scripts.data_acquisition.fetch_meteofrance_precipitation --data-root datasets
python -m scripts.data_acquisition.fetch_previous_runs        # real ECMWF forecast archive

Build and validate the reach graph:

python -m scripts.graph_building.build_reach_graphs --data-root datasets
python -m scripts.graph_building.enrich_reach_graph --data-root datasets   # --skip-climate if that step hangs
python -m scripts.graph_building.compute_cumulative_catchment --data-root datasets
python -m scripts.graph_building.compute_edge_width --data-root datasets
python -m scripts.graph_building.diagnose_confluences --data-root datasets --basin eure
python -m scripts.graph_building.diagnose_confluences --data-root datasets --basin risle

Build genuine [n_nodes, T] dynamic tensors and verify the physics-loss wiring:

python -m scripts.graph_building.build_dynamic_tensors --data-root datasets --basin risle
python -m scripts.graph_building.build_dynamic_tensors --data-root datasets --basin eure

Landcover / NDVI, real gauges only (needs rasterio, and geopandas for NDVI's tile lookup):

python -m scripts.data_acquisition.fetch_landcover --check
python -m scripts.data_acquisition.fetch_landcover
python -m scripts.data_acquisition.fetch_worldcover_ndvi --check
python -m scripts.data_acquisition.fetch_worldcover_ndvi

Real historical flood events (ground truth, needed before evaluation):

python -m scripts.evaluation.identify_flood_events --data-root datasets

Train the model — with real forecast data, use the shifted split dates (§3.13, §4.5) so training and validation actually overlap real forecast coverage:

python -m scripts.training.train_spatiotemporal_gnn --data-root datasets \
    --train-end 2025-06-30 --val-end 2025-12-31

Evaluate the trained model against real historical flood events:

python -m scripts.evaluation.evaluate_flood_detection --data-root datasets

Validate the graph itself against whatever's actually in datasets/:

python -m src.test_build_graph --data-root datasets

Run the explorer:

streamlit run src/app.py -- --data-root datasets