Extend compute_damage to accept rainfall-driven inundation

README.md

@@ -76,7 +76,7 @@ print(f"Total damaged built-up: {damage.values.sum():,.0f} m²")
 | `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) + 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 |
+| `floodpath.damage` | JRC Huizinga 2017 + DEM/HAND/GHSL/routing | Per-cell flood depth and damage in m² of built-up surface — accepts either a static water-level scenario or a rainfall-driven flood from `floodpath.routing` |
 
 ## Depth-damage curves
 

floodpath/damage/compute.py

@@ -1,13 +1,18 @@
 """Combine inundation depth + exposure + depth-damage curve -> DamageMap."""
 
+from typing import Union
+
 import numpy as np
 
 from floodpath.exposure.models import ExposureGrid
+from floodpath.routing.flood import RainfallInundationDepth
 
 from .constants import RATIO_TOLERANCE
 from .curves import DepthDamageCurve
 from .models import DamageMap, InundationDepth
 
+InundationLike = Union[InundationDepth, RainfallInundationDepth]
+
 
 def _upsample_preserving_total(
     values: np.ndarray,
@@ -26,7 +31,7 @@ def _upsample_preserving_total(
 
 
 def compute_damage(
-    depth: InundationDepth,
+    depth: InundationLike,
     exposure: ExposureGrid,
     curve: DepthDamageCurve,
 ) -> DamageMap:
@@ -38,8 +43,19 @@ def compute_damage(
     built-up area is preserved. Per-cell damage =
     depth_damage_fraction(depth) * upsampled_built_up_area.
 
+    Accepts either:
+      * `InundationDepth` from `compute_inundation_depth(hand, water_level=L)`
+        (static scenario, single user-supplied L); produces
+        `scenario = "static water level L m"`.
+      * `RainfallInundationDepth` from `compute_rainfall_inundation(...)`
+        (rainfall-driven, per-cell water depth varies with the upstream
+        Manning solution); produces a scenario string capturing the
+        precip source, event duration, and hydraulic method, with
+        `water_level` set to the per-cell maximum depth as a scalar
+        summary.
+
     Args:
-        depth: InundationDepth from compute_inundation_depth.
+        depth: Either kind of inundation depth grid.
         exposure: ExposureGrid carrying the per-cell built-up surface (m^2).
         curve: Depth-damage curve to apply (e.g. JRC_AFRICA_RESIDENTIAL).
 
@@ -88,12 +104,26 @@ def compute_damage(
     damage_fraction = curve(depth.values)
     damage_values = damage_fraction * upsampled_exposure
 
+    # Scenario metadata: differs between static and rainfall-driven inputs.
+    if isinstance(depth, RainfallInundationDepth):
+        # Use the per-cell max depth as a one-number summary so DamageMaps
+        # remain comparable via `water_level`.
+        scenario_water_level = float(depth.values.max())
+        scenario = (
+            f"rainfall-driven ({depth.discharge_source}, "
+            f"T={depth.duration_s / 3600:.1f} hr, {depth.method})"
+        )
+    else:
+        scenario_water_level = depth.water_level
+        scenario = f"static water level {depth.water_level:.2f} m"
+
     return DamageMap(
         values=damage_values,
         transform=depth.transform,
         crs=depth.crs,
         bbox=depth.bbox,
-        water_level=depth.water_level,
+        water_level=scenario_water_level,
         curve_name=curve.name,
         units=f"damaged {exposure.units}",
+        scenario=scenario,
     )

floodpath/damage/models.py

@@ -34,9 +34,15 @@ class DamageMap:
     """Per-cell flood damage raster with metadata about how it was produced.
 
     `values` carries the per-cell damage in the unit reported by `units`
-    (typically m^2 of damaged built-up surface). `water_level` and
-    `curve_name` document the scenario assumption so two DamageMaps can
+    (typically m^2 of damaged built-up surface). `water_level`, `scenario`,
+    and `curve_name` document the input assumptions so two DamageMaps can
     be compared without ambiguity.
+
+    For static-flood scenarios `water_level` is the scalar input that drove
+    the inundation; for rainfall-driven scenarios it's the maximum
+    per-cell depth observed (a useful one-number summary even when the
+    field varies spatially), and `scenario` records the precip source,
+    event duration, and hydraulic method.
     """
 
     values: np.ndarray
@@ -46,6 +52,7 @@ class DamageMap:
     water_level: float
     curve_name: str
     units: str
+    scenario: str = "static"
 
     @property
     def shape(self) -> tuple[int, int]:

tests/damage/test_rainfall_damage.py

@@ -0,0 +1,196 @@
+"""Tests for compute_damage on the rainfall-driven inundation path.
+
+The static-water-level damage tests live in test_compute.py and continue
+to pass — this module only exercises the new RainfallInundationDepth
+branch and the scenario metadata.
+"""
+
+import numpy as np
+import pytest
+from rasterio.transform import Affine
+
+from floodpath.dem.models import BoundingBox
+from floodpath.exposure.models import ExposureGrid
+from floodpath.damage import (
+    JRC_AFRICA_RESIDENTIAL,
+    compute_damage,
+)
+from floodpath.routing.flood import RainfallInundationDepth
+
+
+def _toy_rainfall_inundation(
+    values: np.ndarray,
+    duration_s: float = 21600.0,
+    discharge_source: str = "uniform",
+) -> RainfallInundationDepth:
+    return RainfallInundationDepth(
+        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=duration_s,
+        discharge_source=discharge_source,
+    )
+
+
+def _toy_exposure(values: np.ndarray) -> ExposureGrid:
+    return ExposureGrid(
+        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]),
+        epoch=2020,
+        units="m^2 of built-up surface per cell",
+    )
+
+
+class TestRainfallDamageScenarioMetadata:
+    def test_scenario_string_includes_source_and_duration(self) -> None:
+        # Hand-pick depths so JRC_AFRICA_RESIDENTIAL evaluates at known
+        # points without ambiguity.
+        depth = _toy_rainfall_inundation(
+            np.array([[0.0, 1.0, 2.0]]),
+            duration_s=3600.0,
+            discharge_source="ERA5",
+        )
+        exposure = _toy_exposure(np.full((1, 3), 100.0))
+        damage = compute_damage(depth, exposure, JRC_AFRICA_RESIDENTIAL)
+        assert "rainfall-driven" in damage.scenario
+        assert "ERA5" in damage.scenario
+        assert "T=1.0 hr" in damage.scenario
+        assert "manning" in damage.scenario.lower()
+
+    def test_water_level_set_to_max_depth(self) -> None:
+        # For rainfall scenarios, water_level should be the per-cell max
+        # depth (a one-number summary), not zero or the input duration.
+        depth = _toy_rainfall_inundation(np.array([[0.5, 3.7, 1.0]]))
+        exposure = _toy_exposure(np.full((1, 3), 100.0))
+        damage = compute_damage(depth, exposure, JRC_AFRICA_RESIDENTIAL)
+        assert damage.water_level == pytest.approx(3.7, abs=1e-5)
+
+    def test_curve_name_propagated(self) -> None:
+        depth = _toy_rainfall_inundation(np.array([[1.0]]))
+        exposure = _toy_exposure(np.full((1, 1), 100.0))
+        damage = compute_damage(depth, exposure, JRC_AFRICA_RESIDENTIAL)
+        assert damage.curve_name == JRC_AFRICA_RESIDENTIAL.name
+
+
+class TestRainfallDamageNumerics:
+    def test_zero_depth_yields_zero_damage(self) -> None:
+        depth = _toy_rainfall_inundation(np.zeros((3, 3)))
+        exposure = _toy_exposure(np.full((3, 3), 200.0))
+        damage = compute_damage(depth, exposure, JRC_AFRICA_RESIDENTIAL)
+        assert damage.values.sum() == pytest.approx(0.0)
+
+    def test_aligned_grids_no_upsample(self) -> None:
+        # depth and exposure share the same grid -> no upsampling needed.
+        # JRC_AFRICA_RESIDENTIAL at 1.0 m depth on a 100 m^2 cell: damage =
+        # fraction(1m) * 100. This must equal what the static path returns
+        # for the same depth.
+        depth = _toy_rainfall_inundation(np.array([[1.0]]))
+        exposure = _toy_exposure(np.full((1, 1), 100.0))
+        damage = compute_damage(depth, exposure, JRC_AFRICA_RESIDENTIAL)
+        expected_fraction = float(JRC_AFRICA_RESIDENTIAL(np.array([1.0]))[0])
+        assert damage.values[0, 0] == pytest.approx(expected_fraction * 100.0)
+
+    def test_upsampled_grids_preserve_total_exposure(self) -> None:
+        # depth grid 3x finer than exposure (factor=3 upsampling), but
+        # both must cover the same geographic bbox. Exposure: 1x1 cell
+        # of 900 m^2; depth: 3x3 cells, each at the 1.0 m depth.
+        # After upsampling, each fine cell holds 900/9 = 100 m^2.
+        # Total damage = fraction(1m) * 900 m^2.
+        bbox = BoundingBox(0.0, 0.0, 0.003, 0.003)
+        # Depth: 3x3 with pixel 0.001
+        depth = RainfallInundationDepth(
+            values=np.full((3, 3), 1.0, dtype=np.float32),
+            transform=Affine.translation(0, 0) * Affine.scale(0.001, -0.001),
+            crs="EPSG:4326",
+            bbox=bbox,
+            duration_s=21600.0,
+            discharge_source="uniform",
+        )
+        # Exposure: 1x1 spanning the same bbox (pixel = 0.003).
+        exposure = ExposureGrid(
+            values=np.array([[900.0]], dtype=np.float32),
+            transform=Affine.translation(0, 0) * Affine.scale(0.003, -0.003),
+            crs="EPSG:4326",
+            bbox=bbox,
+            epoch=2020,
+            units="m^2 of built-up surface per cell",
+        )
+        damage = compute_damage(depth, exposure, JRC_AFRICA_RESIDENTIAL)
+        # Each of 9 fine cells: fraction * (900/9). Total damage = fraction * 900.
+        f = float(JRC_AFRICA_RESIDENTIAL(np.array([1.0]))[0])
+        assert damage.values.sum() == pytest.approx(f * 900.0, rel=1e-5)
+
+
+class TestRobitBataRainfallDamage:
+    """End-to-end pinned values on the Robit Bata fixtures."""
+
+    @pytest.fixture
+    def rainfall_damage(
+        self,
+        robit_bata_dem,
+        robit_bata_flow_grid,
+        robit_bata_streams,
+        robit_bata_hand,
+        robit_bata_worldcover,
+        robit_bata_soil,
+        robit_bata_ghsl,
+    ):
+        from floodpath.landuse import landuse_to_roughness
+        from floodpath.precip import uniform_precip_like
+        from floodpath.routing import (
+            accumulate_runoff, peak_discharge,
+            compute_water_level, compute_rainfall_inundation,
+        )
+        from floodpath.runoff import apply_scs_cn, compute_curve_number
+        from floodpath.soil import texture_to_hsg
+
+        hsg = texture_to_hsg(robit_bata_soil)
+        cn = compute_curve_number(robit_bata_worldcover, hsg)
+        roughness = landuse_to_roughness(robit_bata_worldcover)
+        precip = uniform_precip_like(cn, depth_mm=100.0)
+        runoff = apply_scs_cn(cn, precip)
+        acc = accumulate_runoff(runoff, robit_bata_flow_grid)
+        discharge = peak_discharge(acc, duration_s=6 * 3600.0)
+        water_level = compute_water_level(
+            discharge, roughness, robit_bata_flow_grid,
+            robit_bata_streams, robit_bata_dem,
+        )
+        flood = compute_rainfall_inundation(
+            water_level, robit_bata_hand,
+            robit_bata_flow_grid, robit_bata_streams,
+        )
+        return compute_damage(flood, robit_bata_ghsl, JRC_AFRICA_RESIDENTIAL)
+
+    def test_total_damage_pinned(self, rainfall_damage) -> None:
+        # Pinned: 100 mm × 6 hr storm yields ~7,355 m^2 of damaged
+        # built-up surface across the Robit Bata patch. Significantly
+        # smaller than the 69,199 m^2 from the hypothetical static
+        # 5 m water level (the rainfall flood is a narrow ribbon).
+        assert rainfall_damage.values.sum() == pytest.approx(7355.0, abs=50.0)
+
+    def test_smaller_than_static_5m_scenario(
+        self,
+        rainfall_damage,
+        robit_bata_hand,
+        robit_bata_ghsl,
+    ) -> None:
+        # Sanity: 100 mm rainfall-driven damage should be much less
+        # than the hypothetical 5 m static water level scenario, since
+        # 5 m everywhere floods broad valleys while 100 mm rain only
+        # gives ~10 m at the outlet decaying upstream.
+        from floodpath.damage import compute_inundation_depth
+        static_depth = compute_inundation_depth(robit_bata_hand, water_level=5.0)
+        static_damage = compute_damage(static_depth, robit_bata_ghsl, JRC_AFRICA_RESIDENTIAL)
+        assert rainfall_damage.values.sum() < static_damage.values.sum()
+        # Roughly an order of magnitude smaller — flood extent is ~10x narrower.
+        assert rainfall_damage.values.sum() < 0.2 * static_damage.values.sum()
+
+    def test_scenario_metadata_pinned(self, rainfall_damage) -> None:
+        assert "rainfall-driven" in rainfall_damage.scenario
+        assert "uniform" in rainfall_damage.scenario
+        assert "T=6.0 hr" in rainfall_damage.scenario
+        # water_level set to per-cell max depth from rainfall flood.
+        assert rainfall_damage.water_level == pytest.approx(10.07, abs=0.10)