Add Manning normal-depth + rainfall-driven HAND flood (2/2)

README.md

@@ -75,7 +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.routing` | runoff + flow direction (pyflwdir) + roughness + HAND | Hydrologic routing (accumulation + peak discharge) + hydraulic closure (Manning normal-depth at streams, Leopold-Maddock width) + rainfall-driven HAND flood depth |
 | `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

@@ -1,32 +1,57 @@
-"""Routing module: hydrologic + (later) hydraulic routing of runoff.
+"""Routing module: hydrologic + 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.
+Provides:
+  * accumulate_runoff(...) — flow-accumulation of SCS-CN runoff into
+    upstream-volume and peak-discharge fields (PR1: hydrologic side).
+  * compute_water_level(...) — Manning normal-depth at every stream
+    cell, given discharge + roughness + slope + Leopold-Maddock width
+    (PR2: hydraulic closure).
+  * compute_rainfall_inundation(...) — final step: broadcast per-stream
+    water levels upstream via flow direction and apply the existing HAND
+    machinery to produce a rainfall-driven flood-depth raster.
 """
 
 from .accumulate import accumulate_runoff
 from .constants import (
     DEFAULT_EVENT_DURATION_S,
     EARTH_RADIUS_M,
+    LEOPOLD_MADDOCK_A,
+    LEOPOLD_MADDOCK_B,
+    MANNING_MIN_DEPTH,
+    MANNING_MIN_SLOPE,
 )
 from .discharge import peak_discharge
+from .flood import RainfallInundationDepth, compute_rainfall_inundation
+from .manning import (
+    WaterLevelGrid,
+    compute_water_level,
+    leopold_maddock_width,
+    manning_normal_depth,
+)
 from .models import (
     AccumulatedRunoffGrid,
     DischargeGrid,
 )
-from .utils import cell_areas_m2
+from .utils import cell_areas_m2, pixel_sizes_m, slope_along_flow
 
 __all__ = [
     "AccumulatedRunoffGrid",
     "DEFAULT_EVENT_DURATION_S",
     "DischargeGrid",
     "EARTH_RADIUS_M",
+    "LEOPOLD_MADDOCK_A",
+    "LEOPOLD_MADDOCK_B",
+    "MANNING_MIN_DEPTH",
+    "MANNING_MIN_SLOPE",
+    "RainfallInundationDepth",
+    "WaterLevelGrid",
     "accumulate_runoff",
     "cell_areas_m2",
+    "compute_rainfall_inundation",
+    "compute_water_level",
+    "leopold_maddock_width",
+    "manning_normal_depth",
     "peak_discharge",
+    "pixel_sizes_m",
+    "slope_along_flow",
 ]

floodpath/routing/constants.py

@@ -21,3 +21,28 @@
 # when produced by these helpers. Future routing schemes
 # (unit hydrograph, kinematic wave) will use distinct strings.
 ROUTING_METHOD: str = "steady-state-flow-accumulation"
+
+
+# --- Channel hydraulics (Manning normal-depth + Leopold-Maddock width) -----
+#
+# Leopold & Maddock (1953) "The hydraulic geometry of stream channels and
+# some physiographic implications" — empirical relation between bankfull
+# discharge and channel width:
+#
+#     w = a * Q^b
+#
+# where a, b are calibrated regional coefficients. The values below are the
+# global natural-channel averages cited in Leopold & Maddock and refined by
+# Andreadis et al. (2013) "A simple global river bankfull width and depth
+# database" — appropriate when no local data are available. Override per
+# basin / climate when calibration data exist.
+LEOPOLD_MADDOCK_A: float = 2.5  # m * (m^3/s)^(-b)
+LEOPOLD_MADDOCK_B: float = 0.5  # dimensionless
+
+# Numerical floors. Manning's equation has a sqrt(slope) and divisions by
+# width and slope, so zero values would blow up. These caps are
+# conservatively tight — the right answer is "no flow" when slopes are
+# this small, but flagging the cell at all gives smoother boundary output
+# than a hard NaN.
+MANNING_MIN_SLOPE: float = 1e-5    # m/m — flat-as-a-pond
+MANNING_MIN_DEPTH: float = 0.05    # m — 5 cm channel-bottom default

floodpath/routing/flood.py

@@ -0,0 +1,154 @@
+"""Rainfall-driven HAND inundation: combine per-stream water levels with HAND."""
+
+from dataclasses import dataclass
+
+import numpy as np
+from rasterio.transform import Affine
+
+from floodpath.dem.models import BoundingBox
+from floodpath.hydrology.models import FlowGrid, Hand, StreamNetwork
+
+from .manning import WaterLevelGrid
+
+
+@dataclass
+class RainfallInundationDepth:
+    """Per-cell flood depth driven by per-stream-cell channel water levels.
+
+    Replaces the flat-water-level inundation used by `damage`'s static
+    `compute_inundation_depth(hand, water_level=5.0)`. The water level
+    here varies across the basin: each off-stream cell's depth is
+
+        depth = max(h_at_drainage_stream - HAND, 0)
+
+    where `h_at_drainage_stream` is the Manning normal-depth at the
+    stream cell that the off-stream cell drains to (via the flow
+    direction grid).
+    """
+
+    values: np.ndarray
+    transform: Affine
+    crs: str
+    bbox: BoundingBox
+    duration_s: float
+    discharge_source: str
+    method: str = "hand-driven-by-manning-water-levels"
+    units: str = "m (flood depth, rainfall-driven)"
+
+    @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 wet cells (depth > 0)."""
+        wet = self.values[self.values > 0]
+        if wet.size == 0:
+            return {
+                "min": 0.0, "max": 0.0, "mean": 0.0, "median": 0.0,
+                "wet_cells": 0,
+            }
+        return {
+            "min": float(wet.min()),
+            "max": float(wet.max()),
+            "mean": float(wet.mean()),
+            "median": float(np.median(wet)),
+            "wet_cells": int(wet.size),
+        }
+
+    def flooded_fraction(self) -> float:
+        """Fraction of the patch with depth > 0."""
+        if self.values.size == 0:
+            return 0.0
+        return float((self.values > 0).sum() / self.values.size)
+
+
+def _broadcast_stream_levels(
+    flow_grid: FlowGrid,
+    streams: StreamNetwork,
+    water_level: WaterLevelGrid,
+) -> np.ndarray:
+    """Project per-stream-cell water levels onto every cell that drains to them.
+
+    Uses `pyflwdir.FlwdirRaster.basins(idxs=stream_indices)` to label
+    every cell with the 1-based index of its drainage stream cell;
+    cells outside any stream's catchment (e.g. ridge cells that flow
+    off the bbox edge before hitting a stream) get basin_id=0 and water
+    level NaN.
+    """
+    flwdir = flow_grid.flwdir_raster
+    stream_flat_idxs = np.flatnonzero(streams.mask.flatten())
+
+    basins = flwdir.basins(idxs=stream_flat_idxs)  # 0 = background, else 1..N
+
+    # Build a lookup table: water level for each basin id (1-based).
+    # Position 0 is reserved for "out of basin" -> NaN.
+    h_at_streams = water_level.values.flatten()
+    h_lookup = np.full(stream_flat_idxs.size + 1, np.nan, dtype=np.float32)
+    h_lookup[1:] = h_at_streams[stream_flat_idxs]
+
+    return h_lookup[basins]
+
+
+def compute_rainfall_inundation(
+    water_level: WaterLevelGrid,
+    hand: Hand,
+    flow_grid: FlowGrid,
+    streams: StreamNetwork,
+) -> RainfallInundationDepth:
+    """Produce a rainfall-driven flood depth raster from the full chain.
+
+    For every cell:
+
+        depth = max(h_at_drainage - HAND, 0)
+
+    where `h_at_drainage` is the Manning normal-depth at the cell's
+    associated stream cell (the first downstream stream cell along the
+    flow direction). Cells outside any stream's catchment, or with HAND
+    nodata, stay dry (depth = 0).
+
+    Args:
+        water_level: WaterLevelGrid from `compute_water_level`.
+        hand: HAND raster from `compute_hand`.
+        flow_grid: Same flow grid used for HAND and routing.
+        streams: Same StreamNetwork used for HAND and water-level computation.
+
+    Returns:
+        RainfallInundationDepth carrying per-cell flood depth in m.
+
+    Raises:
+        ValueError: any of the inputs has a mismatched shape.
+    """
+    if water_level.shape != flow_grid.shape:
+        raise ValueError(
+            f"water_level shape {water_level.shape} does not match flow_grid "
+            f"shape {flow_grid.shape}"
+        )
+    if hand.shape != flow_grid.shape:
+        raise ValueError(
+            f"hand shape {hand.shape} does not match flow_grid shape "
+            f"{flow_grid.shape}"
+        )
+
+    h_field = _broadcast_stream_levels(flow_grid, streams, water_level)
+
+    # depth = max(h - HAND, 0). HAND nodata is a large negative value
+    # from pyflwdir; treat anything < 0 as nodata -> dry.
+    valid_hand = hand.values >= 0
+    valid_h = ~np.isnan(h_field)
+    diff = np.where(valid_hand & valid_h, h_field - hand.values, 0.0)
+    depth = np.maximum(diff, 0.0).astype(np.float32)
+
+    return RainfallInundationDepth(
+        values=depth,
+        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,
+        ),
+        duration_s=water_level.duration_s,
+        discharge_source=water_level.discharge_source,
+    )