Add soil module (SoilGrids 2.0) with NEH 630 Ch7 hydrologic soil group

README.md

@@ -71,7 +71,8 @@ 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.landuse` | ESA WorldCover (10 m, AWS Open Data, COG) | 11-class land-cover raster (2020 v100, 2021 v200), Manning's roughness derivation |
+| `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.damage` | JRC Huizinga 2017 + DEM/HAND/GHSL | Per-cell flood depth and damage in m² of built-up surface |
 
 ## Depth-damage curves

floodpath/soil/__init__.py

@@ -0,0 +1,32 @@
+"""Soil module: hydrologic soil properties from SoilGrids 2.0.
+
+Provides the SOIL half of the SCS-CN runoff parameterisation. Pair with
+`floodpath.landuse` to derive a Curve Number raster.
+"""
+
+from .constants import (
+    HSG_CLASSES,
+    HSG_TO_INT,
+    INT_TO_HSG,
+    SOILGRIDS_DEFAULT_DEPTHS,
+    USDA_TEXTURE_CLASSES,
+    USDA_TEXTURE_TO_HSG,
+)
+from .hsg import texture_to_hsg
+from .models import HSGGrid, TextureGrid
+from .soilgrids import get_soilgrids_texture
+from .utils import usda_texture_class
+
+__all__ = [
+    "HSGGrid",
+    "HSG_CLASSES",
+    "HSG_TO_INT",
+    "INT_TO_HSG",
+    "SOILGRIDS_DEFAULT_DEPTHS",
+    "TextureGrid",
+    "USDA_TEXTURE_CLASSES",
+    "USDA_TEXTURE_TO_HSG",
+    "get_soilgrids_texture",
+    "texture_to_hsg",
+    "usda_texture_class",
+]

floodpath/soil/constants.py

@@ -0,0 +1,122 @@
+"""Constants for the soil module."""
+
+from types import MappingProxyType
+
+# --- ISRIC SoilGrids 2.0 (250 m global, COG via WebDAV) --------------------
+#
+# VRT files reference the underlying COG tiles in Goode Homolosine; GDAL's
+# /vsicurl/ + a WarpedVRT projects them to WGS84 on the fly. The bucket
+# honours HTTP range requests, so a window read is cheap.
+
+SOILGRIDS_BASE_URL: str = "https://files.isric.org/soilgrids/latest/data"
+SOILGRIDS_NATIVE_CRS: str = (
+    'PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["unknown",'
+    'ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]]],'
+    'PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]]],'
+    'CONVERSION["unknown",METHOD["Interrupted Goode Homolosine"],'
+    'PARAMETER["Longitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433]],'
+    'PARAMETER["False easting",0,LENGTHUNIT["metre",1]],'
+    'PARAMETER["False northing",0,LENGTHUNIT["metre",1]]],'
+    'CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1]],'
+    'AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1]]]'
+)
+SOILGRIDS_TARGET_CRS: str = "EPSG:4326"
+
+# SoilGrids stores sand/silt/clay as integer g/kg (0-1000). Multiply by 0.1
+# to get percent by mass (0-100).
+SOILGRIDS_TEXTURE_SCALE: float = 0.1
+
+# SoilGrids nodata sentinel (int16). Cells with this value have no
+# prediction (oceans, ice, glaciers, no-cover land). Must be masked
+# *before* applying the scale factor or it pollutes downstream stats.
+SOILGRIDS_NODATA: int = -32768
+
+# Standard depth slices (in cm) published by SoilGrids 2.0.
+SOILGRIDS_DEPTHS: tuple[str, ...] = (
+    "0-5cm",
+    "5-15cm",
+    "15-30cm",
+    "30-60cm",
+    "60-100cm",
+    "100-200cm",
+)
+# Default depth-weighted topsoil window (0-30 cm). Closest available proxy
+# for the 0-50 cm reference depth used in NEH 630 Ch7 Table 7-1.
+SOILGRIDS_DEFAULT_DEPTHS: tuple[str, ...] = ("0-5cm", "5-15cm", "15-30cm")
+# Width in cm of each depth slice — used as weights when averaging.
+SOILGRIDS_DEPTH_WIDTHS_CM: MappingProxyType = MappingProxyType({
+    "0-5cm": 5,
+    "5-15cm": 10,
+    "15-30cm": 15,
+    "30-60cm": 30,
+    "60-100cm": 40,
+    "100-200cm": 100,
+})
+
+SOILGRIDS_BUFFER_DEG: float = 0.1
+SOILGRIDS_TEXTURE_UNITS: str = "% by mass"
+
+
+# --- USDA Soil Texture Triangle (12 classes) -------------------------------
+#
+# Sand / silt / clay in % by mass. Boundaries follow the Soil Survey
+# Manual (1993) and Soil Taxonomy 2nd ed. (1999).
+
+USDA_TEXTURE_CLASSES: tuple[str, ...] = (
+    "sand",
+    "loamy_sand",
+    "sandy_loam",
+    "loam",
+    "silt_loam",
+    "silt",
+    "sandy_clay_loam",
+    "clay_loam",
+    "silty_clay_loam",
+    "sandy_clay",
+    "silty_clay",
+    "clay",
+)
+
+
+# --- USDA Texture class -> NEH 630 Ch7 Hydrologic Soil Group ---------------
+#
+# Source: NRCS National Engineering Handbook Part 630 Chapter 7
+# (Hydrologic Soil Groups, January 2009), §630.0701 textural criteria
+# in conjunction with §630.0701 prose ("typical" assignment for soils
+# with no shallow impermeable layer or high water table).
+#
+# §630.0701 group definitions:
+#   A  <10% clay AND >90% sand or gravel; gravel/sand textures.
+#   B  10-20% clay AND 50-90% sand; loamy sand or sandy loam textures.
+#   C  20-40% clay AND <50% sand; loam, silt loam, sandy clay loam,
+#      clay loam, silty clay loam.
+#   D  >40% clay AND <50% sand; clayey textures (clay, silty clay,
+#      sandy clay).
+#
+# A handful of texture classes straddle the prose ranges; we follow the
+# most-cited consensus mapping used in SWAT, HEC-HMS, and the Kalyanapu /
+# Yu / et-al. global SCS-CN literature.
+
+USDA_TEXTURE_TO_HSG: MappingProxyType = MappingProxyType({
+    "sand": "A",
+    "loamy_sand": "A",
+    "sandy_loam": "B",
+    "loam": "B",
+    "silt_loam": "B",
+    "silt": "B",
+    "sandy_clay_loam": "C",
+    "clay_loam": "C",
+    "silty_clay_loam": "C",
+    "sandy_clay": "D",
+    "silty_clay": "D",
+    "clay": "D",
+})
+
+HSG_CLASSES: tuple[str, ...] = ("A", "B", "C", "D")
+
+# Categorical encoding for raster storage. 0 = nodata.
+HSG_TO_INT: MappingProxyType = MappingProxyType({"A": 1, "B": 2, "C": 3, "D": 4})
+INT_TO_HSG: MappingProxyType = MappingProxyType({v: k for k, v in HSG_TO_INT.items()})
+
+HSG_NODATA: int = 0
+HSG_UNITS: str = "NEH 630 Ch7 hydrologic soil group (1=A, 2=B, 3=C, 4=D)"

floodpath/soil/hsg.py

@@ -0,0 +1,72 @@
+"""Map a TextureGrid to a Hydrologic Soil Group (HSG) raster.
+
+Per USDA NRCS NEH Part 630 Chapter 7 §630.0701 (January 2009): each cell's
+USDA texture class is looked up in `USDA_TEXTURE_TO_HSG` to assign A/B/C/D.
+This is the textural shortcut for the full Table 7-1 criteria (which would
+require Ksat, depth-to-impermeable-layer, and water-table data we do not
+have globally).
+
+The shortcut assumes "deep" soils with no shallow restrictive layer and
+no high water table — appropriate for global SCS-CN runoff modelling at
+~250 m. Users with site-calibrated Ksat data should compute HSG directly
+from Table 7-1.
+"""
+
+import numpy as np
+
+from .constants import (
+    HSG_NODATA,
+    HSG_TO_INT,
+    USDA_TEXTURE_TO_HSG,
+)
+from .models import HSGGrid, TextureGrid
+from .utils import usda_texture_class
+
+
+def texture_to_hsg(
+    texture: TextureGrid,
+    mapping: dict[str, str] | None = None,
+) -> HSGGrid:
+    """Convert a soil texture grid to a hydrologic soil group raster.
+
+    Args:
+        texture: Per-cell sand/silt/clay percentages (TextureGrid).
+        mapping: Optional override of the USDA-class -> HSG-letter table.
+            Defaults to USDA_TEXTURE_TO_HSG (NEH 630 Ch7 §630.0701).
+
+    Returns:
+        HSGGrid with int-encoded HSG classes (1=A, 2=B, 3=C, 4=D, 0=nodata).
+    """
+    table = dict(USDA_TEXTURE_TO_HSG) if mapping is None else dict(mapping)
+
+    # Cells with NaN in any texture component (SoilGrids nodata) cannot be
+    # classified — keep them at HSG_NODATA. Pass safe placeholders to the
+    # classifier and mask the output afterwards.
+    sand = np.where(np.isnan(texture.sand), 0.0, texture.sand)
+    silt = np.where(np.isnan(texture.silt), 0.0, texture.silt)
+    clay = np.where(np.isnan(texture.clay), 0.0, texture.clay)
+    nodata_mask = (
+        np.isnan(texture.sand)
+        | np.isnan(texture.silt)
+        | np.isnan(texture.clay)
+    )
+
+    classes = usda_texture_class(sand, silt, clay)
+    # `classes` is an object array of texture-class names; map each unique
+    # class to its HSG int and apply via fancy indexing.
+    hsg_int = np.full(classes.shape, HSG_NODATA, dtype=np.uint8)
+    for cls_name, letter in table.items():
+        if letter not in HSG_TO_INT:
+            raise ValueError(
+                f"mapping value {letter!r} is not a valid HSG letter"
+            )
+        hsg_int[classes == cls_name] = HSG_TO_INT[letter]
+    hsg_int[nodata_mask] = HSG_NODATA
+
+    return HSGGrid(
+        values=hsg_int,
+        transform=texture.transform,
+        crs=texture.crs,
+        bbox=texture.bbox,
+        depths=texture.depths,
+    )

floodpath/soil/models.py

@@ -0,0 +1,100 @@
+"""Dataclasses for the soil module."""
+
+from dataclasses import dataclass
+
+import numpy as np
+from rasterio.transform import Affine
+
+from floodpath.dem.models import BoundingBox
+
+from .constants import (
+    HSG_CLASSES,
+    HSG_TO_INT,
+    INT_TO_HSG,
+    SOILGRIDS_TEXTURE_UNITS,
+)
+
+
+@dataclass
+class TextureGrid:
+    """Per-cell soil texture composition (sand / silt / clay) as percentages.
+
+    Each of `sand`, `silt`, `clay` is a float32 array in `% by mass`. They
+    should sum to ~100 at every cell (SoilGrids predictions can be off by
+    1-2 due to independent regression — the texture-triangle classifier is
+    robust to that).
+
+    `depths` records which SoilGrids depth slice(s) the values represent;
+    when more than one was averaged, the list reflects the depth-weighted
+    composition.
+    """
+
+    sand: np.ndarray
+    silt: np.ndarray
+    clay: np.ndarray
+    transform: Affine
+    crs: str
+    bbox: BoundingBox
+    depths: tuple[str, ...]
+    units: str = SOILGRIDS_TEXTURE_UNITS
+
+    @property
+    def shape(self) -> tuple[int, int]:
+        """Return the (rows, cols) shape of the underlying rasters."""
+        return self.sand.shape  # type: ignore[return-value]
+
+    def closure_residual(self) -> np.ndarray:
+        """Return per-cell |sand+silt+clay - 100|, useful for diagnostics."""
+        return np.abs(self.sand + self.silt + self.clay - 100.0)
+
+
+@dataclass
+class HSGGrid:
+    """Per-cell NEH 630 Ch7 hydrologic soil group (A/B/C/D).
+
+    Stored as an integer raster with the encoding in HSG_TO_INT (1=A, 2=B,
+    3=C, 4=D, 0=nodata). The `epoch`-equivalent for soil is depth — recorded
+    in `depths` so consumers can introspect the assumption.
+    """
+
+    values: np.ndarray
+    transform: Affine
+    crs: str
+    bbox: BoundingBox
+    depths: tuple[str, ...]
+    units: str = "NEH 630 Ch7 HSG (1=A, 2=B, 3=C, 4=D)"
+
+    @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[str, int]:
+        """Count cells per HSG letter present in the grid (excludes nodata)."""
+        out: dict[str, int] = {}
+        for letter in HSG_CLASSES:
+            n = int(np.count_nonzero(self.values == HSG_TO_INT[letter]))
+            if n > 0:
+                out[letter] = n
+        return out
+
+    def fraction(self, letter: str) -> float:
+        """Return the fraction of cells in the given HSG letter."""
+        if self.values.size == 0:
+            return 0.0
+        n = int(np.count_nonzero(self.values == HSG_TO_INT[letter]))
+        return float(n / self.values.size)
+
+    def dominant_class(self) -> str:
+        """Return the HSG letter with the most cells; raises if grid is empty."""
+        counts = self.class_counts()
+        if not counts:
+            raise ValueError("HSGGrid contains no classified cells")
+        return max(counts, key=counts.get)  # type: ignore[arg-type]
+
+    def as_letter_grid(self) -> np.ndarray:
+        """Return a US1-string array with HSG letters; nodata becomes empty string."""
+        out = np.empty(self.values.shape, dtype="U1")
+        for code, letter in INT_TO_HSG.items():
+            out[self.values == code] = letter
+        return out