From 9b4ef400a15684192eab6f2762679305413e79ea Mon Sep 17 00:00:00 2001 From: rehsani Date: Mon, 4 May 2026 11:10:04 -0700 Subject: [PATCH] Add runoff module with SCS Curve Number derivation MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit floodpath.runoff combines ESA WorldCover landuse with NEH 630 Ch7 hydrologic soil group into the per-cell SCS Curve Number raster needed by the SCS-CN runoff equation Q = (P - 0.2*S)^2 / (P + 0.8*S), S = 25400/CN - 254 (mm). This is the second-to-last piece of the rainfall-runoff chain — only the rainfall fetcher and the equation itself remain. - floodpath/runoff/curve_number.py: compute_curve_number(landuse, hsg) walks each (WorldCover code, HSG letter) pair and applies the table in WORLDCOVER_NEH_CN. SoilGrids HSG (~250 m) is upsampled to the WorldCover grid (10 m) via rasterio.warp.reproject with nearest- neighbour resampling (categorical data — averaging would be wrong). - floodpath/runoff/constants.py: WORLDCOVER_NEH_CN — the 11x4 lookup table (NEH 630 Ch9 Table 9-1, AMC II) mapping every WorldCover class to a CN per HSG letter. Mappings cite the chosen NEH cover type (Woods/Brush/Pasture/Row crops/Residential 1/4 acre/Fallow/...). - floodpath/runoff/models.py: CurveNumberGrid carries the uint8 CN raster, stats() that ignore nodata, and potential_retention_mm() returning S = 25400/CN - 254 ready to plug into the SCS-CN equation. - Tests: 29 new (constants table invariants, monotonicity vs HSG, S formula, georef inheritance, fallback behaviour, custom mapping). Pinned values against the existing Robit Bata WorldCover + soil fixtures (no new fixture needed): mean CN = 87.72, dominant class CN = 89 (cropland on D), min/max 77/100. --- README.md | 1 + floodpath/runoff/__init__.py | 25 +++++ floodpath/runoff/constants.py | 62 ++++++++++ floodpath/runoff/curve_number.py | 114 +++++++++++++++++++ floodpath/runoff/models.py | 67 +++++++++++ tests/runoff/__init__.py | 0 tests/runoff/test_constants.py | 67 +++++++++++ tests/runoff/test_curve_number.py | 180 ++++++++++++++++++++++++++++++ tests/runoff/test_models.py | 72 ++++++++++++ 9 files changed, 588 insertions(+) create mode 100644 floodpath/runoff/__init__.py create mode 100644 floodpath/runoff/constants.py create mode 100644 floodpath/runoff/curve_number.py create mode 100644 floodpath/runoff/models.py create mode 100644 tests/runoff/__init__.py create mode 100644 tests/runoff/test_constants.py create mode 100644 tests/runoff/test_curve_number.py create mode 100644 tests/runoff/test_models.py diff --git a/README.md b/README.md index e3fe758..fdf7fe4 100644 --- a/README.md +++ b/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 diff --git a/floodpath/runoff/__init__.py b/floodpath/runoff/__init__.py new file mode 100644 index 0000000..840913f --- /dev/null +++ b/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", +] diff --git a/floodpath/runoff/constants.py b/floodpath/runoff/constants.py new file mode 100644 index 0000000..31afe0d --- /dev/null +++ b/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)" diff --git a/floodpath/runoff/curve_number.py b/floodpath/runoff/curve_number.py new file mode 100644 index 0000000..8495e4e --- /dev/null +++ b/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, + ), + ) diff --git a/floodpath/runoff/models.py b/floodpath/runoff/models.py new file mode 100644 index 0000000..42b29c9 --- /dev/null +++ b/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 diff --git a/tests/runoff/__init__.py b/tests/runoff/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/tests/runoff/test_constants.py b/tests/runoff/test_constants.py new file mode 100644 index 0000000..04d10b5 --- /dev/null +++ b/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 diff --git a/tests/runoff/test_curve_number.py b/tests/runoff/test_curve_number.py new file mode 100644 index 0000000..4b494dc --- /dev/null +++ b/tests/runoff/test_curve_number.py @@ -0,0 +1,180 @@ +"""Unit tests for runoff.curve_number.""" + +import numpy as np +import pytest +from rasterio.transform import Affine + +from floodpath.dem.models import BoundingBox +from floodpath.landuse import WORLDCOVER_CLASSES +from floodpath.landuse.models import LanduseGrid +from floodpath.runoff import ( + CN_DEFAULT_FALLBACK, + WORLDCOVER_NEH_CN, + compute_curve_number, +) +from floodpath.soil import HSG_TO_INT, texture_to_hsg +from floodpath.soil.models import HSGGrid, TextureGrid + + +def _aligned_landuse_hsg_pair( + lu_values: np.ndarray, + hsg_values: np.ndarray, +) -> tuple[LanduseGrid, HSGGrid]: + """Build a co-located LanduseGrid + HSGGrid sharing the same bbox/grid. + + Used to test compute_curve_number with deterministic inputs (no SoilGrids + upsampling happens because the grids are already aligned). + """ + transform = Affine.translation(0, 0) * Affine.scale(0.001, -0.001) + bbox = BoundingBox(0.0, 0.0, 0.001 * lu_values.shape[1], 0.001 * lu_values.shape[0]) + landuse = LanduseGrid( + values=lu_values.astype(np.uint8), + transform=transform, + crs="EPSG:4326", + bbox=bbox, + epoch=2021, + classes=dict(WORLDCOVER_CLASSES), + units="ESA WorldCover class code", + ) + hsg = HSGGrid( + values=hsg_values.astype(np.uint8), + transform=transform, + crs="EPSG:4326", + bbox=bbox, + depths=("0-5cm", "5-15cm", "15-30cm"), + ) + return landuse, hsg + + +class TestComputeCurveNumberAligned: + def test_canonical_quartet(self) -> None: + # 2x2 grid with one (landuse, HSG) pair per cell. Each cell should + # pick up exactly the corresponding entry from WORLDCOVER_NEH_CN. + # WC 40 (cropland) x HSG D -> 89 + # WC 10 (tree) x HSG A -> 30 + # WC 30 (grass) x HSG B -> 61 + # WC 50 (built) x HSG C -> 83 + landuse, hsg = _aligned_landuse_hsg_pair( + lu_values=np.array([[40, 10], [30, 50]]), + hsg_values=np.array([ + [HSG_TO_INT["D"], HSG_TO_INT["A"]], + [HSG_TO_INT["B"], HSG_TO_INT["C"]], + ]), + ) + cn = compute_curve_number(landuse, hsg) + np.testing.assert_array_equal( + cn.values, + np.array([[89, 30], [61, 83]], dtype=np.uint8), + ) + + def test_nodata_hsg_yields_cn_zero(self) -> None: + # Cells where HSG is nodata (0) must end up CN=0, NOT the fallback. + landuse, hsg = _aligned_landuse_hsg_pair( + lu_values=np.array([[40, 40]]), + hsg_values=np.array([[HSG_TO_INT["D"], 0]]), + ) + cn = compute_curve_number(landuse, hsg) + assert cn.values[0, 0] == 89 + assert cn.values[0, 1] == 0 + + def test_unmapped_landuse_uses_fallback(self) -> None: + # Landuse class 200 is not in WORLDCOVER_NEH_CN -> fallback. + landuse, hsg = _aligned_landuse_hsg_pair( + lu_values=np.array([[200]]), + hsg_values=np.array([[HSG_TO_INT["B"]]]), + ) + cn = compute_curve_number(landuse, hsg) + assert cn.values[0, 0] == CN_DEFAULT_FALLBACK + + def test_custom_fallback_applied(self) -> None: + landuse, hsg = _aligned_landuse_hsg_pair( + lu_values=np.array([[200]]), + hsg_values=np.array([[HSG_TO_INT["B"]]]), + ) + cn = compute_curve_number(landuse, hsg, fallback=42) + assert cn.values[0, 0] == 42 + + def test_custom_mapping(self) -> None: + # Override cropland on HSG B to a non-default value. + landuse, hsg = _aligned_landuse_hsg_pair( + lu_values=np.array([[40]]), + hsg_values=np.array([[HSG_TO_INT["B"]]]), + ) + cn = compute_curve_number( + landuse, hsg, + mapping={40: {"A": 50, "B": 99, "C": 60, "D": 70}}, + ) + assert cn.values[0, 0] == 99 + + def test_invalid_hsg_letter_in_mapping_raises(self) -> None: + landuse, hsg = _aligned_landuse_hsg_pair( + lu_values=np.array([[40]]), + hsg_values=np.array([[HSG_TO_INT["D"]]]), + ) + with pytest.raises(ValueError, match="not one of A/B/C/D"): + compute_curve_number( + landuse, hsg, + mapping={40: {"X": 50}}, + ) + + def test_georef_inherits_landuse(self) -> None: + landuse, hsg = _aligned_landuse_hsg_pair( + lu_values=np.array([[40]]), + hsg_values=np.array([[HSG_TO_INT["D"]]]), + ) + cn = compute_curve_number(landuse, hsg) + assert cn.transform == landuse.transform + assert cn.crs == landuse.crs + + +class TestRobitBataCurveNumber: + """Pinned values from running compute_curve_number on the committed + Robit Bata WorldCover + soil fixtures.""" + + @pytest.fixture + def robit_bata_cn( + self, + robit_bata_worldcover: LanduseGrid, + robit_bata_soil: TextureGrid, + ) -> "CurveNumberGrid": + from floodpath.runoff import CurveNumberGrid as _CN # noqa: F401 + hsg = texture_to_hsg(robit_bata_soil) + return compute_curve_number(robit_bata_worldcover, hsg) + + def test_shape_matches_landuse(self, robit_bata_cn) -> None: + # CN raster inherits the WorldCover (10 m) resolution. + assert robit_bata_cn.shape == (900, 900) + + def test_cn_histogram_pinned(self, robit_bata_cn) -> None: + # Pinned: with the Robit Bata fixtures' all-D HSG and dominant + # cropland cover, the CN values follow WORLDCOVER_NEH_CN[*]['D'] + # plus a few hundred nodata cells from the single HSG nodata + # pixel projected onto the 27x finer landuse grid. + unique, counts = np.unique(robit_bata_cn.values, return_counts=True) + hist = {int(u): int(c) for u, c in zip(unique, counts)} + # Total cells: 900 x 900 = 810,000 + assert sum(hist.values()) == 810_000 + # CN values must all be the HSG-D column from the default mapping + # (or 0 nodata, or 100 for water). + d_column = {WORLDCOVER_NEH_CN[wc]["D"] for wc in WORLDCOVER_NEH_CN} + observed_cn = set(hist.keys()) - {0} + assert observed_cn.issubset(d_column) + # The dominant value must be cropland-on-D = 89. + dominant = max(hist, key=hist.get) + assert dominant == 89 + + def test_cn_stats_pinned(self, robit_bata_cn) -> None: + # Pinned mean (over classified cells only): ~87.7 — high CN + # consistent with clay-rich soils + cropland-dominant cover. + s = robit_bata_cn.stats() + assert s["mean"] == pytest.approx(87.72, abs=0.05) + assert s["min"] == 77 # tree_cover on D + assert s["max"] == 100 # the single water cell + assert s["median"] == 89 # cropland on D + + def test_potential_retention_makes_sense(self, robit_bata_cn) -> None: + # S range follows from CN range. CN=100 -> S=0; CN=77 -> S~76 mm. + S = robit_bata_cn.potential_retention_mm() + valid_S = S[~np.isnan(S)] + assert valid_S.min() == pytest.approx(0.0, abs=0.01) + assert valid_S.max() == pytest.approx(75.9, abs=0.1) diff --git a/tests/runoff/test_models.py b/tests/runoff/test_models.py new file mode 100644 index 0000000..7a559cc --- /dev/null +++ b/tests/runoff/test_models.py @@ -0,0 +1,72 @@ +"""Unit tests for runoff.models (CurveNumberGrid).""" + +import numpy as np +import pytest +from rasterio.transform import Affine + +from floodpath.dem.models import BoundingBox +from floodpath.runoff import CurveNumberGrid + + +def _toy_cn(values: np.ndarray) -> CurveNumberGrid: + return CurveNumberGrid( + values=values.astype(np.uint8), + 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]), + ) + + +class TestCurveNumberGrid: + def test_shape(self) -> None: + cn = _toy_cn(np.full((4, 4), 80, dtype=np.uint8)) + assert cn.shape == (4, 4) + + def test_stats_excludes_nodata(self) -> None: + # Mix of CN=80 (valid) and CN=0 (nodata). + cn = _toy_cn(np.array([[80, 80, 0], [80, 0, 0]])) + s = cn.stats() + # Nodata must not pull mean down — only the three CN=80 cells count. + assert s["min"] == pytest.approx(80) + assert s["max"] == pytest.approx(80) + assert s["mean"] == pytest.approx(80) + + def test_stats_empty_when_all_nodata(self) -> None: + cn = _toy_cn(np.zeros((3, 3), dtype=np.uint8)) + s = cn.stats() + assert s == {"min": 0.0, "max": 0.0, "mean": 0.0, "median": 0.0} + + +class TestPotentialRetention: + def test_cn_100_gives_zero_retention(self) -> None: + # CN=100 (open water, fully impervious): S = 25400/100 - 254 = 0. + cn = _toy_cn(np.array([[100]])) + S = cn.potential_retention_mm() + assert S[0, 0] == pytest.approx(0.0) + + def test_cn_77_gives_75_87_mm(self) -> None: + # CN=77 (Woods on D): S = 25400/77 - 254 = 75.87 mm. Pinned to the + # exact NEH value so a refactor that mis-orders 254/25400 will + # blow up here. + cn = _toy_cn(np.array([[77]])) + S = cn.potential_retention_mm() + assert S[0, 0] == pytest.approx(75.87, abs=0.01) + + def test_cn_50_gives_254_mm(self) -> None: + # CN=50 (mid-table): S = 25400/50 - 254 = 254 mm exactly. + cn = _toy_cn(np.array([[50]])) + S = cn.potential_retention_mm() + assert S[0, 0] == pytest.approx(254.0, abs=0.001) + + def test_nodata_is_nan(self) -> None: + cn = _toy_cn(np.array([[80, 0]])) + S = cn.potential_retention_mm() + assert S[0, 0] == pytest.approx(63.5, abs=0.01) + assert np.isnan(S[0, 1]) + + def test_higher_cn_gives_lower_retention(self) -> None: + # Monotone: as CN increases, S decreases (less retention = + # more runoff). + cn = _toy_cn(np.array([[60, 80, 95]])) + S = cn.potential_retention_mm() + assert S[0, 0] > S[0, 1] > S[0, 2]