From 6f515c400b5e94e444ea554c9de6eea0f91a366b Mon Sep 17 00:00:00 2001 From: rehsani Date: Mon, 4 May 2026 13:56:36 -0700 Subject: [PATCH] Extend compute_damage to accept rainfall-driven inundation MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Damage now consumes either kind of inundation depth — the static InundationDepth (user-supplied scalar water level) or the new RainfallInundationDepth (per-cell depth field driven by Manning's normal-depth + HAND broadcast). Same numerics either way; the DamageMap's metadata records which scenario produced it via a new `scenario` string field plus the existing `water_level`, which now holds the per-cell maximum depth as a one-number summary for the rainfall path. - floodpath/damage/models.py: DamageMap gains `scenario: str = "static"`, defaulted so the existing 35 damage tests pass unchanged. - floodpath/damage/compute.py: `depth` argument typed as `InundationDepth | RainfallInundationDepth`. Inside, branch on type to set scenario string ("static water level X m" vs "rainfall-driven (source, T, method)") and water_level (input scalar vs per-cell max). Tests: 9 new for the rainfall path: - scenario string content (precip source, duration, method) - water_level == max(values) - numerical equivalence between rainfall and static paths at same depth - upsample-preserving-total still works with rainfall depth + coarser exposure - end-to-end pinned: 100 mm × 6 hr on Robit Bata fixtures yields 7,355 m^2 damaged built-up — much smaller than the hypothetical static 5 m scenario's 69,199 m^2 (consistent with rainfall flood being a narrow channel ribbon vs static water level flooding broad valleys) Smoke test gains stage 17: applies JRC residential to the stage-16 rainfall flood, prints both totals side-by-side for direct comparison with the existing stage 8 static damage. --- README.md | 2 +- floodpath/damage/compute.py | 36 ++++- floodpath/damage/models.py | 11 +- tests/damage/test_rainfall_damage.py | 196 +++++++++++++++++++++++++++ 4 files changed, 239 insertions(+), 6 deletions(-) create mode 100644 tests/damage/test_rainfall_damage.py diff --git a/README.md b/README.md index 6ddaa7b..0ef1186 100644 --- a/README.md +++ b/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 diff --git a/floodpath/damage/compute.py b/floodpath/damage/compute.py index af886aa..03581e9 100644 --- a/floodpath/damage/compute.py +++ b/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, ) diff --git a/floodpath/damage/models.py b/floodpath/damage/models.py index 1abbb50..21195b0 100644 --- a/floodpath/damage/models.py +++ b/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]: diff --git a/tests/damage/test_rainfall_damage.py b/tests/damage/test_rainfall_damage.py new file mode 100644 index 0000000..02ccdde --- /dev/null +++ b/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)