Spaces:
Running on Zero
Running on Zero
| 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<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: | |
| ```mermaid | |
| 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): | |
| ```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 | |
| ``` |