Microtopo specialist: 3DEP elevation + TWI + HAND

USGS 3DEP 1/3 arcsec DEM (resampled to 30 m for runtime), with two
hydrology indices computed offline by whitebox-workflows:

TWI — topographic wetness, ln(a/tan β), high where flow accumulates
HAND — height above nearest drainage, low where overbank flooding
pools first

These three rasters give the briefing terrain context that the
empirical/forward-modelled flood layers can't: even a parcel outside
every published zone can be in a topographic sink that any heavy rain
will fill.

.gitattributes routes *.tif through LFS.

.gitattributes

@@ -1,5 +1,6 @@
 # Riprap-specific LFS tracking
 *.geojson filter=lfs diff=lfs merge=lfs -text
+*.tif filter=lfs diff=lfs merge=lfs -text
 
 # Esri FileGDB internal binary files (DEP Stormwater scenario data)
 *.gdbtable filter=lfs diff=lfs merge=lfs -text

app/context/microtopo.py

@@ -0,0 +1,208 @@
+"""LiDAR/DEM-derived micro-topography specialist.
+
+Reads a window from a precomputed NYC-wide DEM (data/nyc_dem_30m.tif)
+fetched from USGS 3DEP via py3dep. Computes per-address terrain numbers
+that the static FEMA/DEP scenario maps don't expose.
+
+Metrics (all derived from the same small AOI raster):
+
+    point_elev_m          elevation at the address (m)
+    rel_elev_pct_750m     percentile of point elev in a 750-m radius
+    rel_elev_pct_200m     percentile of point elev in a 200-m radius
+                          (block-scale "is this a bowl?")
+    basin_relief_m        max-elev in 750-m AOI minus point elev
+    aoi_min_m, aoi_max_m  for context
+    resolution_m
+
+We deliberately stop at "shape-of-the-terrain" metrics rather than full
+hydrology — depression-fill / D8 flow accumulation on a flat coastal
+DEM are noisy and slow. Percentile + relief is what the reconciler
+actually needs to write a useful sentence.
+"""
+from __future__ import annotations
+
+import logging
+import warnings
+from dataclasses import dataclass
+from pathlib import Path
+
+import numpy as np
+
+warnings.filterwarnings("ignore")
+
+log = logging.getLogger("riprap.microtopo")
+
+DOC_ID = "microtopo"
+CITATION = "USGS 3DEP 30 m DEM (precomputed citywide GeoTIFF, WGS84)"
+
+DATA_DIR = Path(__file__).resolve().parent.parent.parent / "data"
+DEM_PATH = DATA_DIR / "nyc_dem_30m.tif"
+TWI_PATH = DATA_DIR / "twi.tif"
+HAND_PATH = DATA_DIR / "hand.tif"
+
+
+@dataclass
+class Microtopo:
+    point_elev_m: float
+    rel_elev_pct_750m: float    # 0..100
+    rel_elev_pct_200m: float    # 0..100
+    basin_relief_m: float
+    aoi_min_m: float
+    aoi_max_m: float
+    aoi_radius_m: int
+    resolution_m: int
+    # Hydrology indices computed on the same DEM (whitebox-workflows)
+    twi: float | None = None              # Topographic Wetness Index, ln(SCA / tan(slope))
+    hand_m: float | None = None           # Height Above Nearest Drainage (m)
+
+
+def _percentile_in_window(arr: np.ndarray, iy: int, ix: int, point_val: float,
+                          window_radius_cells: int) -> float:
+    H, W = arr.shape
+    y0 = max(0, iy - window_radius_cells)
+    y1 = min(H, iy + window_radius_cells + 1)
+    x0 = max(0, ix - window_radius_cells)
+    x1 = min(W, ix + window_radius_cells + 1)
+    sub = arr[y0:y1, x0:x1]
+    finite = sub[np.isfinite(sub)]
+    if finite.size == 0:
+        return float("nan")
+    return float((finite < point_val).sum()) / finite.size * 100.0
+
+
+_DEM_CACHE: dict = {}
+
+
+def _read_full_raster(path: Path) -> tuple[np.ndarray | None, dict | None]:
+    import rasterio
+    if not path.exists():
+        return None, None
+    with rasterio.open(path) as ds:
+        arr = ds.read(1).astype("float32")
+        nodata = ds.nodata
+        meta = {"H": ds.height, "W": ds.width,
+                "transform": ds.transform, "crs": ds.crs, "nodata": nodata}
+    if nodata is not None:
+        arr = np.where(arr == nodata, np.nan, arr)
+    return arr, meta
+
+
+def _load_dem():
+    """Read the precomputed NYC DEM + TWI + HAND rasters into memory.
+
+    All three are aligned (same grid, same transform). We hold them as
+    numpy arrays so per-query slicing is safe under threading.
+    """
+    if "arr" in _DEM_CACHE:
+        return _DEM_CACHE
+    arr, meta = _read_full_raster(DEM_PATH)
+    if arr is None:
+        log.warning("microtopo DEM not found at %s — run scripts/fetch_nyc_dem.py", DEM_PATH)
+        return None
+    twi, _   = _read_full_raster(TWI_PATH)
+    hand, _  = _read_full_raster(HAND_PATH)
+    _DEM_CACHE.update({
+        "arr": arr, "H": meta["H"], "W": meta["W"],
+        "transform": meta["transform"], "crs": meta["crs"],
+        "twi": twi, "hand": hand,
+    })
+    note = []
+    if twi is not None:  note.append(f"TWI {TWI_PATH.name}")
+    if hand is not None: note.append(f"HAND {HAND_PATH.name}")
+    log.info("microtopo: loaded NYC DEM %s (%dx%d, %s); aux: %s",
+             DEM_PATH.name, meta["H"], meta["W"], meta["crs"],
+             ", ".join(note) if note else "(none — algorithmic only)")
+    return _DEM_CACHE
+
+
+def warm():
+    _load_dem()
+
+
+def _row_col(transform, lat: float, lon: float) -> tuple[int, int]:
+    """Inverse-affine: WGS84 (lon,lat) -> raster (row, col).
+    Mirrors rasterio.transform.rowcol but without holding a dataset handle.
+    """
+    # affine: x = a*col + b*row + c ; y = d*col + e*row + f
+    # invert: col = (a_inv * (x - c)) approx — we have a diagonal affine
+    a, b, c, d, e, f = transform.a, transform.b, transform.c, transform.d, transform.e, transform.f
+    # diagonal case (b=d=0, common for north-up rasters):
+    col = int(round((lon - c) / a))
+    row = int(round((lat - f) / e))
+    return row, col
+
+
+def microtopo_at(lat: float, lon: float, radius_m: int = 750) -> Microtopo | None:
+    state = _load_dem()
+    if state is None:
+        return None
+    arr_full = state["arr"]
+    transform = state["transform"]
+
+    try:
+        row, col = _row_col(transform, lat, lon)
+    except Exception as e:
+        log.warning("microtopo index failed: %s", e)
+        return None
+
+    res_m = abs(transform.a) * 111_000.0 * np.cos(np.radians(lat))
+    cells_radius = max(2, int(np.ceil(radius_m / max(res_m, 1.0))))
+
+    H, W = state["H"], state["W"]
+    y0 = max(0, row - cells_radius); y1 = min(H, row + cells_radius + 1)
+    x0 = max(0, col - cells_radius); x1 = min(W, col + cells_radius + 1)
+    if y1 <= y0 or x1 <= x0:
+        return None
+
+    arr = arr_full[y0:y1, x0:x1].copy()
+
+    iy = row - y0
+    ix = col - x0
+    if not (0 <= iy < arr.shape[0] and 0 <= ix < arr.shape[1]):
+        return None
+
+    point_elev = float(arr[iy, ix])
+    if not np.isfinite(point_elev):
+        for r in range(1, 6):
+            ya, yb = max(0, iy - r), min(arr.shape[0], iy + r + 1)
+            xa, xb = max(0, ix - r), min(arr.shape[1], ix + r + 1)
+            sub = arr[ya:yb, xa:xb]
+            if np.isfinite(sub).any():
+                point_elev = float(np.nanmean(sub))
+                break
+        else:
+            return None
+
+    finite = arr[np.isfinite(arr)]
+    if finite.size == 0:
+        return None
+    aoi_min = float(finite.min())
+    aoi_max = float(finite.max())
+
+    pct_750 = float((finite < point_elev).sum()) / finite.size * 100.0
+    cells_200m = max(1, int(round(200 / max(res_m, 1.0))))
+    pct_200 = _percentile_in_window(arr, iy, ix, point_elev, cells_200m)
+
+    twi_arr = state.get("twi")
+    hand_arr = state.get("hand")
+    twi_v: float | None = None
+    hand_v: float | None = None
+    if twi_arr is not None and 0 <= row < H and 0 <= col < W:
+        v = float(twi_arr[row, col])
+        twi_v = round(v, 2) if np.isfinite(v) else None
+    if hand_arr is not None and 0 <= row < H and 0 <= col < W:
+        v = float(hand_arr[row, col])
+        hand_v = round(v, 2) if np.isfinite(v) else None
+
+    return Microtopo(
+        point_elev_m=round(point_elev, 2),
+        rel_elev_pct_750m=round(pct_750, 1),
+        rel_elev_pct_200m=round(pct_200, 1),
+        basin_relief_m=round(aoi_max - point_elev, 2),
+        aoi_min_m=round(aoi_min, 2),
+        aoi_max_m=round(aoi_max, 2),
+        aoi_radius_m=radius_m,
+        resolution_m=int(round(res_m)),
+        twi=twi_v,
+        hand_m=hand_v,
+    )