Add landuse module (ESA WorldCover) with Manning roughness

.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/

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

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",
+]

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",
+})

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"]

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),
+    )

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"
+    )