From ae87002594d33578f473487d1f5e2ccbd36db931 Mon Sep 17 00:00:00 2001 From: rehsani Date: Mon, 4 May 2026 11:55:20 -0700 Subject: [PATCH] Add routing module: steady-state runoff accumulation + peak discharge MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit floodpath.routing closes the rainfall-runoff loop's hydrologic-routing half. Per-cell SCS-CN runoff (mm) is reprojected onto the DEM-derived flow direction grid, converted to volume (m^3) using lat-aware cell areas, and accumulated downstream via pyflwdir's accuflux. Dividing the accumulated volume by an event duration yields per-cell peak discharge (m^3/s). PR1 of a planned two-PR routing sequence; the follow-up adds Manning normal-depth at stream cells to close the hydraulic side. - floodpath/routing/utils.py: cell_areas_m2(transform, shape, crs) handles WGS84 cos(lat) longitude convergence for geographic CRSes and falls back to abs(a*e) for projected CRSes. - floodpath/routing/accumulate.py: accumulate_runoff(runoff, flow_grid) -> AccumulatedRunoffGrid. Reprojects Q via Resampling.average (correct for area-weighted depth aggregation), zero-fills NaN before accuflux so the SCS-CN nodata patch doesn't poison downstream cells. - floodpath/routing/discharge.py: peak_discharge(accumulated, duration_s) -> DischargeGrid with Q_peak = V_acc / T per cell. - AccumulatedRunoffGrid carries total_volume_m3() (=outlet, NOT sum of cells which would double-count); DischargeGrid carries outlet_peak(). - Tests: 30 new — cell_areas_m2 sanity at equator vs 60 deg N (verifies cos(lat) shrinkage), peak-discharge V/T arithmetic with NRCS-default 6-hour duration, robust pinned values against the existing Robit Bata fixtures: outlet V_acc = 1.43 Mm^3, peak Q at outlet = 66.2 m^3/s for 100 mm uniform rain over 6-hour event. Smoke test gains stage 15 (two-panel: log-scale V_acc + Q_peak with stream overlay) and a new --event-duration-h CLI flag (default 6.0). --- README.md | 1 + floodpath/routing/__init__.py | 32 ++++++ floodpath/routing/accumulate.py | 101 +++++++++++++++++++ floodpath/routing/constants.py | 23 +++++ floodpath/routing/discharge.py | 51 ++++++++++ floodpath/routing/models.py | 98 +++++++++++++++++++ floodpath/routing/utils.py | 60 ++++++++++++ tests/routing/__init__.py | 0 tests/routing/test_accumulate.py | 162 +++++++++++++++++++++++++++++++ tests/routing/test_discharge.py | 78 +++++++++++++++ tests/routing/test_models.py | 72 ++++++++++++++ tests/routing/test_utils.py | 79 +++++++++++++++ 12 files changed, 757 insertions(+) create mode 100644 floodpath/routing/__init__.py create mode 100644 floodpath/routing/accumulate.py create mode 100644 floodpath/routing/constants.py create mode 100644 floodpath/routing/discharge.py create mode 100644 floodpath/routing/models.py create mode 100644 floodpath/routing/utils.py create mode 100644 tests/routing/__init__.py create mode 100644 tests/routing/test_accumulate.py create mode 100644 tests/routing/test_discharge.py create mode 100644 tests/routing/test_models.py create mode 100644 tests/routing/test_utils.py diff --git a/README.md b/README.md index 4e1b328..bd62e04 100644 --- a/README.md +++ b/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 diff --git a/floodpath/routing/__init__.py b/floodpath/routing/__init__.py new file mode 100644 index 0000000..a2a3915 --- /dev/null +++ b/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", +] diff --git a/floodpath/routing/accumulate.py b/floodpath/routing/accumulate.py new file mode 100644 index 0000000..579f313 --- /dev/null +++ b/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, + ) diff --git a/floodpath/routing/constants.py b/floodpath/routing/constants.py new file mode 100644 index 0000000..03d75b4 --- /dev/null +++ b/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" diff --git a/floodpath/routing/discharge.py b/floodpath/routing/discharge.py new file mode 100644 index 0000000..a793c19 --- /dev/null +++ b/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, + ) diff --git a/floodpath/routing/models.py b/floodpath/routing/models.py new file mode 100644 index 0000000..4599975 --- /dev/null +++ b/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)) diff --git a/floodpath/routing/utils.py b/floodpath/routing/utils.py new file mode 100644 index 0000000..bed261e --- /dev/null +++ b/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() diff --git a/tests/routing/__init__.py b/tests/routing/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/tests/routing/test_accumulate.py b/tests/routing/test_accumulate.py new file mode 100644 index 0000000..bf8236b --- /dev/null +++ b/tests/routing/test_accumulate.py @@ -0,0 +1,162 @@ +"""Unit + fixture tests for routing.accumulate.accumulate_runoff.""" + +import numpy as np +import pytest +from rasterio.transform import Affine + +from floodpath.dem.models import BoundingBox +from floodpath.hydrology.models import FlowGrid +from floodpath.landuse.models import LanduseGrid +from floodpath.precip import uniform_precip_like +from floodpath.routing import accumulate_runoff +from floodpath.runoff import apply_scs_cn, compute_curve_number +from floodpath.runoff.models import CurveNumberGrid, RunoffGrid +from floodpath.soil import texture_to_hsg +from floodpath.soil.models import TextureGrid + + +class TestAccumulateRunoffMetadata: + def test_inherits_flow_grid_geom( + self, + robit_bata_flow_grid: FlowGrid, + robit_bata_worldcover: LanduseGrid, + robit_bata_soil: TextureGrid, + ) -> None: + hsg = texture_to_hsg(robit_bata_soil) + cn = compute_curve_number(robit_bata_worldcover, hsg) + precip = uniform_precip_like(cn, depth_mm=100.0) + runoff = apply_scs_cn(cn, precip) + acc = accumulate_runoff(runoff, robit_bata_flow_grid) + assert acc.shape == robit_bata_flow_grid.shape + assert acc.transform == robit_bata_flow_grid.transform + assert acc.crs == robit_bata_flow_grid.crs + assert acc.precip_source == "uniform" + assert acc.method == "steady-state-flow-accumulation" + + +class TestRobitBataAccumulation: + """Pinned values from running the full Robit Bata pipeline: + DEM -> flow grid + landuse + HSG -> CN + uniform 100 mm -> SCS-CN runoff + runoff -> accumulate_runoff -> volume per cell + """ + + @pytest.fixture + def robit_bata_acc_uniform_100mm( + self, + robit_bata_flow_grid: FlowGrid, + robit_bata_worldcover: LanduseGrid, + robit_bata_soil: TextureGrid, + ) -> "AccumulatedRunoffGrid": + from floodpath.routing.models import AccumulatedRunoffGrid # noqa + hsg = texture_to_hsg(robit_bata_soil) + cn = compute_curve_number(robit_bata_worldcover, hsg) + precip = uniform_precip_like(cn, depth_mm=100.0) + runoff = apply_scs_cn(cn, precip) + return accumulate_runoff(runoff, robit_bata_flow_grid) + + def test_outlet_volume_pinned(self, robit_bata_acc_uniform_100mm) -> None: + # Pinned: maximum accumulated volume in the patch (the basin + # outlet) is 1.43 Mm^3 for 100 mm uniform rain over the + # Vertisol-rich (mean Q = 67.5 mm), 22,705-cell catchment. + outlet = robit_bata_acc_uniform_100mm.total_volume_m3() + assert outlet == pytest.approx(1_429_000.0, rel=0.01) + + def test_min_is_nonnegative(self, robit_bata_acc_uniform_100mm) -> None: + # Accumulated runoff volume must never be negative. + assert robit_bata_acc_uniform_100mm.values.min() >= 0.0 + + def test_outlet_is_at_basin_outlet_pixel( + self, + robit_bata_acc_uniform_100mm, + robit_bata_flow_grid: FlowGrid, + ) -> None: + # The cell with the highest *runoff* accumulation must coincide + # with the cell of highest *flow* accumulation — they're driven + # by the same flow direction grid. + idx_runoff = np.unravel_index( + np.argmax(robit_bata_acc_uniform_100mm.values), + robit_bata_acc_uniform_100mm.shape, + ) + idx_flow = np.unravel_index( + np.argmax(robit_bata_flow_grid.accumulation), + robit_bata_flow_grid.shape, + ) + assert idx_runoff == idx_flow + + def test_outlet_bounded_by_patch_runoff_total( + self, + robit_bata_acc_uniform_100mm, + robit_bata_flow_grid: FlowGrid, + robit_bata_worldcover: LanduseGrid, + robit_bata_soil: TextureGrid, + ) -> None: + # The outlet's accumulated volume drains its upstream basin only, + # which is a SUBSET of the bbox (most cells drain off the edges + # of a square crop). So the outlet volume must be > 0 and + # <= total per-cell runoff volume across the whole patch. + from floodpath.routing.utils import cell_areas_m2 + import rasterio.warp + from rasterio.enums import Resampling + + hsg = texture_to_hsg(robit_bata_soil) + cn = compute_curve_number(robit_bata_worldcover, hsg) + precip = uniform_precip_like(cn, depth_mm=100.0) + runoff = apply_scs_cn(cn, precip) + + q_at_flow = np.zeros(robit_bata_flow_grid.shape, dtype=np.float32) + rasterio.warp.reproject( + source=runoff.values, + destination=q_at_flow, + src_transform=runoff.transform, + src_crs=runoff.crs, + dst_transform=robit_bata_flow_grid.transform, + dst_crs=robit_bata_flow_grid.crs, + resampling=Resampling.average, + ) + q_at_flow = np.where(np.isnan(q_at_flow), 0.0, q_at_flow) + areas = cell_areas_m2( + robit_bata_flow_grid.transform, + robit_bata_flow_grid.shape, + robit_bata_flow_grid.crs, + ) + total_volume = float(np.sum(q_at_flow * areas / 1000.0)) + + outlet = robit_bata_acc_uniform_100mm.total_volume_m3() + assert outlet > 0 + assert outlet <= total_volume + + def test_outlet_basin_share_pinned( + self, + robit_bata_acc_uniform_100mm, + robit_bata_flow_grid: FlowGrid, + ) -> None: + # Pinned: the outlet's basin is 22,705 cells out of 270*270=72,900 + # patch cells. Equivalently, the outlet's accumulated volume + # divided by the per-cell mean volume should equal 22,705 cells. + # This catches off-by-one / wrong-direction accuflux bugs. + outlet_volume = robit_bata_acc_uniform_100mm.total_volume_m3() + # Mean per-cell volume in the accumulated grid IS the per-cell + # contribution (since accuflux includes the cell's own volume). + # Use the smaller-than-mean cells (pure ridge cells with one + # contribution) as the per-cell volume proxy. + per_cell_volume = float(np.median(robit_bata_acc_uniform_100mm.values[ + robit_bata_acc_uniform_100mm.values > 0 + ])) + # Outlet basin must be at least 1% of the patch. + assert outlet_volume / per_cell_volume > 100 + + def test_zero_precip_yields_zero_accumulation( + self, + robit_bata_flow_grid: FlowGrid, + robit_bata_worldcover: LanduseGrid, + robit_bata_soil: TextureGrid, + ) -> None: + # No rain -> no runoff -> accumulation is zero everywhere. + hsg = texture_to_hsg(robit_bata_soil) + cn = compute_curve_number(robit_bata_worldcover, hsg) + precip = uniform_precip_like(cn, depth_mm=0.0) + runoff = apply_scs_cn(cn, precip) + acc = accumulate_runoff(runoff, robit_bata_flow_grid) + assert acc.values.max() == pytest.approx(0.0, abs=1e-3) diff --git a/tests/routing/test_discharge.py b/tests/routing/test_discharge.py new file mode 100644 index 0000000..04db573 --- /dev/null +++ b/tests/routing/test_discharge.py @@ -0,0 +1,78 @@ +"""Unit tests for routing.discharge.peak_discharge.""" + +import numpy as np +import pytest +from rasterio.transform import Affine + +from floodpath.dem.models import BoundingBox +from floodpath.routing import ( + AccumulatedRunoffGrid, + DEFAULT_EVENT_DURATION_S, + peak_discharge, +) + + +def _toy_acc(values: np.ndarray) -> AccumulatedRunoffGrid: + return AccumulatedRunoffGrid( + 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]), + precip_source="uniform", + method="steady-state-flow-accumulation", + ) + + +class TestPeakDischarge: + def test_basic_division(self) -> None: + # 21,600 m^3 over 6 hours = 1.0 m^3/s. + acc = _toy_acc(np.array([[21600.0]])) + q = peak_discharge(acc, duration_s=21600.0) + assert q.values[0, 0] == pytest.approx(1.0, rel=1e-5) + + def test_default_duration_helper(self) -> None: + # Pass DEFAULT_EVENT_DURATION_S explicitly. + acc = _toy_acc(np.array([[DEFAULT_EVENT_DURATION_S]])) + q = peak_discharge(acc, duration_s=DEFAULT_EVENT_DURATION_S) + assert q.values[0, 0] == pytest.approx(1.0, rel=1e-5) + + def test_negative_duration_raises(self) -> None: + acc = _toy_acc(np.array([[100.0]])) + with pytest.raises(ValueError, match="must be > 0"): + peak_discharge(acc, duration_s=-1.0) + + def test_zero_duration_raises(self) -> None: + acc = _toy_acc(np.array([[100.0]])) + with pytest.raises(ValueError, match="must be > 0"): + peak_discharge(acc, duration_s=0.0) + + def test_metadata_propagated(self) -> None: + acc = _toy_acc(np.array([[100.0]])) + q = peak_discharge(acc, duration_s=3600.0) + assert q.duration_s == 3600.0 + assert q.precip_source == acc.precip_source + assert q.method == acc.method + assert q.transform == acc.transform + assert q.crs == acc.crs + assert q.bbox == acc.bbox + + def test_shape_inherited(self) -> None: + acc = _toy_acc(np.full((4, 5), 100.0)) + q = peak_discharge(acc, duration_s=3600.0) + assert q.shape == (4, 5) + + def test_outlet_peak_matches_v_over_t(self) -> None: + # Outlet (max accumulation) divided by T must equal the outlet + # discharge from the DischargeGrid helper. + values = np.array([[10.0, 50.0, 30.0]]) + acc = _toy_acc(values) + T = 3600.0 + q = peak_discharge(acc, duration_s=T) + assert q.outlet_peak() == pytest.approx(50.0 / T) + + def test_short_duration_gives_higher_peak(self) -> None: + # Halving the event duration doubles the peak discharge. + acc = _toy_acc(np.array([[7200.0]])) # 7200 m^3 accumulated + q1 = peak_discharge(acc, duration_s=7200.0) # 1 m^3/s + q2 = peak_discharge(acc, duration_s=3600.0) # 2 m^3/s + assert q2.values[0, 0] == pytest.approx(2.0 * q1.values[0, 0]) diff --git a/tests/routing/test_models.py b/tests/routing/test_models.py new file mode 100644 index 0000000..b39eb6a --- /dev/null +++ b/tests/routing/test_models.py @@ -0,0 +1,72 @@ +"""Unit tests for routing.models (AccumulatedRunoffGrid, DischargeGrid).""" + +import numpy as np +import pytest +from rasterio.transform import Affine + +from floodpath.dem.models import BoundingBox +from floodpath.routing import AccumulatedRunoffGrid, DischargeGrid + + +def _toy_acc(values: np.ndarray) -> AccumulatedRunoffGrid: + return AccumulatedRunoffGrid( + 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]), + precip_source="uniform", + method="steady-state-flow-accumulation", + ) + + +def _toy_disch(values: np.ndarray) -> DischargeGrid: + return DischargeGrid( + 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]), + duration_s=21600.0, + precip_source="uniform", + method="steady-state-flow-accumulation", + ) + + +class TestAccumulatedRunoffGrid: + def test_shape(self) -> None: + a = _toy_acc(np.array([[1.0, 2.0, 3.0]])) + assert a.shape == (1, 3) + + def test_total_volume_is_max_value(self) -> None: + # `total_volume_m3()` returns the largest accumulated value (the + # outlet) — NOT the sum, which would double-count cells. + a = _toy_acc(np.array([[1.0, 5.0, 2.0]])) + assert a.total_volume_m3() == pytest.approx(5.0) + + def test_stats(self) -> None: + a = _toy_acc(np.array([[1.0, 2.0, 3.0]])) + s = a.stats() + assert s["min"] == 1.0 + assert s["max"] == 3.0 + assert s["mean"] == pytest.approx(2.0) + + +class TestDischargeGrid: + def test_shape(self) -> None: + d = _toy_disch(np.array([[0.5, 1.0, 1.5]])) + assert d.shape == (1, 3) + + def test_outlet_peak(self) -> None: + d = _toy_disch(np.array([[0.5, 10.0, 1.5]])) + assert d.outlet_peak() == pytest.approx(10.0) + + def test_outlet_peak_skips_nan(self) -> None: + d = _toy_disch(np.array([[np.nan, 1.0, 2.0]])) + assert d.outlet_peak() == pytest.approx(2.0) + + def test_units_are_m3_per_s(self) -> None: + d = _toy_disch(np.array([[1.0]])) + assert "m^3/s" in d.units + + def test_duration_s_preserved(self) -> None: + d = _toy_disch(np.array([[1.0]])) + assert d.duration_s == 21600.0 diff --git a/tests/routing/test_utils.py b/tests/routing/test_utils.py new file mode 100644 index 0000000..237d496 --- /dev/null +++ b/tests/routing/test_utils.py @@ -0,0 +1,79 @@ +"""Unit tests for routing.utils.cell_areas_m2.""" + +import numpy as np +import pytest +from rasterio.transform import Affine + +from floodpath.routing import cell_areas_m2 +from floodpath.routing.constants import DEG_LAT_M + + +class TestProjectedCrs: + def test_uniform_metres(self) -> None: + # A 30 m × 30 m projected pixel: every cell area must be 900 m^2. + transform = Affine.translation(0, 0) * Affine.scale(30.0, -30.0) + areas = cell_areas_m2(transform, shape=(4, 4), crs="EPSG:32636") + assert areas.shape == (4, 4) + np.testing.assert_array_equal(areas, np.full((4, 4), 900.0)) + + def test_non_square_pixel(self) -> None: + # 25 × 50 m pixels -> 1250 m^2 each. + transform = Affine.translation(0, 0) * Affine.scale(25.0, -50.0) + areas = cell_areas_m2(transform, shape=(2, 2), crs="EPSG:32636") + np.testing.assert_array_equal(areas, np.full((2, 2), 1250.0)) + + +class TestGeographicCrs: + def test_equator_pixel_area(self) -> None: + # 1 deg × 1 deg pixel at the equator: cos(0) = 1, area = (DEG_LAT_M)^2. + # transform top-edge latitude = 1 deg, e = -1 deg -> centre row at 0.5 deg. + transform = Affine(1.0, 0.0, 0.0, 0.0, -1.0, 1.0) # row 0 centre at lat 0.5 + areas = cell_areas_m2(transform, shape=(1, 1), crs="EPSG:4326") + # cos(0.5 deg) ~ 0.99996; effectively (DEG_LAT_M)^2. + assert areas.shape == (1, 1) + expected = DEG_LAT_M * DEG_LAT_M * np.cos(np.radians(0.5)) + assert areas[0, 0] == pytest.approx(expected, rel=1e-6) + + def test_high_latitude_pixel_smaller(self) -> None: + # At 60 deg N, cos(60) = 0.5 -> pixel area is HALF the equator's. + # Build two grids: one at equator, one at 60 deg N. + eq = Affine(0.01, 0.0, 0.0, 0.0, -0.01, 0.005) # row 0 centre at lat 0 + hi = Affine(0.01, 0.0, 0.0, 0.0, -0.01, 60.005) # row 0 centre at lat 60 + a_eq = cell_areas_m2(eq, shape=(1, 1), crs="EPSG:4326") + a_hi = cell_areas_m2(hi, shape=(1, 1), crs="EPSG:4326") + assert a_hi[0, 0] == pytest.approx(a_eq[0, 0] * 0.5, rel=1e-3) + + def test_southern_hemisphere(self) -> None: + # cos(-30) == cos(30); area at -30 deg should equal area at +30 deg. + north = Affine(0.01, 0.0, 0.0, 0.0, -0.01, 30.005) + south = Affine(0.01, 0.0, 0.0, 0.0, -0.01, -29.995) # row centre = -30 + a_n = cell_areas_m2(north, shape=(1, 1), crs="EPSG:4326") + a_s = cell_areas_m2(south, shape=(1, 1), crs="EPSG:4326") + assert a_n[0, 0] == pytest.approx(a_s[0, 0], rel=1e-6) + + def test_robit_bata_dem_cell_area(self) -> None: + # Sanity: Robit Bata DEM at lat ~11.8, 1 arc-sec pixels. + # Expected area ≈ 30 m × 30 m × cos(11.8°) ≈ 940 m^2. + # transform.a = e = 1/3600 deg. + pixel = 1.0 / 3600.0 + transform = Affine(pixel, 0.0, 37.5625, 0.0, -pixel, 11.8425) + areas = cell_areas_m2(transform, shape=(270, 270), crs="EPSG:4326") + # Centre row of the patch is at lat ~ 11.805. + centre_area = float(areas[135, 0]) + # Expected: pixel * 111319.5 * cos(11.8) = 30.21 m (x) + # pixel * 111319.5 = 30.92 m (y) + # Area ≈ 934 m^2. + assert centre_area == pytest.approx(934.0, abs=2.0) + + +class TestRowVariation: + def test_areas_decrease_northward(self) -> None: + # In the northern hemisphere, area decreases as latitude increases — + # so within a single grid, the *northern* (smaller-row-index) rows + # have *smaller* areas than the southern rows. + # transform.f = top-edge lat = 60, transform.e = -1 -> rows go 60, 59, 58, ... + transform = Affine(1.0, 0.0, 0.0, 0.0, -1.0, 60.0) + areas = cell_areas_m2(transform, shape=(3, 3), crs="EPSG:4326") + # Row 0 centre = 59.5 deg, row 2 centre = 57.5 deg. + # cos(59.5) < cos(57.5), so row 0 area should be smaller. + assert areas[0, 0] < areas[2, 0]