--- 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](#4-the-model) for the architecture and real, current evaluation results. ```mermaid flowchart LR hubeau["Hub'Eau
discharge Β· water level
catchment area"] ades["ADES
groundwater levels"] era5["Copernicus ERA5
climate reanalysis"] otd["Open Topo Data
station elevation"] brgm["BRGM
IDPR Β· BD Charm-50 geology"] bdtopo["IGN BD TOPO
real reach topology + catchment polygons"] bdcav["GΓ©orisques
BDCavitΓ©s (sinkholes)"] wc["ESA WorldCover
landcover Β· NDVI"] mf["MΓ©tΓ©o-France
real historical precipitation"] ecmwf["ECMWF (Open-Meteo)
real forecast archive"] bdtopo --> brg["build_reach_graph.py
real confluences, splits/rejoins,
gauge snapping"] brg --> brgs["build_reach_graphs.py
~4,500 nodes/basin"] hubeau --> nf ades --> nf era5 --> nf otd --> nf brgm --> nf bdcav --> nf wc --> nf brgs --> nf["node_features.py /
enrich_reach_graph.py
date-filtered 2013-2026"] bdtopo --> cc["compute_cumulative_catchment.py
graph-wide catchment area"] cc --> nf nf --> pyg["build_pyg_graph
x_static / x_dynamic split"] pyg --> phys["physics_losses.py
confluence Β· split-rejoin Β·
routing Β· manning Β· water balance"] mf --> train ecmwf --> train phys --> train["train_spatiotemporal_gnn.py
quantile regression + physics losses,
100-node subgraph"] train --> model["model.pt
trained checkpoint"] model --> eval["evaluate_flood_detection.py
vs. real historical flood events"] pyg --> app["src/app.py
Streamlit explorer +
network validation view"] pyg --> testsuite["test_build_graph.py
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: ```mermaid flowchart TD raw["Enriched node table
(node_features.py)"] raw --> struct{"structural /
graph-role column?"} struct -->|"is_gauged, is_confluence,
is_split_point, is_rejoin_point,
braid_id, snap_distance_km"| structout["data.is_gauged, data.is_confluence, ...
own Data attribute β€” never in x"] raw --> tgt{"target_* column?"} tgt -->|"target_discharge_m3s_*
target_waterlevel_mm_*"| y["data.y
never in x β€” label leakage otherwise"] raw --> feat{"real model input"} feat -->|"static: elevation_m, idpr_*,
catchment_area_km2, landcover_*,
geology_*, cavites distance/count"| xstatic["data.x_static"] feat -->|"dynamic: climate_*,
avg_groundwater_*, ndvi_*
(period-aggregate, not a real series yet)"| xdynamic["data.x_dynamic"] xstatic --> x["data.x β€” full combined tensor,
z-scored"] xdynamic --> x edges["Edge table
(build_reach_graph_tables)"] --> eattr{"numeric edge
attribute?"} eattr -->|"distance_km,
elevation_drop_m,
verified_continuous,
width_m"| edgeattr["data.edge_attr
[n_edges, 4]"] eattr -->|"toponym, cleabs
(diagnostic metadata)"| meta["not used by build_pyg_graph β€”
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. ![La Risle reach graph β€” network validation view, real confluences as diamonds, gauges elevation-colored](docs/images/risle_reach_graph.png) ![La Eure reach graph, same view](docs/images/eure_reach_graph.png) ![Explore view: clicking along the river resolves to the nearest gauge and shows its real discharge/water-level/rating-curve plots](docs/images/explore_view.png) --- ## 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): ```bash 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: ```bash 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: ```bash 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): ```bash 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): ```bash 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: ```bash 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: ```bash python -m scripts.evaluation.evaluate_flood_detection --data-root datasets ``` Validate the graph itself against whatever's actually in `datasets/`: ```bash python -m src.test_build_graph --data-root datasets ``` Run the explorer: ```bash streamlit run src/app.py -- --data-root datasets ```