Add routing module: steady-state runoff accumulation + peak discharge (1/2)

README.md

@@ -75,6 +75,7 @@ print(f"Total damaged built-up: {damage.values.sum():,.0f} m²")
 | `floodpath.soil` | ISRIC SoilGrids 2.0 (250 m, COG) | Sand/silt/clay topsoil composition + USDA texture-triangle classification + NEH 630 Ch7 hydrologic soil group (A/B/C/D) |
 | `floodpath.precip` | Synthetic uniform (real fetchers later: ERA5 / IMERG / CHIRPS) | Precipitation depth raster (mm) — pluggable input to the runoff equation |
 | `floodpath.runoff` | NEH 630 Ch9 + Ch10 + landuse + HSG + precip | SCS Curve Number raster + SCS-CN runoff equation `Q = (P-0.2S)²/(P+0.8S)` |
+| `floodpath.routing` | runoff + flow direction (pyflwdir) | Steady-state hydrologic routing: accumulated runoff volume + peak discharge |
 | `floodpath.damage` | JRC Huizinga 2017 + DEM/HAND/GHSL | Per-cell flood depth and damage in m² of built-up surface |
 
 ## Depth-damage curves

floodpath/routing/__init__.py

@@ -0,0 +1,32 @@
+"""Routing module: hydrologic + (later) hydraulic routing of runoff.
+
+Currently provides steady-state hydrologic routing — accumulation of
+SCS-CN runoff along the DEM-derived flow direction grid, plus a
+volume-to-peak-discharge conversion given an event duration.
+
+Manning normal-depth (the hydraulic closure that converts discharge to
+water level at stream cells) is the natural next addition; it lives in
+this module's `manning.py` once shipped.
+"""
+
+from .accumulate import accumulate_runoff
+from .constants import (
+    DEFAULT_EVENT_DURATION_S,
+    EARTH_RADIUS_M,
+)
+from .discharge import peak_discharge
+from .models import (
+    AccumulatedRunoffGrid,
+    DischargeGrid,
+)
+from .utils import cell_areas_m2
+
+__all__ = [
+    "AccumulatedRunoffGrid",
+    "DEFAULT_EVENT_DURATION_S",
+    "DischargeGrid",
+    "EARTH_RADIUS_M",
+    "accumulate_runoff",
+    "cell_areas_m2",
+    "peak_discharge",
+]

floodpath/routing/accumulate.py

@@ -0,0 +1,101 @@
+"""Accumulate runoff (mm depth) along a flow direction grid into volume (m^3)."""
+
+import numpy as np
+import rasterio.warp
+from rasterio.enums import Resampling
+
+from floodpath.dem.models import BoundingBox
+from floodpath.hydrology.models import FlowGrid
+from floodpath.runoff.models import RunoffGrid
+
+from .constants import ROUTING_METHOD
+from .models import AccumulatedRunoffGrid
+from .utils import cell_areas_m2
+
+
+def _runoff_at_flow_grid(
+    runoff: RunoffGrid,
+    flow_grid: FlowGrid,
+) -> np.ndarray:
+    """Reproject runoff depth (mm) onto the flow grid via area-weighted average.
+
+    Runoff is a depth (mm) — averaging fine cells into a coarse cell is
+    correct (the coarse cell's depth is the area-weighted mean of its
+    fine cells'). Volume aggregation comes later, after multiplying by
+    the coarse-grid cell area.
+    """
+    out = np.zeros(flow_grid.shape, dtype=np.float32)
+    rasterio.warp.reproject(
+        source=runoff.values,
+        destination=out,
+        src_transform=runoff.transform,
+        src_crs=runoff.crs,
+        dst_transform=flow_grid.transform,
+        dst_crs=flow_grid.crs,
+        resampling=Resampling.average,
+    )
+    return out
+
+
+def accumulate_runoff(
+    runoff: RunoffGrid,
+    flow_grid: FlowGrid,
+) -> AccumulatedRunoffGrid:
+    """Accumulate runoff downstream along the flow direction grid.
+
+    Three steps:
+
+    1. **Reproject** runoff depth (mm) from its native grid (typically
+       WorldCover 10 m) to the flow grid (typically DEM 30 m) via
+       area-weighted average. Each coarse cell now carries the
+       area-weighted mean depth of its fine cells.
+    2. **Convert** depth (mm) at each cell to volume (m^3) by multiplying
+       by the cell's area in m^2 (computed from the transform via
+       `cell_areas_m2`, which handles WGS84 cos(lat) convergence).
+    3. **Accumulate** the per-cell volume downstream using
+       `pyflwdir.FlwdirRaster.accuflux`. Each cell's value becomes its
+       own volume plus the sum of all upstream cells'.
+
+    NaN cells in the input runoff grid (typically CN nodata at the HSG
+    reprojection edge) are zero-filled before accumulation — they
+    contribute no runoff downstream rather than poisoning the entire
+    accumulation.
+
+    Args:
+        runoff: Per-cell runoff depth (mm) from `apply_scs_cn`.
+        flow_grid: Flow direction + accumulation grid from `build_flow_grid`.
+
+    Returns:
+        AccumulatedRunoffGrid with per-cell upstream runoff volume (m^3),
+        sharing the flow grid's georef.
+    """
+    # 1. Reproject runoff depth (mm) onto the flow grid.
+    q_at_flow = _runoff_at_flow_grid(runoff, flow_grid)
+    # NaN -> 0 (unknown -> contributes nothing). Without this, accuflux
+    # would propagate NaN through the entire downstream network.
+    q_at_flow = np.where(np.isnan(q_at_flow), 0.0, q_at_flow)
+
+    # 2. Volume per cell: depth (mm) * area (m^2) / 1000 mm/m -> m^3.
+    areas = cell_areas_m2(
+        transform=flow_grid.transform,
+        shape=flow_grid.shape,
+        crs=flow_grid.crs,
+    )
+    volume_per_cell = (q_at_flow * areas / 1000.0).astype(np.float32)
+
+    # 3. Accumulate downstream via pyflwdir.
+    accumulated = flow_grid.flwdir_raster.accuflux(volume_per_cell)
+
+    return AccumulatedRunoffGrid(
+        values=accumulated.astype(np.float32),
+        transform=flow_grid.transform,
+        crs=flow_grid.crs,
+        bbox=BoundingBox(
+            min_lon=flow_grid.bbox.min_lon,
+            min_lat=flow_grid.bbox.min_lat,
+            max_lon=flow_grid.bbox.max_lon,
+            max_lat=flow_grid.bbox.max_lat,
+        ),
+        precip_source=runoff.precip_source,
+        method=ROUTING_METHOD,
+    )

floodpath/routing/constants.py

@@ -0,0 +1,23 @@
+"""Constants for the routing module."""
+
+# Approximate radius of the Earth (m). Used to convert WGS84 pixel sizes
+# in degrees to areas in square metres. Mean of equatorial and polar
+# radii is ~6,371,000 m; for the small-bbox use cases we care about the
+# residual error is well below one percent.
+EARTH_RADIUS_M: float = 6_371_000.0
+
+# Length of one degree of latitude at the surface, in metres. Constant
+# across the globe to first order (latitude lines are circles around the
+# polar axis only at the equator; everywhere else they shorten by cos(lat)).
+DEG_LAT_M: float = 111_319.5  # = 2 * pi * EARTH_RADIUS_M / 360 * (something close)
+
+# Default event duration for the steady-state peak-discharge conversion.
+# 6 hours is the standard NRCS / FEMA "design storm" duration for small
+# rural watersheds; users with their own event windows pass `duration_s`
+# explicitly.
+DEFAULT_EVENT_DURATION_S: float = 6 * 3600.0  # 21,600 s
+
+# Source-attribution string for the AccumulatedRunoffGrid / DischargeGrid
+# when produced by these helpers. Future routing schemes
+# (unit hydrograph, kinematic wave) will use distinct strings.
+ROUTING_METHOD: str = "steady-state-flow-accumulation"

floodpath/routing/discharge.py

@@ -0,0 +1,51 @@
+"""Convert accumulated runoff volume to peak discharge."""
+
+import numpy as np
+
+from .models import AccumulatedRunoffGrid, DischargeGrid
+
+
+def peak_discharge(
+    accumulated: AccumulatedRunoffGrid,
+    duration_s: float,
+) -> DischargeGrid:
+    """Convert accumulated upstream runoff volume to peak discharge.
+
+    Steady-state estimate: per-cell peak discharge equals the total
+    upstream runoff volume divided by the event duration.
+
+        Q_peak (m^3/s) = V_acc (m^3) / T (s)
+
+    This treats the entire upstream runoff as arriving over the same
+    event window — accurate for storms much shorter than the basin
+    time-of-concentration (small mountain basins, intense convective
+    storms), increasingly biased high for larger basins where peak
+    attenuation along the channel matters. A future kinematic-wave or
+    unit-hydrograph router will replace this for those cases.
+
+    Args:
+        accumulated: Output of `accumulate_runoff`.
+        duration_s: Event duration in seconds. Use 6*3600 for the
+            standard NRCS 6-hour design storm; pass `DEFAULT_EVENT_DURATION_S`
+            from constants for a typed default.
+
+    Returns:
+        DischargeGrid with per-cell peak discharge in m^3/s.
+
+    Raises:
+        ValueError: `duration_s` is not strictly positive.
+    """
+    if duration_s <= 0:
+        raise ValueError(f"duration_s must be > 0, got {duration_s}")
+
+    q_peak = (accumulated.values / float(duration_s)).astype(np.float32)
+
+    return DischargeGrid(
+        values=q_peak,
+        transform=accumulated.transform,
+        crs=accumulated.crs,
+        bbox=accumulated.bbox,
+        duration_s=float(duration_s),
+        precip_source=accumulated.precip_source,
+        method=accumulated.method,
+    )

floodpath/routing/models.py

@@ -0,0 +1,98 @@
+"""Dataclasses for the routing module."""
+
+from dataclasses import dataclass
+
+import numpy as np
+from rasterio.transform import Affine
+
+from floodpath.dem.models import BoundingBox
+
+
+@dataclass
+class AccumulatedRunoffGrid:
+    """Accumulated upstream runoff volume (m^3) at each cell.
+
+    `values[i, j]` is the total volume of water that has flowed into cell
+    (i, j) from itself and all upstream cells in the DEM-derived flow
+    direction grid, given the per-cell runoff depth from a RunoffGrid.
+    Stream cells thus have the largest values; ridge cells have the
+    smallest (just their own per-cell volume).
+
+    Output of `accumulate_runoff(runoff, flow_grid)`.
+    """
+
+    values: np.ndarray
+    transform: Affine
+    crs: str
+    bbox: BoundingBox
+    precip_source: str
+    method: str
+    units: str = "m^3 (accumulated runoff volume)"
+
+    @property
+    def shape(self) -> tuple[int, int]:
+        """Return the (rows, cols) shape of the underlying raster."""
+        return self.values.shape  # type: ignore[return-value]
+
+    def total_volume_m3(self) -> float:
+        """Bulk volume integrated over the patch — *not* the sum of
+        accumulated values (which double-counts each cell along its
+        downstream path). Use this when reporting water-balance totals.
+        """
+        return float(self.values.max())
+
+    def stats(self) -> dict[str, float]:
+        """Summary statistics over finite cells."""
+        valid = self.values[~np.isnan(self.values)]
+        if valid.size == 0:
+            return {"min": 0.0, "max": 0.0, "mean": 0.0, "median": 0.0}
+        return {
+            "min": float(valid.min()),
+            "max": float(valid.max()),
+            "mean": float(valid.mean()),
+            "median": float(np.median(valid)),
+        }
+
+
+@dataclass
+class DischargeGrid:
+    """Per-cell peak discharge Q_peak (m^3/s) under steady-state routing.
+
+    Computed as accumulated upstream runoff volume divided by event
+    duration: Q_peak = V_acc / T. This is the simplest possible discharge
+    estimate and assumes the entire upstream runoff arrives at each cell
+    over the same window T — true only for events much shorter than the
+    basin time-of-concentration (small basins, intense storms). Larger
+    basins will see peak attenuation that this estimator does NOT model.
+    """
+
+    values: np.ndarray
+    transform: Affine
+    crs: str
+    bbox: BoundingBox
+    duration_s: float
+    precip_source: str
+    method: str
+    units: str = "m^3/s (peak discharge)"
+
+    @property
+    def shape(self) -> tuple[int, int]:
+        """Return the (rows, cols) shape of the underlying raster."""
+        return self.values.shape  # type: ignore[return-value]
+
+    def stats(self) -> dict[str, float]:
+        """Summary statistics over finite cells."""
+        valid = self.values[~np.isnan(self.values)]
+        if valid.size == 0:
+            return {"min": 0.0, "max": 0.0, "mean": 0.0, "median": 0.0}
+        return {
+            "min": float(valid.min()),
+            "max": float(valid.max()),
+            "mean": float(valid.mean()),
+            "median": float(np.median(valid)),
+        }
+
+    def outlet_peak(self) -> float:
+        """Peak discharge at the outlet (the cell with the largest accumulated
+        volume — typically the basin outlet)."""
+        return float(np.nanmax(self.values))

floodpath/routing/utils.py

@@ -0,0 +1,60 @@
+"""Geometric helpers for the routing module."""
+
+import math
+
+import numpy as np
+from rasterio.transform import Affine
+
+from .constants import DEG_LAT_M
+
+
+def cell_areas_m2(
+    transform: Affine,
+    shape: tuple[int, int],
+    crs: str,
+) -> np.ndarray:
+    """Return per-cell areas in square metres for a raster grid.
+
+    Supports two cases:
+
+    * Geographic CRS in degrees (EPSG:4326 etc.): each row uses a
+      latitude-dependent cosine factor — pixels shrink at higher
+      latitudes. Pixel size on the X axis is multiplied by
+      `cos(latitude)` because longitude lines converge toward the poles.
+    * Projected CRS in metres: pixel area is just the absolute product
+      of the transform's scale components.
+
+    Args:
+        transform: Rasterio Affine transform of the grid.
+        shape: (rows, cols) of the raster.
+        crs: CRS of the grid as a string (e.g. ``"EPSG:4326"``).
+
+    Returns:
+        Float64 array of shape `shape` with each cell's area in m^2.
+    """
+    pixel_x = abs(transform.a)
+    pixel_y = abs(transform.e)
+    rows, _ = shape
+
+    is_geographic = ("4326" in crs) or ("WGS 84" in crs) or ("WGS84" in crs)
+    if not is_geographic:
+        # Projected CRS: pixel size is already in metres.
+        return np.full(shape, pixel_x * pixel_y, dtype=np.float64)
+
+    # Geographic CRS in degrees. Compute the latitude at every row centre
+    # using the affine transform and apply the cos(lat) longitude
+    # convergence factor.
+    row_indices = np.arange(rows)
+    # transform.f = top-edge latitude; transform.e = -pixel_y for north-up
+    # rasters. Centre latitude of row r is f + (r + 0.5) * e.
+    lat_centres = transform.f + (row_indices + 0.5) * transform.e
+    cos_lat = np.cos(np.radians(lat_centres))
+
+    cell_x_m = pixel_x * DEG_LAT_M * cos_lat  # m, per row
+    cell_y_m = pixel_y * DEG_LAT_M             # m, constant
+    areas_per_row = cell_x_m * cell_y_m        # m^2, per row
+
+    return np.broadcast_to(
+        areas_per_row[:, None],
+        shape,
+    ).astype(np.float64).copy()