Add precip module + SCS-CN equation: end-to-end rainfall->runoff pipeline

README.md

@@ -73,7 +73,8 @@ 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.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.damage` | JRC Huizinga 2017 + DEM/HAND/GHSL | Per-cell flood depth and damage in m² of built-up surface |
 
 ## Depth-damage curves

floodpath/precip/__init__.py

@@ -0,0 +1,16 @@
+"""Precipitation module: precipitation depth rasters from any source.
+
+For now, only `uniform_precip` (synthetic spatially-flat depth) is exposed —
+the consumer (`floodpath.runoff.apply_scs_cn`) doesn't care whether the
+precip came from synthetic input, ERA5, IMERG, CHIRPS, or anywhere else.
+Future fetchers slot in here.
+"""
+
+from .models import PrecipGrid
+from .uniform import uniform_precip, uniform_precip_like
+
+__all__ = [
+    "PrecipGrid",
+    "uniform_precip",
+    "uniform_precip_like",
+]

floodpath/precip/models.py

@@ -0,0 +1,54 @@
+"""Dataclasses for the precip module."""
+
+from dataclasses import dataclass
+
+import numpy as np
+from rasterio.transform import Affine
+
+from floodpath.dem.models import BoundingBox
+
+PRECIP_UNITS: str = "mm (precipitation depth)"
+
+
+@dataclass
+class PrecipGrid:
+    """Per-cell precipitation depth (mm) over an event window.
+
+    `values` is `mm` of liquid-equivalent precipitation accumulated over
+    the period described by `(start, end)`. Synthetic uniform inputs
+    leave the period as `None` (event-agnostic). Real fetchers
+    (ERA5/IMERG/CHIRPS) record the actual UTC window.
+
+    `source` is a free-form string identifying where the data came from
+    so downstream consumers can attribute / cite — e.g. `"uniform"`,
+    `"ERA5"`, `"IMERG-Final"`.
+    """
+
+    values: np.ndarray
+    transform: Affine
+    crs: str
+    bbox: BoundingBox
+    source: str
+    start: str | None = None
+    end: str | None = None
+    units: str = PRECIP_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 total_volume_m3(self, cell_area_m2: float) -> float:
+        """Bulk volume of liquid water over the patch, in m^3.
+
+        Args:
+            cell_area_m2: Area of one cell in m^2 (the caller knows their
+                grid's cell area; we don't infer it from the WGS84
+                transform because that would require a per-cell cosine
+                latitude correction).
+
+        Returns:
+            Total water volume in cubic metres.
+        """
+        # mm * m^2 = L, divide by 1000 -> m^3.
+        return float(np.nansum(self.values) * cell_area_m2 / 1000.0)

floodpath/precip/uniform.py

@@ -0,0 +1,84 @@
+"""Synthetic uniform-depth precipitation grids."""
+
+from typing import Protocol
+
+import numpy as np
+from rasterio.transform import Affine
+
+from floodpath.dem.models import BoundingBox
+
+from .models import PrecipGrid
+
+
+class _Griddable(Protocol):
+    """Anything with a `shape`, `transform`, `crs`, `bbox` quartet — every
+    floodpath grid dataclass satisfies this."""
+
+    shape: tuple[int, int]
+    transform: Affine
+    crs: str
+    bbox: BoundingBox
+
+
+def uniform_precip(
+    bbox: BoundingBox,
+    transform: Affine,
+    shape: tuple[int, int],
+    depth_mm: float,
+    crs: str = "EPSG:4326",
+) -> PrecipGrid:
+    """Build a spatially-uniform precipitation grid of given depth.
+
+    Useful as the simplest possible input to test the rainfall-runoff
+    pipeline before plugging in a real fetcher (ERA5 / IMERG / CHIRPS).
+
+    Args:
+        bbox: Geographic bounding box matching the target grid.
+        transform: Rasterio Affine transform of the target grid.
+        shape: (rows, cols) of the target grid.
+        depth_mm: Constant precipitation depth applied to every cell, in mm.
+        crs: CRS of the grid; defaults to WGS84.
+
+    Returns:
+        PrecipGrid filled with `depth_mm` everywhere.
+
+    Raises:
+        ValueError: `depth_mm` is negative.
+    """
+    if depth_mm < 0:
+        raise ValueError(f"depth_mm must be >= 0, got {depth_mm}")
+    values = np.full(shape, float(depth_mm), dtype=np.float32)
+    return PrecipGrid(
+        values=values,
+        transform=transform,
+        crs=crs,
+        bbox=bbox,
+        source="uniform",
+    )
+
+
+def uniform_precip_like(
+    grid: _Griddable,
+    depth_mm: float,
+) -> PrecipGrid:
+    """Build a uniform PrecipGrid that overlays another grid exactly.
+
+    Convenience wrapper — when you already have a CN/HSG/landuse grid and
+    want a precip grid that aligns one-to-one, this avoids passing the
+    georef arguments by hand.
+
+    Args:
+        grid: Any grid dataclass exposing `shape`, `transform`, `crs`,
+            `bbox` (CurveNumberGrid, LanduseGrid, HSGGrid, ...).
+        depth_mm: Constant precipitation depth applied everywhere, in mm.
+
+    Returns:
+        PrecipGrid sharing the input grid's georef.
+    """
+    return uniform_precip(
+        bbox=grid.bbox,
+        transform=grid.transform,
+        shape=grid.shape,
+        depth_mm=depth_mm,
+        crs=grid.crs,
+    )

floodpath/runoff/__init__.py

@@ -1,9 +1,11 @@
 """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).
+Provides:
+  * compute_curve_number(landuse, hsg) -> CurveNumberGrid (NEH 630 Ch9)
+  * apply_scs_cn(cn, precip) -> RunoffGrid (the SCS-CN equation itself)
+
+These together turn a precipitation depth raster (any source — synthetic
+uniform, ERA5, IMERG, ...) into per-cell runoff depth, ready for routing.
 """
 
 from .constants import (
@@ -13,13 +15,16 @@
     WORLDCOVER_NEH_CN,
 )
 from .curve_number import compute_curve_number
-from .models import CurveNumberGrid
+from .models import CurveNumberGrid, RunoffGrid
+from .runoff import apply_scs_cn
 
 __all__ = [
     "CN_DEFAULT_FALLBACK",
     "CN_MAX",
     "CN_MIN",
     "CurveNumberGrid",
+    "RunoffGrid",
     "WORLDCOVER_NEH_CN",
+    "apply_scs_cn",
     "compute_curve_number",
 ]

floodpath/runoff/models.py

@@ -9,6 +9,8 @@
 
 from .constants import CN_NODATA, CN_UNITS
 
+RUNOFF_UNITS: str = "mm (runoff depth, SCS-CN)"
+
 
 @dataclass
 class CurveNumberGrid:
@@ -65,3 +67,41 @@ def potential_retention_mm(self) -> np.ndarray:
         # is already excluded by the mask but be defensive about float math.
         out[mask] = 25400.0 / cn - 254.0
         return out
+
+
+@dataclass
+class RunoffGrid:
+    """Per-cell runoff depth Q (mm) from the SCS-CN equation.
+
+    Pair-companion to PrecipGrid. The same grid that carried `mm` of
+    rainfall now carries `mm` of effective runoff depth — the share of
+    precipitation that doesn't infiltrate. Cells where CN was nodata
+    inherit NaN.
+
+    `precip_source` records where the input precipitation came from
+    (`"uniform"`, `"ERA5"`, ...) so consumers can attribute results.
+    """
+
+    values: np.ndarray
+    transform: Affine
+    crs: str
+    bbox: BoundingBox
+    precip_source: str
+    units: str = RUNOFF_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]:
+        """Summary statistics over cells with a finite Q."""
+        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)),
+        }

floodpath/runoff/runoff.py

@@ -0,0 +1,68 @@
+"""Apply the SCS Curve Number runoff equation to a precip + CN pair."""
+
+import numpy as np
+
+from floodpath.precip.models import PrecipGrid
+
+from .models import CurveNumberGrid, RunoffGrid
+
+
+def apply_scs_cn(
+    cn: CurveNumberGrid,
+    precip: PrecipGrid,
+    initial_abstraction_ratio: float = 0.2,
+) -> RunoffGrid:
+    """Compute per-cell runoff depth Q (mm) via the SCS-CN equation.
+
+    For each cell:
+
+        S  = 25400 / CN - 254          (mm, potential maximum retention)
+        Ia = lambda * S                (mm, initial abstraction)
+        Q  = (P - Ia)^2 / (P - Ia + S) if P > Ia, else 0
+
+    where `lambda` is the initial-abstraction ratio (NEH 630 Ch10 standard
+    value: 0.2). Cells with CN nodata propagate NaN.
+
+    Args:
+        cn: Curve Number raster (any source).
+        precip: Precipitation depth raster, in mm. Must share the CN grid's
+            shape (call uniform_precip_like(cn, depth_mm) for the
+            simplest case).
+        initial_abstraction_ratio: lambda in Ia = lambda * S. Defaults to
+            0.2 per NEH 630 Ch10. Some authors recommend 0.05 for newer
+            calibrations — pass it explicitly if you want that.
+
+    Returns:
+        RunoffGrid with per-cell Q in mm. Same georef as the CN raster.
+
+    Raises:
+        ValueError: precip and CN have mismatched shapes.
+    """
+    if precip.shape != cn.shape:
+        raise ValueError(
+            f"precip shape {precip.shape} does not match CN shape {cn.shape}; "
+            "use uniform_precip_like(cn, depth_mm) or reproject the precip grid."
+        )
+
+    S = cn.potential_retention_mm()  # mm; NaN where CN is nodata
+    P = precip.values.astype(np.float32)
+    Ia = float(initial_abstraction_ratio) * S
+
+    with np.errstate(invalid="ignore", divide="ignore"):
+        runoff = np.where(
+            P > Ia,
+            (P - Ia) ** 2 / (P - Ia + S),
+            0.0,
+        )
+
+    # Anywhere S is NaN (CN nodata), preserve NaN in Q so downstream
+    # consumers can tell unmodelled cells apart from "no runoff".
+    runoff = np.where(np.isnan(S), np.nan, runoff)
+
+    return RunoffGrid(
+        values=runoff.astype(np.float32),
+        transform=cn.transform,
+        crs=cn.crs,
+        bbox=cn.bbox,
+        precip_source=precip.source,
+    )

tests/precip/test_models.py

@@ -0,0 +1,37 @@
+"""Unit tests for precip.models (PrecipGrid)."""
+
+import numpy as np
+import pytest
+from rasterio.transform import Affine
+
+from floodpath.dem.models import BoundingBox
+from floodpath.precip import PrecipGrid
+
+
+def _toy(values: np.ndarray) -> PrecipGrid:
+    return PrecipGrid(
+        values=values.astype(np.float32),
+        transform=Affine.translation(0, 0) * Affine.scale(0.001, -0.001),
+        crs="EPSG:4326",
+        bbox=BoundingBox(0.0, 0.0, 0.001 * values.shape[1], 0.001 * values.shape[0]),
+        source="test",
+    )
+
+
+class TestPrecipGrid:
+    def test_shape(self) -> None:
+        p = _toy(np.full((4, 4), 10.0))
+        assert p.shape == (4, 4)
+
+    def test_total_volume_m3_uniform(self) -> None:
+        # 100 mm over 1 m^2 -> 0.1 m^3 = 100 L per cell.
+        # 4 cells = 0.4 m^3 total.
+        p = _toy(np.full((2, 2), 100.0))
+        assert p.total_volume_m3(cell_area_m2=1.0) == pytest.approx(0.4)
+
+    def test_total_volume_m3_nan_excluded(self) -> None:
+        # nansum should ignore NaN cells.
+        values = np.array([[100.0, np.nan], [50.0, 100.0]])
+        p = _toy(values)
+        # 250 mm-cells * 1 m^2 / 1000 = 0.25 m^3
+        assert p.total_volume_m3(cell_area_m2=1.0) == pytest.approx(0.25)