README.md

---
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.

![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
```