diff --git a/.gitignore b/.gitignore index f35b67f..f4924fe 100644 --- a/.gitignore +++ b/.gitignore @@ -32,6 +32,9 @@ smoke_test.py # Per-machine Claude Code settings. .claude/ +# PyPI upload credentials. NEVER commit. +.pypirc + # Local virtual envs (the project uses a named conda env, not a local one, # but guard against accidental commits if a contributor sets up differently). .venv/ diff --git a/README.md b/README.md index 0891230..c57a883 100644 --- a/README.md +++ b/README.md @@ -71,6 +71,7 @@ print(f"Total damaged built-up: {damage.values.sum():,.0f} m²") | `floodpath.dem` | Copernicus GLO-30 (AWS Open Data, COG) | Elevation patch around any (lat, lon) | | `floodpath.hydrology` | derived from DEM via `pyflwdir` | Flow direction + accumulation, stream networks (with Strahler order), basin delineation, HAND | | `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) | | `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/landuse/__init__.py b/floodpath/landuse/__init__.py new file mode 100644 index 0000000..99e774c --- /dev/null +++ b/floodpath/landuse/__init__.py @@ -0,0 +1,19 @@ +"""Landuse module: categorical land-cover rasters answering 'what is the surface here'.""" + +from .constants import WORLDCOVER_CLASSES +from .models import LanduseGrid +from .roughness import ( + WORLDCOVER_MANNING_N, + RoughnessGrid, + landuse_to_roughness, +) +from .worldcover import get_worldcover_landuse + +__all__ = [ + "LanduseGrid", + "RoughnessGrid", + "WORLDCOVER_CLASSES", + "WORLDCOVER_MANNING_N", + "get_worldcover_landuse", + "landuse_to_roughness", +] diff --git a/floodpath/landuse/constants.py b/floodpath/landuse/constants.py new file mode 100644 index 0000000..db91186 --- /dev/null +++ b/floodpath/landuse/constants.py @@ -0,0 +1,46 @@ +"""Constants for the landuse module.""" + +from types import MappingProxyType + +# --- ESA WorldCover (10 m global land cover, AWS Open Data) ---------------- +# +# ESA publishes two annual maps as Cloud-Optimized GeoTIFFs on a public S3 +# bucket. Each year ships under its own product version directory. +# +# v100 -> 2020 release +# v200 -> 2021 release +# +# Globe is tiled at 3 deg x 3 deg; tile names encode the SW corner. + +WORLDCOVER_BASE_URL: str = ( + "https://esa-worldcover.s3.eu-central-1.amazonaws.com" +) +WORLDCOVER_VERSIONS: MappingProxyType = MappingProxyType({ + 2020: "v100", + 2021: "v200", +}) +WORLDCOVER_AVAILABLE_YEARS: tuple[int, ...] = tuple(WORLDCOVER_VERSIONS.keys()) +WORLDCOVER_DEFAULT_YEAR: int = 2021 + +WORLDCOVER_TILE_DEG: int = 3 +WORLDCOVER_BUFFER_DEG: float = 0.1 +WORLDCOVER_CRS: str = "EPSG:4326" +WORLDCOVER_NODATA: int = 0 +WORLDCOVER_UNITS: str = "ESA WorldCover class code" + +# Classes from the ESA WorldCover 2021 v200 product user manual. +# 2020 v100 uses the same set minus class 95 (Mangroves) at the original +# release; the per-cell codes themselves are stable. +WORLDCOVER_CLASSES: MappingProxyType = MappingProxyType({ + 10: "tree_cover", + 20: "shrubland", + 30: "grassland", + 40: "cropland", + 50: "built_up", + 60: "bare_or_sparse_vegetation", + 70: "snow_and_ice", + 80: "permanent_water_bodies", + 90: "herbaceous_wetland", + 95: "mangroves", + 100: "moss_and_lichen", +}) diff --git a/floodpath/landuse/models.py b/floodpath/landuse/models.py new file mode 100644 index 0000000..9989e0b --- /dev/null +++ b/floodpath/landuse/models.py @@ -0,0 +1,59 @@ +"""Dataclasses for the landuse module.""" + +from dataclasses import dataclass + +import numpy as np +from rasterio.transform import Affine + +from floodpath.dem.models import BoundingBox + +from .constants import WORLDCOVER_CLASSES + + +@dataclass +class LanduseGrid: + """A categorical land-cover raster clipped to a bbox. + + `values` carries per-cell class codes (uint8) drawn from a discrete + code-to-name mapping (e.g. WORLDCOVER_CLASSES). `epoch` is the year of + the source product. Use `class_counts()` to summarise composition; use + the `classes` mapping to translate codes to human-readable names. + """ + + values: np.ndarray + transform: Affine + crs: str + bbox: BoundingBox + epoch: int + classes: dict[int, str] + units: str + + @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 class_counts(self) -> dict[int, int]: + """Count cells per class code present in the grid. + + Returns: + Mapping from class code to cell count. Codes absent from the + grid are omitted. + """ + codes, counts = np.unique(self.values, return_counts=True) + return {int(c): int(n) for c, n in zip(codes, counts)} + + def class_name(self, code: int) -> str: + """Return the human-readable name for a class code, or 'unknown'.""" + return self.classes.get(code, "unknown") + + def fraction(self, code: int) -> float: + """Return the fraction of cells with the given class code.""" + if self.values.size == 0: + return 0.0 + return float(np.count_nonzero(self.values == code) / self.values.size) + + +# Re-export so callers can do `from floodpath.landuse import WORLDCOVER_CLASSES` +# while also seeing it through the dataclass `classes` field. +__all__ = ["LanduseGrid", "WORLDCOVER_CLASSES"] diff --git a/floodpath/landuse/roughness.py b/floodpath/landuse/roughness.py new file mode 100644 index 0000000..64f9a4c --- /dev/null +++ b/floodpath/landuse/roughness.py @@ -0,0 +1,112 @@ +"""Map ESA WorldCover land-cover classes to Manning's roughness coefficient n. + +This is a pure transformation of a categorical landuse raster into a +continuous roughness raster. Default values come from Chow 1959 Table 5-6 +('Open-channel hydraulics', floodplain values) cross-checked against +USACE EM 1110-2-1417 and FEMA flood-mapping guidance — i.e. n suitable +for floodplain/channel routing, not overland-flow runoff (the latter +takes much higher n; see Kalyanapu et al. 2009). + +Users with site-calibrated values should pass a custom mapping to +`landuse_to_roughness`. +""" + +from dataclasses import dataclass +from types import MappingProxyType + +import numpy as np +from rasterio.transform import Affine + +from floodpath.dem.models import BoundingBox + +from .models import LanduseGrid + +# Manning's n by ESA WorldCover class code (s/m^(1/3)). +# Channel/floodplain values from Chow 1959 Table 5-6, USACE EM 1110-2-1417. +WORLDCOVER_MANNING_N: MappingProxyType = MappingProxyType({ + 10: 0.150, # tree_cover — dense forest with undergrowth + 20: 0.080, # shrubland — scattered brush, heavy weeds + 30: 0.035, # grassland — short grass / pasture + 40: 0.040, # cropland — cultivated, mature row crops + 50: 0.060, # built_up — mixed urban + 60: 0.025, # bare_or_sparse_vegetation — bare earth + 70: 0.030, # snow_and_ice — smooth surface + 80: 0.025, # permanent_water_bodies — open water + 90: 0.080, # herbaceous_wetland — emergent wetland + 95: 0.100, # mangroves — woody wetland with roots + 100: 0.030, # moss_and_lichen +}) + +ROUGHNESS_UNITS: str = "Manning's n (s/m^(1/3))" + + +@dataclass +class RoughnessGrid: + """Per-cell Manning's roughness coefficient n. + + `values` carries the dimensionless n at each cell; `mapping` records + which class-to-n table was used so downstream consumers (e.g. a + Manning normal-depth solver) can introspect the assumptions. `fallback` + is the n value applied to landuse codes absent from `mapping`, + including the WorldCover nodata value (0). + """ + + values: np.ndarray + transform: Affine + crs: str + bbox: BoundingBox + mapping: dict[int, float] + fallback: float + units: str = ROUGHNESS_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 (min, max, mean, median) of n.""" + return { + "min": float(self.values.min()), + "max": float(self.values.max()), + "mean": float(self.values.mean()), + "median": float(np.median(self.values)), + } + + +def landuse_to_roughness( + landuse: LanduseGrid, + mapping: dict[int, float] | None = None, + fallback: float = 0.035, +) -> RoughnessGrid: + """Convert a categorical landuse raster to a Manning's roughness raster. + + Each cell's class code is looked up in `mapping`; codes absent from + the mapping (including WorldCover nodata=0) get `fallback`. + + Args: + landuse: Categorical land-cover grid (e.g. ESA WorldCover). + mapping: Class-code -> Manning's n. Defaults to WORLDCOVER_MANNING_N + (Chow 1959 / USACE-aligned floodplain values). + fallback: n applied to codes not present in `mapping`. Defaults to + 0.035 (cropland-like, a reasonable middle of the floodplain + range when nothing else is known). + + Returns: + RoughnessGrid sharing georef with the input landuse grid. + """ + table: dict[int, float] = ( + dict(WORLDCOVER_MANNING_N) if mapping is None else dict(mapping) + ) + n_values = np.full(landuse.values.shape, float(fallback), dtype=np.float32) + for code, n in table.items(): + n_values[landuse.values == code] = float(n) + + return RoughnessGrid( + values=n_values, + transform=landuse.transform, + crs=landuse.crs, + bbox=landuse.bbox, + mapping=table, + fallback=float(fallback), + ) diff --git a/floodpath/landuse/utils.py b/floodpath/landuse/utils.py new file mode 100644 index 0000000..84ecddd --- /dev/null +++ b/floodpath/landuse/utils.py @@ -0,0 +1,95 @@ +"""Pure helper functions for ESA WorldCover tile addressing.""" + +import math + +from floodpath.dem.models import BoundingBox + +from .constants import ( + WORLDCOVER_BASE_URL, + WORLDCOVER_TILE_DEG, + WORLDCOVER_VERSIONS, +) + + +def worldcover_tile_id( + lat: float, + lon: float, +) -> tuple[int, int]: + """Return the SW-corner (lat, lon) of the 3x3 deg WorldCover tile containing the point. + + ESA WorldCover tiles align to multiples of 3 degrees and are named by + their southwest corner. For example, (lat=11.8, lon=37.5) lies in the + tile rooted at (9, 36) -> "N09E036", covering [9, 12)N x [36, 39)E. + + Args: + lat: Latitude in decimal degrees. + lon: Longitude in decimal degrees. + + Returns: + Tuple `(lat_sw, lon_sw)` snapped down to the nearest 3-degree multiple. + """ + lat_sw = int(math.floor(lat / WORLDCOVER_TILE_DEG) * WORLDCOVER_TILE_DEG) + lon_sw = int(math.floor(lon / WORLDCOVER_TILE_DEG) * WORLDCOVER_TILE_DEG) + return lat_sw, lon_sw + + +def tiles_for_bbox(bbox: BoundingBox) -> list[tuple[int, int]]: + """Enumerate the WorldCover 3x3 deg tile SW-corners covering a bbox. + + Walks every interior tile between the bbox corners — a bbox spanning + more than two tiles in either dimension must include the interior + tiles too. + + Args: + bbox: Geographic bounding box in WGS84. + + Returns: + List of `(lat_sw, lon_sw)` tile corners snapped to 3-degree multiples. + """ + lat_sw_min, lon_sw_min = worldcover_tile_id(bbox.min_lat, bbox.min_lon) + lat_sw_max, lon_sw_max = worldcover_tile_id(bbox.max_lat, bbox.max_lon) + + # When a bbox edge sits exactly on a tile boundary, the max corner is the + # next tile up — exclude it. + if bbox.max_lat == lat_sw_max and lat_sw_max > lat_sw_min: + lat_sw_max -= WORLDCOVER_TILE_DEG + if bbox.max_lon == lon_sw_max and lon_sw_max > lon_sw_min: + lon_sw_max -= WORLDCOVER_TILE_DEG + + return [ + (lat, lon) + for lat in range(lat_sw_min, lat_sw_max + 1, WORLDCOVER_TILE_DEG) + for lon in range(lon_sw_min, lon_sw_max + 1, WORLDCOVER_TILE_DEG) + ] + + +def worldcover_tile_name( + lat_sw: int, + lon_sw: int, +) -> str: + """Build the ESA WorldCover tile name from a SW-corner. + + Tiles are named like 'N09E036' or 'S12W045': two-digit zero-padded + latitude, three-digit zero-padded longitude, with hemispheric letters. + """ + ns = "N" if lat_sw >= 0 else "S" + ew = "E" if lon_sw >= 0 else "W" + return f"{ns}{abs(lat_sw):02d}{ew}{abs(lon_sw):03d}" + + +def worldcover_tile_url( + lat_sw: int, + lon_sw: int, + year: int, +) -> str: + """Build the public HTTPS URL for an ESA WorldCover tile COG.""" + if year not in WORLDCOVER_VERSIONS: + raise ValueError( + f"year must be one of {tuple(WORLDCOVER_VERSIONS.keys())}, got {year}" + ) + version = WORLDCOVER_VERSIONS[year] + tile = worldcover_tile_name(lat_sw, lon_sw) + return ( + f"{WORLDCOVER_BASE_URL}/{version}/{year}/map/" + f"ESA_WorldCover_10m_{year}_{version}_{tile}_Map.tif" + ) diff --git a/floodpath/landuse/worldcover.py b/floodpath/landuse/worldcover.py new file mode 100644 index 0000000..e183fc3 --- /dev/null +++ b/floodpath/landuse/worldcover.py @@ -0,0 +1,91 @@ +"""Fetch ESA WorldCover 10 m land-cover patches by streaming COGs over HTTPS.""" + +import rasterio +from rasterio.env import Env +from rasterio.errors import RasterioIOError +from rasterio.merge import merge + +from floodpath.dem.constants import GDAL_VSICURL_ENV +from floodpath.dem.models import BoundingBox +from floodpath.dem.utils import bbox_from_point + +from .constants import ( + WORLDCOVER_AVAILABLE_YEARS, + WORLDCOVER_BUFFER_DEG, + WORLDCOVER_CLASSES, + WORLDCOVER_CRS, + WORLDCOVER_DEFAULT_YEAR, + WORLDCOVER_UNITS, +) +from .models import LanduseGrid +from .utils import tiles_for_bbox, worldcover_tile_url + + +def get_worldcover_landuse( + lat: float, + lon: float, + buffer_deg: float = WORLDCOVER_BUFFER_DEG, + year: int = WORLDCOVER_DEFAULT_YEAR, +) -> LanduseGrid: + """Return ESA WorldCover land-cover classes for the bbox around a point. + + ESA WorldCover ships as Cloud-Optimized GeoTIFFs on AWS Open Data + (eu-central-1). The bucket honours HTTP range requests, so this reads + only the bbox window via GDAL's /vsicurl/ driver — no full tile + download. Each tile is 3 deg x 3 deg at 10 m (~36k x 36k cells), so + streaming is meaningfully cheaper than caching. + + Args: + lat: Latitude of the region center, decimal degrees. + lon: Longitude of the region center, decimal degrees. + buffer_deg: Half-width of the square bbox around the point, degrees. + year: 2020 (v100) or 2021 (v200). Defaults to 2021. + + Returns: + LanduseGrid with per-cell class codes, georef, and the active + class-name mapping. + + Raises: + ValueError: `year` is not one of the available WorldCover releases. + FileNotFoundError: No WorldCover tile covers the requested bbox + (e.g. fully oceanic bboxes outside ESA's land mask). + """ + if year not in WORLDCOVER_AVAILABLE_YEARS: + raise ValueError( + f"year must be one of {WORLDCOVER_AVAILABLE_YEARS}, got {year}" + ) + + bbox = bbox_from_point(lat, lon, buffer_deg) + tile_corners = tiles_for_bbox(bbox) + urls = [worldcover_tile_url(lat_sw, lon_sw, year=year) for lat_sw, lon_sw in tile_corners] + + with Env(**GDAL_VSICURL_ENV): + datasets = [] + try: + for url in urls: + try: + datasets.append(rasterio.open(f"/vsicurl/{url}")) + except RasterioIOError: + # Fully oceanic tiles are absent from the bucket. + continue + + if not datasets: + raise FileNotFoundError( + f"No WorldCover {year} tiles found for bbox {bbox.bounds}. " + "The region may be ocean or outside ESA's land mask." + ) + + values, transform = merge(datasets, bounds=bbox.bounds) + finally: + for ds in datasets: + ds.close() + + return LanduseGrid( + values=values[0], + transform=transform, + crs=WORLDCOVER_CRS, + bbox=BoundingBox(*bbox.bounds), + epoch=year, + classes=dict(WORLDCOVER_CLASSES), + units=WORLDCOVER_UNITS, + ) diff --git a/tests/conftest.py b/tests/conftest.py index 7019248..0748bd6 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -14,6 +14,8 @@ from floodpath.exposure.models import BuildingCollection, ExposureGrid from floodpath.hydrology import build_flow_grid, compute_hand, extract_streams from floodpath.hydrology.models import FlowGrid, Hand, StreamNetwork +from floodpath.landuse import WORLDCOVER_CLASSES +from floodpath.landuse.models import LanduseGrid ROBIT_BATA_STREAM_THRESHOLD: int = 200 ROBIT_BATA_GHSL_FIXTURE: Path = ( @@ -25,6 +27,9 @@ ROBIT_BATA_OSM_FIXTURE: Path = ( Path(__file__).resolve().parent / "fixtures" / "robit_bata_osm.json" ) +ROBIT_BATA_WORLDCOVER_FIXTURE: Path = ( + Path(__file__).resolve().parent / "fixtures" / "robit_bata_worldcover.tif" +) ROBIT_BATA_FIXTURE: Path = Path(__file__).resolve().parent / "fixtures" / "robit_bata.tif" @@ -157,6 +162,35 @@ def robit_bata_osm() -> BuildingCollection: return parse_overpass_response(payload) +@pytest.fixture(scope="session") +def robit_bata_worldcover() -> LanduseGrid: + """Load the committed Robit Bata WorldCover fixture (900x900 patch, year 2021). + + Same bbox as the DEM/GHSL/WorldPop fixtures, sampled at ESA WorldCover + 10 m (v200, year 2021). Regenerate via + `python tests/fixtures/_generate_robit_bata_worldcover.py`. + """ + with rasterio.open(ROBIT_BATA_WORLDCOVER_FIXTURE) as src: + values = src.read(1) + transform = src.transform + crs = str(src.crs) + bounds = src.bounds + return LanduseGrid( + values=values, + transform=transform, + crs=crs, + bbox=BoundingBox( + min_lon=bounds.left, + min_lat=bounds.bottom, + max_lon=bounds.right, + max_lat=bounds.top, + ), + epoch=2021, + classes=dict(WORLDCOVER_CLASSES), + units="ESA WorldCover class code", + ) + + @pytest.fixture(scope="session") def ne_brazil_dem() -> DEM: """Live DEM at (-5.0, -39.0) with the default buffer. diff --git a/tests/fixtures/_generate_robit_bata_worldcover.py b/tests/fixtures/_generate_robit_bata_worldcover.py new file mode 100644 index 0000000..ac7675f --- /dev/null +++ b/tests/fixtures/_generate_robit_bata_worldcover.py @@ -0,0 +1,62 @@ +"""Regenerate the Robit Bata ESA WorldCover test fixture. + +Run this when ESA's tile changes or we move test bboxes: + + python tests/fixtures/_generate_robit_bata_worldcover.py + +The fixture is the same 0.075 deg patch used by the DEM/GHSL fixtures, +sampled at WorldCover 10 m (year 2021, v200). +""" + +import sys +from pathlib import Path + +PROJECT_ROOT = Path(__file__).resolve().parents[2] +if str(PROJECT_ROOT) not in sys.path: + sys.path.insert(0, str(PROJECT_ROOT)) + +import rasterio +from rasterio.profiles import default_gtiff_profile + +from floodpath.landuse import get_worldcover_landuse + +FIXTURE_PATH: Path = Path(__file__).resolve().parent / "robit_bata_worldcover.tif" +FIXTURE_LAT: float = 11.805 +FIXTURE_LON: float = 37.5625 +FIXTURE_BUFFER_DEG: float = 0.0375 +FIXTURE_YEAR: int = 2021 + + +def regenerate() -> None: + """Fetch the live WorldCover patch and overwrite the on-disk fixture.""" + grid = get_worldcover_landuse( + lat=FIXTURE_LAT, + lon=FIXTURE_LON, + buffer_deg=FIXTURE_BUFFER_DEG, + year=FIXTURE_YEAR, + ) + profile = default_gtiff_profile.copy() + profile.update( + height=grid.values.shape[0], + width=grid.values.shape[1], + count=1, + dtype=grid.values.dtype, + crs=grid.crs, + transform=grid.transform, + compress="lzw", + nodata=0, + ) + with rasterio.open(FIXTURE_PATH, "w", **profile) as dst: + dst.write(grid.values, 1) + size_kb = FIXTURE_PATH.stat().st_size / 1024 + print( + f"Wrote {FIXTURE_PATH} shape={grid.values.shape} " + f"year={grid.epoch} size={size_kb:.1f} KB" + ) + print("Class counts:") + for code, count in sorted(grid.class_counts().items()): + print(f" {code:>3} {grid.class_name(code):<28} {count}") + + +if __name__ == "__main__": + regenerate() diff --git a/tests/fixtures/robit_bata_worldcover.tif b/tests/fixtures/robit_bata_worldcover.tif new file mode 100644 index 0000000..91518b2 Binary files /dev/null and b/tests/fixtures/robit_bata_worldcover.tif differ diff --git a/tests/landuse/__init__.py b/tests/landuse/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/tests/landuse/test_fixture_regression.py b/tests/landuse/test_fixture_regression.py new file mode 100644 index 0000000..4c2bead --- /dev/null +++ b/tests/landuse/test_fixture_regression.py @@ -0,0 +1,35 @@ +"""Verify the committed WorldCover fixture still matches a live ESA fetch. + +Mirrors tests/exposure/test_fixture_regression.py: the single network-bound +contract guarding the offline landuse tests. If ESA reprocesses the +2021 v200 tile and the fixture drifts, this test fails first. +""" + +import numpy as np +import pytest + +from floodpath.landuse import get_worldcover_landuse +from floodpath.landuse.models import LanduseGrid +from tests.fixtures._generate_robit_bata_worldcover import ( + FIXTURE_BUFFER_DEG, + FIXTURE_LAT, + FIXTURE_LON, + FIXTURE_YEAR, +) + + +@pytest.mark.integration +class TestFixtureRegression: + def test_live_fetch_matches_committed_fixture( + self, + robit_bata_worldcover: LanduseGrid, + ) -> None: + live = get_worldcover_landuse( + lat=FIXTURE_LAT, + lon=FIXTURE_LON, + buffer_deg=FIXTURE_BUFFER_DEG, + year=FIXTURE_YEAR, + ) + assert live.shape == robit_bata_worldcover.shape + assert live.epoch == robit_bata_worldcover.epoch + np.testing.assert_array_equal(live.values, robit_bata_worldcover.values) diff --git a/tests/landuse/test_models.py b/tests/landuse/test_models.py new file mode 100644 index 0000000..6d09cf2 --- /dev/null +++ b/tests/landuse/test_models.py @@ -0,0 +1,86 @@ +"""Unit tests for landuse.models.""" + +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 + + +def _toy_grid() -> LanduseGrid: + """A 4x4 grid with three classes for arithmetic checks.""" + values = np.array( + [ + [10, 10, 40, 40], + [10, 40, 40, 50], + [10, 40, 40, 50], + [10, 40, 40, 50], + ], + dtype=np.uint8, + ) + return LanduseGrid( + values=values, + transform=Affine.translation(0, 0) * Affine.scale(0.001, -0.001), + crs="EPSG:4326", + bbox=BoundingBox(0.0, 0.0, 0.004, 0.004), + epoch=2021, + classes=dict(WORLDCOVER_CLASSES), + units="ESA WorldCover class code", + ) + + +class TestLanduseGrid: + def test_shape(self) -> None: + assert _toy_grid().shape == (4, 4) + + def test_class_counts(self) -> None: + assert _toy_grid().class_counts() == {10: 5, 40: 8, 50: 3} + + def test_class_name_known(self) -> None: + assert _toy_grid().class_name(40) == "cropland" + + def test_class_name_unknown(self) -> None: + assert _toy_grid().class_name(99) == "unknown" + + def test_fraction_present_class(self) -> None: + assert _toy_grid().fraction(40) == pytest.approx(8 / 16) + + def test_fraction_absent_class(self) -> None: + assert _toy_grid().fraction(20) == 0.0 + + +class TestRobitBataWorldCoverFixture: + """Pinned values from the committed Robit Bata WorldCover fixture (2021 v200).""" + + def test_shape_and_metadata(self, robit_bata_worldcover: LanduseGrid) -> None: + # 0.075 deg bbox at 1/12000 deg/pixel -> 900 cells per side. + assert robit_bata_worldcover.shape == (900, 900) + assert robit_bata_worldcover.epoch == 2021 + assert robit_bata_worldcover.crs == "EPSG:4326" + assert "WorldCover" in robit_bata_worldcover.units + + def test_class_counts_pinned(self, robit_bata_worldcover: LanduseGrid) -> None: + # Pinned: ESA WorldCover 2021 v200 class breakdown for the Robit + # Bata 0.075 deg patch (810,000 cells total). + assert robit_bata_worldcover.class_counts() == { + 10: 34877, # tree_cover + 20: 15902, # shrubland + 30: 46291, # grassland + 40: 698366, # cropland (dominant) + 50: 11774, # built_up + 60: 2774, # bare_or_sparse_vegetation + 80: 1, # permanent_water_bodies + 90: 15, # herbaceous_wetland + } + + def test_dominant_class_is_cropland( + self, + robit_bata_worldcover: LanduseGrid, + ) -> None: + # The Robit Bata watershed is overwhelmingly agricultural — cropland + # should dominate by a large margin. + assert robit_bata_worldcover.fraction(40) == pytest.approx( + 698366 / 810000, abs=1e-6 + ) diff --git a/tests/landuse/test_roughness.py b/tests/landuse/test_roughness.py new file mode 100644 index 0000000..9e78653 --- /dev/null +++ b/tests/landuse/test_roughness.py @@ -0,0 +1,164 @@ +"""Unit tests for landuse.roughness.""" + +import numpy as np +import pytest +from rasterio.transform import Affine + +from floodpath.dem.models import BoundingBox +from floodpath.landuse import ( + WORLDCOVER_CLASSES, + WORLDCOVER_MANNING_N, + LanduseGrid, + RoughnessGrid, + landuse_to_roughness, +) + + +def _toy_grid(values: np.ndarray) -> LanduseGrid: + """A LanduseGrid wrapper around a small numpy array for unit tests.""" + return LanduseGrid( + values=values, + 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=2021, + classes=dict(WORLDCOVER_CLASSES), + units="ESA WorldCover class code", + ) + + +class TestDefaultMapping: + def test_covers_all_worldcover_classes(self) -> None: + # Every WorldCover class code should have a default n. + assert set(WORLDCOVER_MANNING_N.keys()) == set(WORLDCOVER_CLASSES.keys()) + + def test_n_values_in_floodplain_range(self) -> None: + # Chow 1959 floodplain values fall in roughly [0.02, 0.20]. + for code, n in WORLDCOVER_MANNING_N.items(): + assert 0.02 <= n <= 0.20, f"class {code} has n={n} outside floodplain range" + + def test_tree_cover_is_roughest(self) -> None: + # Dense forest with undergrowth should be the highest-friction class. + max_code = max(WORLDCOVER_MANNING_N, key=WORLDCOVER_MANNING_N.get) + assert max_code == 10 # tree_cover + + def test_water_and_bare_are_smoothest(self) -> None: + # Open water and bare soil share the lowest n in the default table. + min_n = min(WORLDCOVER_MANNING_N.values()) + smoothest = {c for c, n in WORLDCOVER_MANNING_N.items() if n == min_n} + assert smoothest == {60, 80} # bare_or_sparse_vegetation, permanent_water + + +class TestLanduseToRoughness: + def test_known_codes_get_default_mapping(self) -> None: + # Tree cover (10) and cropland (40) -> 0.150 and 0.040 respectively. + toy = _toy_grid(np.array([[10, 40], [40, 10]], dtype=np.uint8)) + r = landuse_to_roughness(toy) + np.testing.assert_array_equal( + r.values, + np.array([[0.150, 0.040], [0.040, 0.150]], dtype=np.float32), + ) + + def test_unknown_code_gets_fallback(self) -> None: + # Code 200 is not in the default mapping -> default fallback 0.035. + toy = _toy_grid(np.array([[10, 200]], dtype=np.uint8)) + r = landuse_to_roughness(toy) + assert r.values[0, 0] == pytest.approx(0.150) + assert r.values[0, 1] == pytest.approx(0.035) + + def test_nodata_zero_gets_fallback(self) -> None: + # WorldCover nodata=0 is not a class -> fallback applied. + toy = _toy_grid(np.array([[0, 40]], dtype=np.uint8)) + r = landuse_to_roughness(toy, fallback=0.099) + assert r.values[0, 0] == pytest.approx(0.099) + assert r.values[0, 1] == pytest.approx(0.040) + + def test_custom_mapping_overrides_defaults(self) -> None: + # Site-specific calibration: bump cropland up, leave tree cover at default. + toy = _toy_grid(np.array([[10, 40]], dtype=np.uint8)) + r = landuse_to_roughness(toy, mapping={10: 0.150, 40: 0.060}) + assert r.values[0, 0] == pytest.approx(0.150) + assert r.values[0, 1] == pytest.approx(0.060) + # The mapping field records what was used. + assert r.mapping == {10: 0.150, 40: 0.060} + + def test_georef_is_preserved(self) -> None: + toy = _toy_grid(np.array([[10, 40]], dtype=np.uint8)) + r = landuse_to_roughness(toy) + assert r.transform == toy.transform + assert r.crs == toy.crs + assert r.bbox == toy.bbox + + def test_output_dtype_is_float32(self) -> None: + toy = _toy_grid(np.array([[10]], dtype=np.uint8)) + r = landuse_to_roughness(toy) + assert r.values.dtype == np.float32 + + def test_default_fallback_value(self) -> None: + toy = _toy_grid(np.array([[200]], dtype=np.uint8)) + r = landuse_to_roughness(toy) + # Cropland-like default for unknown classes. + assert r.fallback == pytest.approx(0.035) + assert r.values[0, 0] == pytest.approx(0.035) + + +class TestRoughnessGrid: + def test_shape(self) -> None: + toy = _toy_grid(np.array([[10, 40], [40, 10]], dtype=np.uint8)) + r = landuse_to_roughness(toy) + assert r.shape == (2, 2) + + def test_stats(self) -> None: + toy = _toy_grid(np.array([[10, 40], [40, 10]], dtype=np.uint8)) + r = landuse_to_roughness(toy) + s = r.stats() + assert s["min"] == pytest.approx(0.040) + assert s["max"] == pytest.approx(0.150) + assert s["mean"] == pytest.approx((0.150 + 0.040) / 2) + assert s["median"] == pytest.approx((0.150 + 0.040) / 2) + + def test_units(self) -> None: + toy = _toy_grid(np.array([[10]], dtype=np.uint8)) + r = landuse_to_roughness(toy) + assert "Manning" in r.units + + +class TestRobitBataRoughness: + """Pinned values from applying the default mapping to the committed + Robit Bata WorldCover fixture (810,000 cells, 2021 v200).""" + + def test_shape(self, robit_bata_worldcover: LanduseGrid) -> None: + r = landuse_to_roughness(robit_bata_worldcover) + assert r.shape == (900, 900) + + def test_stats_pinned(self, robit_bata_worldcover: LanduseGrid) -> None: + # Pinned: applying WORLDCOVER_MANNING_N (Chow/USACE) to the + # cropland-dominant Robit Bata patch yields: + # min 0.025 (water + bare/sparse) + # max 0.150 (tree_cover) + # mean 0.045476 (cropland-weighted) + # median 0.040 (cropland) + r = landuse_to_roughness(robit_bata_worldcover) + s = r.stats() + assert s["min"] == pytest.approx(0.025) + assert s["max"] == pytest.approx(0.150) + assert s["mean"] == pytest.approx(0.045476, abs=1e-5) + assert s["median"] == pytest.approx(0.040) + + def test_area_weighted_mean_matches_class_breakdown( + self, + robit_bata_worldcover: LanduseGrid, + ) -> None: + # Independent recomputation: sum_c (count_c * n_c) / total_cells. + # Validates that landuse_to_roughness applies the mapping cell-by-cell + # consistent with class composition. + counts = robit_bata_worldcover.class_counts() + weighted = sum( + cnt * WORLDCOVER_MANNING_N[code] + for code, cnt in counts.items() + if code in WORLDCOVER_MANNING_N + ) + expected_mean = weighted / sum(counts.values()) + + r = landuse_to_roughness(robit_bata_worldcover) + assert r.values.mean() == pytest.approx(expected_mean, abs=1e-5) diff --git a/tests/landuse/test_utils.py b/tests/landuse/test_utils.py new file mode 100644 index 0000000..cf0ece8 --- /dev/null +++ b/tests/landuse/test_utils.py @@ -0,0 +1,94 @@ +"""Unit tests for landuse.utils.""" + +import pytest + +from floodpath.dem.models import BoundingBox +from floodpath.landuse.constants import WORLDCOVER_BASE_URL +from floodpath.landuse.utils import ( + tiles_for_bbox, + worldcover_tile_id, + worldcover_tile_name, + worldcover_tile_url, +) + + +class TestWorldcoverTileId: + def test_robit_bata(self) -> None: + # Robit Bata at (11.805, 37.5625) snaps to the (9, 36) 3x3 tile. + assert worldcover_tile_id(11.805, 37.5625) == (9, 36) + + def test_negative_lat(self) -> None: + # (-5.5, 12.0) -> floor(-5.5/3)*3 = -6 lat, floor(12/3)*3 = 12 lon. + assert worldcover_tile_id(-5.5, 12.0) == (-6, 12) + + def test_negative_lon(self) -> None: + # (10.0, -39.0) -> 9 lat, -39 lon (already a 3-deg multiple). + assert worldcover_tile_id(10.0, -39.0) == (9, -39) + + def test_negative_lon_off_boundary(self) -> None: + # (10.0, -38.5) -> floor(-38.5/3)*3 = -39. + assert worldcover_tile_id(10.0, -38.5) == (9, -39) + + def test_exact_corner(self) -> None: + # The SW-corner (9, 36) belongs to its own tile, not the southern one. + assert worldcover_tile_id(9.0, 36.0) == (9, 36) + + +class TestWorldcoverTileName: + def test_northern_eastern(self) -> None: + assert worldcover_tile_name(9, 36) == "N09E036" + + def test_southern_western(self) -> None: + assert worldcover_tile_name(-12, -45) == "S12W045" + + def test_zero_pads_lat_two_lon_three(self) -> None: + assert worldcover_tile_name(3, 6) == "N03E006" + + +class TestTilesForBbox: + def test_single_tile(self) -> None: + bbox = BoundingBox(min_lon=37.4, min_lat=11.4, max_lon=37.6, max_lat=11.6) + assert tiles_for_bbox(bbox) == [(9, 36)] + + def test_corner_spanning_four_tiles(self) -> None: + # A 1-deg bbox crossing both a 3-deg lat boundary (12) and a + # 3-deg lon boundary (39) hits four neighbouring tiles. + bbox = BoundingBox(min_lon=38.5, min_lat=11.5, max_lon=39.5, max_lat=12.5) + assert sorted(tiles_for_bbox(bbox)) == [ + (9, 36), (9, 39), (12, 36), (12, 39), + ] + + def test_bbox_spanning_three_tile_rows_includes_interior_row(self) -> None: + # 7-deg lat span (covering 9, 12, 15) with 7-deg lon span + # (covering 36, 39, 42) -> 3x3 = 9 tiles, including interior (12, 39). + bbox = BoundingBox(min_lon=36.5, min_lat=9.5, max_lon=42.5, max_lat=15.5) + assert sorted(tiles_for_bbox(bbox)) == [ + (9, 36), (9, 39), (9, 42), + (12, 36), (12, 39), (12, 42), + (15, 36), (15, 39), (15, 42), + ] + + def test_bbox_aligned_exactly_to_tile_boundary(self) -> None: + # max edge sits exactly on a tile boundary -> exclude the tile beyond it. + bbox = BoundingBox(min_lon=36.0, min_lat=9.0, max_lon=39.0, max_lat=12.0) + assert tiles_for_bbox(bbox) == [(9, 36)] + + +class TestWorldcoverTileUrl: + def test_2021_url_uses_v200(self) -> None: + url = worldcover_tile_url(9, 36, year=2021) + assert url == ( + f"{WORLDCOVER_BASE_URL}/v200/2021/map/" + "ESA_WorldCover_10m_2021_v200_N09E036_Map.tif" + ) + + def test_2020_url_uses_v100(self) -> None: + url = worldcover_tile_url(9, 36, year=2020) + assert url == ( + f"{WORLDCOVER_BASE_URL}/v100/2020/map/" + "ESA_WorldCover_10m_2020_v100_N09E036_Map.tif" + ) + + def test_unknown_year_raises(self) -> None: + with pytest.raises(ValueError, match="year must be"): + worldcover_tile_url(9, 36, year=2019) diff --git a/tests/landuse/test_worldcover.py b/tests/landuse/test_worldcover.py new file mode 100644 index 0000000..f4f0f8b --- /dev/null +++ b/tests/landuse/test_worldcover.py @@ -0,0 +1,47 @@ +"""Unit tests for landuse.worldcover (offline + integration).""" + +import pytest + +from floodpath.landuse import get_worldcover_landuse + + +class TestGetWorldcoverOffline: + def test_invalid_year_raises(self) -> None: + # 2019 has no WorldCover release. + with pytest.raises(ValueError, match="year must be"): + get_worldcover_landuse(11.805, 37.5625, year=2019) + + +@pytest.mark.integration +class TestGetWorldcoverIntegration: + """Live test against the ESA WorldCover S3 bucket. + + Uses a tiny bbox over Djibouti (the same region exercised by the + WorldPop integration test) so the windowed read transfers a small + payload while still hitting the full /vsicurl/ + merge code path. + """ + + def test_djibouti_2021_returns_landuse(self) -> None: + grid = get_worldcover_landuse( + lat=11.59, + lon=43.15, + buffer_deg=0.05, + year=2021, + ) + assert grid.epoch == 2021 + assert grid.crs == "EPSG:4326" + # Djibouti city is dense — built-up class (50) should be present. + counts = grid.class_counts() + assert 50 in counts + assert counts[50] > 0 + + def test_2020_release_also_works(self) -> None: + grid = get_worldcover_landuse( + lat=11.59, + lon=43.15, + buffer_deg=0.05, + year=2020, + ) + assert grid.epoch == 2020 + # Built-up class still present in 2020 release. + assert grid.class_counts().get(50, 0) > 0