Add runoff module with SCS Curve Number derivation

README.md

@@ -73,6 +73,7 @@ print(f"Total damaged built-up: {damage.values.sum():,.0f} m²")
 | `floodpath.exposure` | GHSL R2023A, WorldPop, OpenStreetMap (Overpass) | Built-up surface, population, building footprints |
 | `floodpath.landuse` | ESA WorldCover (10 m, AWS Open Data, COG) | 11-class land-cover raster (2020 v100, 2021 v200), Manning's roughness derivation |
 | `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.runoff` | NEH 630 Ch9 Table 9-1 + landuse + HSG | SCS Curve Number raster (precursor to the rainfall-runoff equation) |
 | `floodpath.damage` | JRC Huizinga 2017 + DEM/HAND/GHSL | Per-cell flood depth and damage in m² of built-up surface |
 
 ## Depth-damage curves

floodpath/runoff/__init__.py

@@ -0,0 +1,25 @@
+"""Runoff module: rainfall-runoff parameterisations and (future) routing.
+
+Currently provides the SCS Curve Number derivation that combines a
+landuse raster (`floodpath.landuse`) with a hydrologic-soil-group raster
+(`floodpath.soil`) into a per-cell CN suitable for the SCS-CN equation
+Q = (P - 0.2 * S)^2 / (P + 0.8 * S), where S = 25400/CN - 254 (mm).
+"""
+
+from .constants import (
+    CN_DEFAULT_FALLBACK,
+    CN_MAX,
+    CN_MIN,
+    WORLDCOVER_NEH_CN,
+)
+from .curve_number import compute_curve_number
+from .models import CurveNumberGrid
+
+__all__ = [
+    "CN_DEFAULT_FALLBACK",
+    "CN_MAX",
+    "CN_MIN",
+    "CurveNumberGrid",
+    "WORLDCOVER_NEH_CN",
+    "compute_curve_number",
+]

floodpath/runoff/constants.py

@@ -0,0 +1,62 @@
+"""Constants for the runoff module: NEH 630 Ch9 Curve Number tables."""
+
+from types import MappingProxyType
+
+# --- NEH 630 Chapter 9 Table 9-1 -------------------------------------------
+#
+# Per-cell Curve Numbers depend on (cover_type, hydrologic_condition,
+# treatment) AND (Hydrologic Soil Group A/B/C/D). Values below come from
+# NEH 630 Ch9 Table 9-1 (Jan 2009 release) for "AMC II" — average
+# antecedent moisture condition. AMC I/III are documented adjustments
+# that the user can apply on top.
+#
+# References:
+#   USDA NRCS National Engineering Handbook Part 630 Chapter 9
+#   "Hydrologic Soil-Cover Complexes" (Jan 2009).
+#
+# Picked NEH cover types per ESA WorldCover class follow the consensus
+# mapping used in SWAT, HEC-HMS, and the global SCS-CN literature:
+#
+#   WC 10 (tree_cover)          -> Woods, hydrologic condition GOOD
+#   WC 20 (shrubland)           -> Brush-grass, condition FAIR
+#   WC 30 (grassland)           -> Pasture/grassland, condition GOOD
+#   WC 40 (cropland)            -> Row crops, straight row, condition GOOD
+#   WC 50 (built_up)            -> Residential 1/4 acre lots
+#                                   (suburban average; users with dense
+#                                   urban patches should override with
+#                                   Commercial: A=89/B=92/C=94/D=95)
+#   WC 60 (bare_or_sparse)      -> Fallow, bare soil
+#   WC 70 (snow_and_ice)        -> Practical impervious (CN=98)
+#   WC 80 (permanent_water)     -> Open water (CN=100; no infiltration)
+#   WC 90 (herbaceous_wetland)  -> Wetlands; common practice values
+#   WC 95 (mangroves)           -> Wetlands; same as 90
+#   WC 100 (moss_and_lichen)    -> Brush, condition FAIR
+
+WORLDCOVER_NEH_CN: MappingProxyType = MappingProxyType({
+    10:  MappingProxyType({"A": 30,  "B": 55,  "C": 70,  "D": 77}),
+    20:  MappingProxyType({"A": 35,  "B": 56,  "C": 70,  "D": 77}),
+    30:  MappingProxyType({"A": 39,  "B": 61,  "C": 74,  "D": 80}),
+    40:  MappingProxyType({"A": 67,  "B": 78,  "C": 85,  "D": 89}),
+    50:  MappingProxyType({"A": 61,  "B": 75,  "C": 83,  "D": 87}),
+    60:  MappingProxyType({"A": 77,  "B": 86,  "C": 91,  "D": 94}),
+    70:  MappingProxyType({"A": 98,  "B": 98,  "C": 98,  "D": 98}),
+    80:  MappingProxyType({"A": 100, "B": 100, "C": 100, "D": 100}),
+    90:  MappingProxyType({"A": 78,  "B": 84,  "C": 88,  "D": 90}),
+    95:  MappingProxyType({"A": 78,  "B": 84,  "C": 88,  "D": 90}),
+    100: MappingProxyType({"A": 35,  "B": 56,  "C": 70,  "D": 77}),
+})
+
+# Per-cell CN range. NEH and most practitioners cap the lower end at ~30
+# (Woods, HSG A — anything more permeable than that is artificial).
+# CN=100 is the upper cap (open water / fully impervious).
+CN_MIN: int = 30
+CN_MAX: int = 100
+
+# Default value for cells where either landuse or HSG is nodata, or
+# the (landuse, HSG) pair is missing from the mapping. 70 is a typical
+# pasture/grassland B middle-of-the-table value — neutral and easy to
+# notice when it shows up unexpectedly in a CN map.
+CN_DEFAULT_FALLBACK: int = 70
+
+CN_NODATA: int = 0
+CN_UNITS: str = "Curve Number (NEH 630 Ch9, AMC II)"

floodpath/runoff/curve_number.py

@@ -0,0 +1,114 @@
+"""Combine landuse + hydrologic soil group into an SCS Curve Number raster."""
+
+import numpy as np
+import rasterio.warp
+from rasterio.enums import Resampling
+
+from floodpath.dem.models import BoundingBox
+from floodpath.landuse.models import LanduseGrid
+from floodpath.soil.constants import HSG_TO_INT
+from floodpath.soil.models import HSGGrid
+
+from .constants import (
+    CN_DEFAULT_FALLBACK,
+    CN_NODATA,
+    WORLDCOVER_NEH_CN,
+)
+from .models import CurveNumberGrid
+
+
+def _upsample_hsg_to_landuse_grid(
+    hsg: HSGGrid,
+    landuse: LanduseGrid,
+) -> np.ndarray:
+    """Project an HSG raster onto the landuse grid via nearest-neighbour.
+
+    HSG is categorical (1=A, 2=B, 3=C, 4=D, 0=nodata) so any resampling
+    other than nearest-neighbour would produce meaningless intermediate
+    codes. WorldCover (10 m) is much finer than SoilGrids (~250 m), so
+    every fine cell inherits its parent coarse cell's HSG.
+
+    Used by `compute_curve_number`. Both grids are assumed to be in the
+    same CRS and to overlap geographically (the caller is responsible
+    for using the same `lat, lon, buffer_deg` when fetching).
+    """
+    out = np.zeros(landuse.values.shape, dtype=np.uint8)
+    rasterio.warp.reproject(
+        source=hsg.values,
+        destination=out,
+        src_transform=hsg.transform,
+        src_crs=hsg.crs,
+        dst_transform=landuse.transform,
+        dst_crs=landuse.crs,
+        resampling=Resampling.nearest,
+    )
+    return out
+
+
+def compute_curve_number(
+    landuse: LanduseGrid,
+    hsg: HSGGrid,
+    mapping: dict[int, dict[str, int]] | None = None,
+    fallback: int = CN_DEFAULT_FALLBACK,
+) -> CurveNumberGrid:
+    """Derive a per-cell Curve Number raster from landuse + HSG.
+
+    Each cell's CN is read from `mapping[landuse_code][hsg_letter]`,
+    with NEH 630 Ch9 Table 9-1 typical values as the default. Cells
+    where either landuse or HSG is nodata, or where the (landuse, HSG)
+    pair is missing from the mapping, get `fallback`.
+
+    The output raster shares the landuse grid's resolution and georef.
+    HSG (typically ~250 m from SoilGrids) is upsampled to the landuse
+    grid (typically 10 m for WorldCover) via nearest-neighbour — every
+    landuse cell inherits its parent HSG cell's group.
+
+    Args:
+        landuse: Categorical land-cover grid (e.g. `floodpath.landuse`).
+        hsg: NEH 630 Ch7 hydrologic-soil-group grid (e.g. `floodpath.soil`).
+        mapping: Two-level dict landuse_code -> HSG_letter -> CN. Defaults
+            to `WORLDCOVER_NEH_CN` (NEH 630 Ch9 + ESA WorldCover).
+        fallback: CN value to use for cells the mapping doesn't cover.
+            Defaults to 70 (NEH-typical pasture/grassland B).
+
+    Returns:
+        CurveNumberGrid sharing the landuse grid's georef.
+    """
+    table = WORLDCOVER_NEH_CN if mapping is None else mapping
+
+    hsg_at_lu = _upsample_hsg_to_landuse_grid(hsg, landuse)
+
+    cn = np.full(landuse.values.shape, fallback, dtype=np.uint8)
+    # Walk every (landuse_code, HSG_letter) pair in the mapping and
+    # apply its CN to the matching cells.
+    for lu_code, hsg_to_cn in table.items():
+        lu_mask = landuse.values == lu_code
+        if not lu_mask.any():
+            continue
+        for letter, cn_value in hsg_to_cn.items():
+            hsg_int = HSG_TO_INT.get(letter)
+            if hsg_int is None:
+                raise ValueError(
+                    f"mapping refers to HSG letter {letter!r} which is not "
+                    f"one of A/B/C/D"
+                )
+            cell_mask = lu_mask & (hsg_at_lu == hsg_int)
+            cn[cell_mask] = int(cn_value)
+
+    # Cells where either source is nodata: tag CN as nodata (0). Landuse
+    # nodata is implicit (no class code matches). HSG nodata = 0 in the
+    # upsampled raster.
+    nodata_mask = (hsg_at_lu == 0)
+    cn[nodata_mask] = CN_NODATA
+
+    return CurveNumberGrid(
+        values=cn,
+        transform=landuse.transform,
+        crs=landuse.crs,
+        bbox=BoundingBox(
+            min_lon=landuse.bbox.min_lon,
+            min_lat=landuse.bbox.min_lat,
+            max_lon=landuse.bbox.max_lon,
+            max_lat=landuse.bbox.max_lat,
+        ),
+    )

floodpath/runoff/models.py

@@ -0,0 +1,67 @@
+"""Dataclasses for the runoff module."""
+
+from dataclasses import dataclass
+
+import numpy as np
+from rasterio.transform import Affine
+
+from floodpath.dem.models import BoundingBox
+
+from .constants import CN_NODATA, CN_UNITS
+
+
+@dataclass
+class CurveNumberGrid:
+    """Per-cell SCS Curve Number raster.
+
+    Each cell carries an integer CN in [CN_MIN, CN_MAX]; CN_NODATA (0)
+    marks cells where either the landuse or hydrologic-soil-group input
+    was missing.
+
+    The Curve Number couples directly to the SCS-CN runoff equation:
+        S = 25400/CN - 254          (mm, potential maximum retention)
+        Q = (P - 0.2*S)^2 / (P + 0.8*S)  if P > 0.2*S, else 0
+    where P is precipitation depth (mm) and Q is runoff depth (mm).
+    """
+
+    values: np.ndarray
+    transform: Affine
+    crs: str
+    bbox: BoundingBox
+    units: str = CN_UNITS
+
+    @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]:
+        """Return summary statistics over classified cells (CN > 0).
+
+        Cells with CN=0 (nodata) are excluded so a few unmapped pixels
+        don't drag the mean down.
+        """
+        valid = self.values[self.values > 0]
+        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 potential_retention_mm(self) -> np.ndarray:
+        """Return S = 25400/CN - 254 (mm) at every cell.
+
+        S is the potential maximum retention used by the SCS-CN runoff
+        equation. Cells with CN=0 (nodata) get NaN. CN=100 (open water,
+        fully impervious) gives S=0 — runoff equals precipitation.
+        """
+        out = np.full(self.values.shape, np.nan, dtype=np.float32)
+        mask = self.values > 0
+        cn = self.values[mask].astype(np.float32)
+        # Clip to avoid division-by-zero blow-up at CN exactly 0; CN_NODATA
+        # is already excluded by the mask but be defensive about float math.
+        out[mask] = 25400.0 / cn - 254.0
+        return out

tests/runoff/test_constants.py

@@ -0,0 +1,67 @@
+"""Sanity tests for the runoff constants."""
+
+from floodpath.landuse import WORLDCOVER_CLASSES
+from floodpath.runoff import (
+    CN_DEFAULT_FALLBACK,
+    CN_MAX,
+    CN_MIN,
+    WORLDCOVER_NEH_CN,
+)
+
+
+class TestWorldcoverNehCn:
+    def test_covers_every_worldcover_class(self) -> None:
+        # Every published WorldCover class needs a CN entry; otherwise
+        # those cells would silently fall through to the fallback.
+        assert set(WORLDCOVER_NEH_CN.keys()) == set(WORLDCOVER_CLASSES.keys())
+
+    def test_each_entry_has_all_four_hsg_letters(self) -> None:
+        for wc_code, hsg_table in WORLDCOVER_NEH_CN.items():
+            assert set(hsg_table.keys()) == {"A", "B", "C", "D"}, (
+                f"WC {wc_code} missing one of A/B/C/D"
+            )
+
+    def test_cn_values_in_valid_range(self) -> None:
+        for wc_code, hsg_table in WORLDCOVER_NEH_CN.items():
+            for letter, cn in hsg_table.items():
+                assert CN_MIN <= cn <= CN_MAX, (
+                    f"WC {wc_code}/{letter} = {cn} outside [{CN_MIN},{CN_MAX}]"
+                )
+
+    def test_cn_increases_with_hsg_letter_for_each_class(self) -> None:
+        # NEH 630 Ch9: for any cover type, CN(A) <= CN(B) <= CN(C) <= CN(D),
+        # because more clay -> less infiltration -> more runoff -> higher CN.
+        for wc_code, hsg_table in WORLDCOVER_NEH_CN.items():
+            cns = [hsg_table[letter] for letter in ("A", "B", "C", "D")]
+            assert cns == sorted(cns), (
+                f"WC {wc_code} CN not monotone in HSG: {cns}"
+            )
+
+    def test_water_is_cn_100(self) -> None:
+        # Open water (no infiltration capacity) -> CN=100 by NEH convention.
+        assert WORLDCOVER_NEH_CN[80] == {"A": 100, "B": 100, "C": 100, "D": 100}
+
+    def test_snow_ice_is_practical_impervious(self) -> None:
+        # Frozen surfaces don't infiltrate; CN=98 is the typical NEH value.
+        assert WORLDCOVER_NEH_CN[70] == {"A": 98, "B": 98, "C": 98, "D": 98}
+
+    def test_cropland_on_d_is_high(self) -> None:
+        # Cropland (row crops, straight row, good condition) on HSG D is
+        # NEH 630 Ch9 Table 9-1's canonical 89 — the floodpath default
+        # for cropland on the most-runoff-prone soils.
+        assert WORLDCOVER_NEH_CN[40]["D"] == 89
+
+    def test_tree_cover_on_a_is_low(self) -> None:
+        # Woods, good condition, HSG A is NEH 630 Ch9's canonical 30
+        # (the table's lowest non-pavement value).
+        assert WORLDCOVER_NEH_CN[10]["A"] == 30
+
+
+class TestFallback:
+    def test_fallback_is_in_valid_range(self) -> None:
+        assert CN_MIN <= CN_DEFAULT_FALLBACK <= CN_MAX
+
+    def test_fallback_is_neutral_middle(self) -> None:
+        # 70 corresponds to NEH-typical pasture/grassland B — a sensible
+        # middle-of-table neutral when nothing else is known.
+        assert CN_DEFAULT_FALLBACK == 70