diff --git a/examples/user_guide/33_Vegetation_Flood_Modeling.ipynb b/examples/user_guide/33_Vegetation_Flood_Modeling.ipynb new file mode 100644 index 00000000..6d295627 --- /dev/null +++ b/examples/user_guide/33_Vegetation_Flood_Modeling.ipynb @@ -0,0 +1,791 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "title", + "metadata": {}, + "source": [ + "# Vegetation-Aware Flood Modeling\n", + "\n", + "Vegetation controls how water moves across a landscape. Dense forest slows overland flow and holds water deeper locally, while bare or paved surfaces let water move fast and shallow. The SCS curve number method captures the other side: vegetation and soil type together determine how much rainfall infiltrates vs. runs off.\n", + "\n", + "This notebook compares the standard flood analysis tools (uniform roughness, manually assigned curve numbers) with the new vegetation-aware functions that derive hydraulic parameters directly from land cover or NDVI data:\n", + "\n", + "| Standard approach | Vegetation-aware approach |\n", + "|---|---|\n", + "| `travel_time(fl, slope, mannings_n=0.03)` | `vegetation_roughness()` → `travel_time()` |\n", + "| `curve_number_runoff(rain, curve_number=80)` | `vegetation_curve_number()` → `curve_number_runoff()` |\n", + "| `flood_depth(hand, water_level=10)` | `flood_depth_vegetation(hand, slope, n, q)` |\n", + "\n", + "We'll use the same Copernicus DEM tile as the [flood analysis guide](12_Flood_Analysis.ipynb), with synthetic NLCD land cover and soil group rasters to demonstrate the full pipeline." + ] + }, + { + "cell_type": "markdown", + "id": "toc", + "metadata": {}, + "source": [ + "[Setup](#Setup) · [Vegetation roughness](#Vegetation-roughness) · [Travel time comparison](#Travel-time-comparison) · [Curve number from land cover](#Curve-number-from-land-cover) · [Runoff comparison](#Runoff-comparison) · [Vegetation-adjusted flood depth](#Vegetation-adjusted-flood-depth) · [Full pipeline comparison](#Full-pipeline-comparison)" + ] + }, + { + "cell_type": "markdown", + "id": "setup-md", + "metadata": {}, + "source": [ + "## Setup\n", + "\n", + "Load a DEM, compute the hydrology stack (fill, flow direction, accumulation, HAND), and build synthetic land cover / soil rasters." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "imports", + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "import numpy as np\n", + "import xarray as xr\n", + "\n", + "import matplotlib.pyplot as plt\n", + "from matplotlib.colors import ListedColormap, LinearSegmentedColormap, LogNorm, BoundaryNorm\n", + "from matplotlib.patches import Patch\n", + "\n", + "import xrspatial\n", + "from xrspatial.flood import (\n", + " flood_depth, inundation, curve_number_runoff, travel_time,\n", + " vegetation_roughness, vegetation_curve_number, flood_depth_vegetation,\n", + " NLCD_MANNINGS_N, NLCD_CURVE_NUMBER,\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "load-dem", + "metadata": {}, + "outputs": [], + "source": [ + "try:\n", + " import rasterio\n", + " from rasterio.windows import Window\n", + "\n", + " url = (\n", + " \"https://copernicus-dem-30m.s3.amazonaws.com/\"\n", + " \"Copernicus_DSM_COG_10_N46_00_W123_00_DEM/\"\n", + " \"Copernicus_DSM_COG_10_N46_00_W123_00_DEM.tif\"\n", + " )\n", + "\n", + " with rasterio.open(url) as src:\n", + " window = Window(col_off=2400, row_off=2400, width=600, height=600)\n", + " data = src.read(1, window=window).astype(np.float64)\n", + " nodata = src.nodata\n", + "\n", + " if nodata is not None:\n", + " data[data == nodata] = np.nan\n", + "\n", + " H, W = data.shape\n", + " dem = xr.DataArray(data, dims=['y', 'x'], name='elevation',\n", + " attrs={'res': (1, 1)})\n", + " dem['y'] = np.linspace(H - 1, 0, H)\n", + " dem['x'] = np.linspace(0, W - 1, W)\n", + " print(f\"Loaded Copernicus 30m DEM: {dem.shape}\")\n", + "\n", + "except Exception as e:\n", + " print(f\"Remote DEM unavailable ({e}), generating synthetic terrain\")\n", + " H, W = 600, 600\n", + " dem = xr.DataArray(np.zeros((H, W)), dims=['y', 'x'])\n", + " dem = dem.xrs.generate_terrain(seed=10)\n", + " dem.name = 'elevation'\n", + " print(f\"Generated terrain: {dem.shape}\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "hydro-stack", + "metadata": {}, + "outputs": [], + "source": [ + "# Fill, resolve flats, flow direction, accumulation, HAND\n", + "dem_filled = xrspatial.fill(dem)\n", + "rng = np.random.RandomState(42)\n", + "dem_filled.values += rng.uniform(0, 0.001, dem_filled.shape)\n", + "dem_filled = xrspatial.fill(dem_filled)\n", + "rng2 = np.random.RandomState(123)\n", + "dem_filled.values += rng2.uniform(0, 1e-6, dem_filled.shape)\n", + "\n", + "flow_dir = xrspatial.flow_direction(dem_filled)\n", + "flow_accum = xrspatial.flow_accumulation(flow_dir)\n", + "\n", + "threshold = 200\n", + "hand_raster = xrspatial.hand(flow_dir, flow_accum, dem_filled, threshold=threshold)\n", + "slope_raster = xrspatial.slope(dem_filled)\n", + "fl_raster = xrspatial.flow_length(flow_dir)\n", + "hillshade = xrspatial.hillshade(dem)\n", + "\n", + "# Stream network for overlays\n", + "streams = xr.where(flow_accum >= threshold, flow_accum, np.nan)\n", + "streams.name = 'streams'\n", + "\n", + "print(f\"HAND range: {np.nanmin(hand_raster.values):.1f} to {np.nanmax(hand_raster.values):.1f} m\")\n", + "print(f\"Slope range: {np.nanmin(slope_raster.values):.2f} to {np.nanmax(slope_raster.values):.1f} deg\")" + ] + }, + { + "cell_type": "markdown", + "id": "synth-lc-md", + "metadata": {}, + "source": [ + "### Synthetic land cover and soil rasters\n", + "\n", + "We build a plausible NLCD land cover raster from the terrain itself: valley floors get developed classes, mid-slopes get forest, ridgetops get grassland/shrub. This is crude but gives spatially varying roughness and curve numbers that track the terrain in a realistic way.\n", + "\n", + "The soil group raster assigns sandy soil (group A) on ridges and clay (group D) in valleys, reflecting typical catena patterns." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "synth-lc", + "metadata": {}, + "outputs": [], + "source": [ + "hand_vals = hand_raster.fillna(0).values\n", + "\n", + "# NLCD classes based on HAND (proximity to drainage)\n", + "nlcd_data = np.full((H, W), 71, dtype=np.int32) # default: Grassland\n", + "nlcd_data[hand_vals < 5] = 24 # Developed, High Intensity (floodplain)\n", + "nlcd_data[(hand_vals >= 5) & (hand_vals < 15)] = 21 # Developed, Open Space\n", + "nlcd_data[(hand_vals >= 15) & (hand_vals < 40)] = 41 # Deciduous Forest\n", + "nlcd_data[(hand_vals >= 40) & (hand_vals < 80)] = 42 # Evergreen Forest\n", + "nlcd_data[hand_vals >= 80] = 52 # Shrub/Scrub (ridgetops)\n", + "\n", + "nlcd_raster = xr.DataArray(nlcd_data, dims=['y', 'x'], name='nlcd',\n", + " coords=dem.coords, attrs=dem.attrs)\n", + "\n", + "# Soil group: valleys=D (clay), mid=C/B, ridges=A (sandy)\n", + "sg_data = np.full((H, W), 2, dtype=np.int32) # default: B\n", + "sg_data[hand_vals < 10] = 4 # D (clay in valley bottoms)\n", + "sg_data[(hand_vals >= 10) & (hand_vals < 30)] = 3 # C\n", + "sg_data[(hand_vals >= 30) & (hand_vals < 60)] = 2 # B\n", + "sg_data[hand_vals >= 60] = 1 # A (sandy ridgetops)\n", + "\n", + "sg_raster = xr.DataArray(sg_data, dims=['y', 'x'], name='soil_group',\n", + " coords=dem.coords, attrs=dem.attrs)\n", + "\n", + "# Show the land cover map\n", + "nlcd_labels = {24: 'Dev High', 21: 'Dev Open', 41: 'Deciduous',\n", + " 42: 'Evergreen', 52: 'Shrub', 71: 'Grassland'}\n", + "nlcd_colors = {24: '#eb1a1d', 21: '#e8d1aa', 41: '#6ca966',\n", + " 42: '#1d6533', 52: '#ccb879', 71: '#e8e87b'}\n", + "\n", + "codes_in_order = [24, 21, 41, 42, 52, 71]\n", + "cmap_nlcd = ListedColormap([nlcd_colors[c] for c in codes_in_order])\n", + "bounds = [0] + [codes_in_order[i] + 0.5 for i in range(len(codes_in_order))]\n", + "# Map codes to sequential integers for plotting\n", + "nlcd_plot = np.full_like(nlcd_data, np.nan, dtype=np.float64)\n", + "for i, code in enumerate(codes_in_order):\n", + " nlcd_plot[nlcd_data == code] = i\n", + "nlcd_plot_da = xr.DataArray(nlcd_plot, dims=['y', 'x'], coords=dem.coords)\n", + "\n", + "fig, axes = plt.subplots(1, 2, figsize=(16, 6))\n", + "\n", + "# Land cover\n", + "ax = axes[0]\n", + "hillshade.plot.imshow(ax=ax, cmap='gray', add_colorbar=False)\n", + "nlcd_plot_da.plot.imshow(ax=ax, cmap=cmap_nlcd, alpha=0.7, add_colorbar=False,\n", + " vmin=-0.5, vmax=len(codes_in_order) - 0.5)\n", + "ax.legend(handles=[Patch(facecolor=nlcd_colors[c], label=nlcd_labels[c])\n", + " for c in codes_in_order],\n", + " loc='lower right', fontsize=9, framealpha=0.9)\n", + "ax.set_title('Synthetic NLCD land cover')\n", + "ax.set_axis_off()\n", + "\n", + "# Soil group\n", + "ax = axes[1]\n", + "sg_cmap = ListedColormap(['#fee5d9', '#fcae91', '#fb6a4a', '#cb181d'])\n", + "sg_plot = xr.DataArray(sg_data.astype(np.float64), dims=['y', 'x'], coords=dem.coords)\n", + "hillshade.plot.imshow(ax=ax, cmap='gray', add_colorbar=False)\n", + "sg_plot.plot.imshow(ax=ax, cmap=sg_cmap, alpha=0.7, add_colorbar=False,\n", + " vmin=0.5, vmax=4.5)\n", + "ax.legend(handles=[Patch(facecolor=c, label=f'Group {l}')\n", + " for c, l in zip(['#fee5d9', '#fcae91', '#fb6a4a', '#cb181d'],\n", + " ['A (sandy)', 'B', 'C', 'D (clay)'])],\n", + " loc='lower right', fontsize=9, framealpha=0.9)\n", + "ax.set_title('Hydrologic soil group')\n", + "ax.set_axis_off()\n", + "\n", + "plt.tight_layout()" + ] + }, + { + "cell_type": "markdown", + "id": "veg-roughness-md", + "metadata": {}, + "source": [ + "## Vegetation roughness\n", + "\n", + "`vegetation_roughness` converts land cover or NDVI data into a Manning's n raster. Two modes are available:\n", + "\n", + "- **`mode='nlcd'`**: Categorical lookup from NLCD codes. Each code maps to a literature-standard roughness value (Chow 1959, Arcement & Schneider 1989). Unrecognized codes produce NaN.\n", + "- **`mode='ndvi'`**: Piecewise linear interpolation. Bare ground (NDVI < 0.1) gets low n (~0.02), dense canopy (NDVI > 0.6) gets high n (~0.16). Less precise than categorical, but works with any multispectral imagery.\n", + "\n", + "Compare this with the standard approach of picking a single uniform n value." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "veg-roughness-nlcd", + "metadata": {}, + "outputs": [], + "source": [ + "# Derive spatially varying Manning's n from NLCD\n", + "n_raster = vegetation_roughness(nlcd_raster, mode='nlcd')\n", + "\n", + "n_vals = n_raster.values[~np.isnan(n_raster.values)]\n", + "print(f\"Manning's n from NLCD:\")\n", + "print(f\" min: {n_vals.min():.3f}\")\n", + "print(f\" median: {np.median(n_vals):.3f}\")\n", + "print(f\" max: {n_vals.max():.3f}\")\n", + "print()\n", + "\n", + "# Print the lookup table\n", + "print(\"NLCD code -> Manning's n:\")\n", + "for code in sorted(NLCD_MANNINGS_N):\n", + " count = np.sum(nlcd_data == code)\n", + " if count > 0:\n", + " print(f\" {code:2d}: n={NLCD_MANNINGS_N[code]:.3f} ({count:>6d} cells)\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "veg-roughness-plot", + "metadata": {}, + "outputs": [], + "source": [ + "n_cmap = LinearSegmentedColormap.from_list('roughness',\n", + " ['#ffffcc', '#a1dab4', '#41b6c4', '#225ea8'])\n", + "\n", + "fig, axes = plt.subplots(1, 2, figsize=(16, 6))\n", + "\n", + "# Uniform n\n", + "ax = axes[0]\n", + "uniform_n = xr.DataArray(np.full((H, W), 0.03, dtype=np.float64),\n", + " dims=['y', 'x'], coords=dem.coords)\n", + "hillshade.plot.imshow(ax=ax, cmap='gray', add_colorbar=False)\n", + "uniform_n.plot.imshow(ax=ax, cmap=n_cmap, alpha=0.7, add_colorbar=True,\n", + " vmin=0.02, vmax=0.16,\n", + " cbar_kwargs={'label': \"Manning's n\", 'shrink': 0.7})\n", + "ax.set_title('Standard: uniform n = 0.03')\n", + "ax.set_axis_off()\n", + "\n", + "# Vegetation-derived n\n", + "ax = axes[1]\n", + "hillshade.plot.imshow(ax=ax, cmap='gray', add_colorbar=False)\n", + "n_raster.plot.imshow(ax=ax, cmap=n_cmap, alpha=0.7, add_colorbar=True,\n", + " vmin=0.02, vmax=0.16,\n", + " cbar_kwargs={'label': \"Manning's n\", 'shrink': 0.7})\n", + "ax.set_title('Vegetation-aware: n from NLCD')\n", + "ax.set_axis_off()\n", + "\n", + "plt.tight_layout()" + ] + }, + { + "cell_type": "markdown", + "id": "ndvi-mode-md", + "metadata": {}, + "source": [ + "### NDVI mode\n", + "\n", + "When you have multispectral imagery instead of a classified land cover map, use `mode='ndvi'`. We'll generate a synthetic NDVI raster that correlates with elevation (higher = sparser vegetation) to demonstrate." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ndvi-mode", + "metadata": {}, + "outputs": [], + "source": [ + "# Synthetic NDVI: dense vegetation in valleys, sparse on ridges\n", + "elev_norm = (dem_filled.values - np.nanmin(dem_filled.values)) / \\\n", + " (np.nanmax(dem_filled.values) - np.nanmin(dem_filled.values))\n", + "ndvi_data = np.clip(0.7 - 0.5 * elev_norm + np.random.RandomState(99).normal(0, 0.05, (H, W)),\n", + " -0.1, 0.9).astype(np.float64)\n", + "ndvi_raster = xr.DataArray(ndvi_data, dims=['y', 'x'], name='ndvi',\n", + " coords=dem.coords, attrs=dem.attrs)\n", + "\n", + "n_from_ndvi = vegetation_roughness(ndvi_raster, mode='ndvi')\n", + "\n", + "fig, axes = plt.subplots(1, 2, figsize=(16, 6))\n", + "\n", + "ax = axes[0]\n", + "ndvi_cmap = LinearSegmentedColormap.from_list('ndvi', ['#8c510a', '#d8b365', '#f6e8c3',\n", + " '#c7eae5', '#5ab4ac', '#01665e'])\n", + "ndvi_raster.plot.imshow(ax=ax, cmap=ndvi_cmap, add_colorbar=True,\n", + " cbar_kwargs={'label': 'NDVI', 'shrink': 0.7})\n", + "ax.set_title('Synthetic NDVI')\n", + "ax.set_axis_off()\n", + "\n", + "ax = axes[1]\n", + "hillshade.plot.imshow(ax=ax, cmap='gray', add_colorbar=False)\n", + "n_from_ndvi.plot.imshow(ax=ax, cmap=n_cmap, alpha=0.7, add_colorbar=True,\n", + " vmin=0.02, vmax=0.16,\n", + " cbar_kwargs={'label': \"Manning's n\", 'shrink': 0.7})\n", + "ax.set_title(\"Manning's n from NDVI\")\n", + "ax.set_axis_off()\n", + "\n", + "plt.tight_layout()" + ] + }, + { + "cell_type": "markdown", + "id": "tt-compare-md", + "metadata": {}, + "source": [ + "## Travel time comparison\n", + "\n", + "The existing `travel_time` function accepts Manning's n as either a scalar or a DataArray. With `vegetation_roughness`, the n raster plugs in directly. The effect is pronounced: forested slopes slow water down, while developed areas near channels let it move faster.\n", + "\n", + "We compare the time of concentration (maximum travel time across the grid) for three scenarios:\n", + "1. Uniform n = 0.03 (channel-like, standard assumption)\n", + "2. Uniform n = 0.10 (forested, conservative assumption)\n", + "3. Spatially varying n from NLCD (the vegetation-aware approach)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "tt-compare", + "metadata": {}, + "outputs": [], + "source": [ + "tt_uniform_low = travel_time(fl_raster, slope_raster, mannings_n=0.03)\n", + "tt_uniform_high = travel_time(fl_raster, slope_raster, mannings_n=0.10)\n", + "tt_vegetation = travel_time(fl_raster, slope_raster, mannings_n=n_raster)\n", + "\n", + "for label, tt in [('Uniform n=0.03', tt_uniform_low),\n", + " ('Uniform n=0.10', tt_uniform_high),\n", + " ('NLCD vegetation', tt_vegetation)]:\n", + " vals = tt.values[np.isfinite(tt.values)]\n", + " print(f\"{label:>20s}: Tc = {vals.max():.1f}, median = {np.median(vals):.1f}\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "tt-compare-plot", + "metadata": {}, + "outputs": [], + "source": [ + "tt_cmap = LinearSegmentedColormap.from_list('tt', ['#08306b', '#2171b5', '#fee08b', '#d73027'])\n", + "\n", + "# Shared color scale from the vegetation result\n", + "tt_finite = tt_vegetation.values[np.isfinite(tt_vegetation.values)]\n", + "vmin, vmax = np.nanpercentile(tt_finite, [2, 98])\n", + "\n", + "fig, axes = plt.subplots(1, 3, figsize=(20, 6))\n", + "for ax, tt, title in zip(axes,\n", + " [tt_uniform_low, tt_uniform_high, tt_vegetation],\n", + " ['Uniform n=0.03', 'Uniform n=0.10', 'Vegetation n (NLCD)']):\n", + " hillshade.plot.imshow(ax=ax, cmap='gray', add_colorbar=False)\n", + " tt.plot.imshow(ax=ax, cmap=tt_cmap, alpha=0.75, vmin=vmin, vmax=vmax,\n", + " add_colorbar=True,\n", + " cbar_kwargs={'label': 'Travel time', 'shrink': 0.7})\n", + " ax.set_title(title)\n", + " ax.set_axis_off()\n", + "\n", + "plt.tight_layout()" + ] + }, + { + "cell_type": "markdown", + "id": "tt-note", + "metadata": {}, + "source": [ + "The vegetation-aware result splits the difference: developed valley floors move water fast (like n=0.03), while forested mid-slopes slow it down (closer to n=0.10). This gives a more realistic spatial pattern than either uniform assumption." + ] + }, + { + "cell_type": "markdown", + "id": "cn-md", + "metadata": {}, + "source": [ + "## Curve number from land cover\n", + "\n", + "`vegetation_curve_number` automates the NLCD + soil group → CN lookup that the [flood analysis guide](12_Flood_Analysis.ipynb) did manually. The lookup table contains 60 entries (15 NLCD classes × 4 soil groups) from TR-55 and HEC-HMS.\n", + "\n", + "Compare:\n", + "- **Before**: manually build a CN raster with `xr.where` (error-prone, limited to a few classes)\n", + "- **After**: one function call that handles all 15 NLCD classes and 4 soil groups" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cn-compute", + "metadata": {}, + "outputs": [], + "source": [ + "cn_raster = vegetation_curve_number(nlcd_raster, sg_raster)\n", + "\n", + "cn_vals = cn_raster.values[~np.isnan(cn_raster.values)]\n", + "print(f\"Curve number from NLCD + soil group:\")\n", + "print(f\" min: {cn_vals.min():.0f}\")\n", + "print(f\" median: {np.median(cn_vals):.0f}\")\n", + "print(f\" max: {cn_vals.max():.0f}\")\n", + "print()\n", + "\n", + "# Compare with the manual approach from the flood analysis guide\n", + "hand_filled = hand_raster.fillna(0)\n", + "cn_manual = xr.where(hand_filled < 20, 85.0, 55.0).astype(np.float64)\n", + "cn_manual.name = 'curve_number_manual'\n", + "\n", + "print(f\"Manual CN (2-class): {np.unique(cn_manual.values)}\")\n", + "print(f\"Vegetation CN range: {cn_vals.min():.0f} to {cn_vals.max():.0f} ({len(np.unique(cn_vals))} distinct values)\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cn-plot", + "metadata": {}, + "outputs": [], + "source": [ + "cn_cmap = LinearSegmentedColormap.from_list('cn', ['#1a9850', '#fee08b', '#d73027'])\n", + "\n", + "fig, axes = plt.subplots(1, 2, figsize=(16, 6))\n", + "\n", + "ax = axes[0]\n", + "hillshade.plot.imshow(ax=ax, cmap='gray', add_colorbar=False)\n", + "cn_manual_da = xr.DataArray(cn_manual.values, dims=['y', 'x'], coords=dem.coords)\n", + "cn_manual_da.plot.imshow(ax=ax, cmap=cn_cmap, alpha=0.7, vmin=25, vmax=100,\n", + " add_colorbar=True,\n", + " cbar_kwargs={'label': 'Curve Number', 'shrink': 0.7})\n", + "ax.set_title('Manual: 2-class CN (55 / 85)')\n", + "ax.set_axis_off()\n", + "\n", + "ax = axes[1]\n", + "hillshade.plot.imshow(ax=ax, cmap='gray', add_colorbar=False)\n", + "cn_raster.plot.imshow(ax=ax, cmap=cn_cmap, alpha=0.7, vmin=25, vmax=100,\n", + " add_colorbar=True,\n", + " cbar_kwargs={'label': 'Curve Number', 'shrink': 0.7})\n", + "ax.set_title('Vegetation-aware: CN from NLCD + soil')\n", + "ax.set_axis_off()\n", + "\n", + "plt.tight_layout()" + ] + }, + { + "cell_type": "markdown", + "id": "runoff-md", + "metadata": {}, + "source": [ + "## Runoff comparison\n", + "\n", + "Feed both CN rasters into `curve_number_runoff` with a uniform 100 mm storm. The vegetation-aware CN produces a more nuanced runoff pattern: developed areas on clay soils (high CN) generate much more runoff than forested ridges on sandy soil (low CN)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "runoff-compute", + "metadata": {}, + "outputs": [], + "source": [ + "rainfall = xr.DataArray(\n", + " np.full((H, W), 100.0, dtype=np.float64),\n", + " dims=['y', 'x'], coords=dem.coords, attrs=dem.attrs, name='rainfall',\n", + ")\n", + "\n", + "runoff_manual = curve_number_runoff(rainfall, curve_number=cn_manual_da)\n", + "runoff_veg = curve_number_runoff(rainfall, curve_number=cn_raster)\n", + "\n", + "for label, ro in [('Manual 2-class CN', runoff_manual),\n", + " ('Vegetation CN', runoff_veg)]:\n", + " vals = ro.values[~np.isnan(ro.values)]\n", + " print(f\"{label:>20s}: mean={vals.mean():.1f} mm, \"\n", + " f\"min={vals.min():.1f} mm, max={vals.max():.1f} mm\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "runoff-plot", + "metadata": {}, + "outputs": [], + "source": [ + "runoff_cmap = LinearSegmentedColormap.from_list('runoff', ['#ffffcc', '#fd8d3c', '#800026'])\n", + "\n", + "fig, axes = plt.subplots(1, 2, figsize=(16, 6))\n", + "\n", + "ax = axes[0]\n", + "hillshade.plot.imshow(ax=ax, cmap='gray', add_colorbar=False)\n", + "runoff_manual.plot.imshow(ax=ax, cmap=runoff_cmap, alpha=0.8, vmin=0, vmax=80,\n", + " add_colorbar=True,\n", + " cbar_kwargs={'label': 'Runoff (mm)', 'shrink': 0.7})\n", + "ax.set_title('Runoff: manual 2-class CN')\n", + "ax.set_axis_off()\n", + "\n", + "ax = axes[1]\n", + "hillshade.plot.imshow(ax=ax, cmap='gray', add_colorbar=False)\n", + "runoff_veg.plot.imshow(ax=ax, cmap=runoff_cmap, alpha=0.8, vmin=0, vmax=80,\n", + " add_colorbar=True,\n", + " cbar_kwargs={'label': 'Runoff (mm)', 'shrink': 0.7})\n", + "ax.set_title('Runoff: vegetation-aware CN')\n", + "ax.set_axis_off()\n", + "\n", + "plt.tight_layout()" + ] + }, + { + "cell_type": "markdown", + "id": "fdv-md", + "metadata": {}, + "source": [ + "## Vegetation-adjusted flood depth\n", + "\n", + "The standard `flood_depth` uses a uniform water level: every channel gets the same stage. `flood_depth_vegetation` computes depth from Manning's normal depth equation instead, which factors in slope and roughness:\n", + "\n", + "$$h = \\left(\\frac{q \\cdot n}{\\sqrt{\\tan(S)}}\\right)^{3/5}$$\n", + "\n", + "where $q$ is unit discharge (m²/s), $n$ is Manning's roughness, and $S$ is slope in degrees. The flood depth is $h - \\text{HAND}$ where $h > \\text{HAND}$.\n", + "\n", + "The key difference: rougher surfaces (denser vegetation) produce deeper, slower water locally. Steeper slopes produce shallower, faster water. This gives a physically grounded depth estimate that varies with both terrain and land cover." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fdv-compute", + "metadata": {}, + "outputs": [], + "source": [ + "# Standard flood depth at 10 m water level\n", + "depth_standard = flood_depth(hand_raster, water_level=10)\n", + "\n", + "# Vegetation-adjusted flood depth\n", + "# q=2.0 m^2/s is a moderate unit discharge\n", + "depth_veg = flood_depth_vegetation(hand_raster, slope_raster, n_raster,\n", + " unit_discharge=2.0)\n", + "\n", + "for label, d in [('Standard (wl=10m)', depth_standard),\n", + " ('Vegetation-aware', depth_veg)]:\n", + " vals = d.values[~np.isnan(d.values)]\n", + " n_flooded = len(vals)\n", + " pct = 100 * n_flooded / (H * W)\n", + " if len(vals) > 0:\n", + " print(f\"{label:>25s}: {n_flooded:>6d} cells ({pct:.1f}%), \"\n", + " f\"median={np.median(vals):.2f} m, max={vals.max():.2f} m\")\n", + " else:\n", + " print(f\"{label:>25s}: no flooded cells\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fdv-plot", + "metadata": {}, + "outputs": [], + "source": [ + "depth_cmap = LinearSegmentedColormap.from_list('depth', ['#ffffcc', '#fd8d3c', '#800026'])\n", + "\n", + "# Get shared vmax from whichever has deeper water\n", + "d1 = depth_standard.values[~np.isnan(depth_standard.values)]\n", + "d2 = depth_veg.values[~np.isnan(depth_veg.values)]\n", + "depth_vmax = max(np.percentile(d1, 98) if len(d1) else 1,\n", + " np.percentile(d2, 98) if len(d2) else 1)\n", + "\n", + "fig, axes = plt.subplots(1, 2, figsize=(16, 6))\n", + "\n", + "ax = axes[0]\n", + "hillshade.plot.imshow(ax=ax, cmap='gray', add_colorbar=False)\n", + "depth_standard.plot.imshow(ax=ax, cmap=depth_cmap, alpha=0.85, vmin=0, vmax=depth_vmax,\n", + " add_colorbar=True,\n", + " cbar_kwargs={'label': 'Flood depth (m)', 'shrink': 0.7})\n", + "ax.set_title('Standard: uniform water level = 10 m')\n", + "ax.set_axis_off()\n", + "\n", + "ax = axes[1]\n", + "hillshade.plot.imshow(ax=ax, cmap='gray', add_colorbar=False)\n", + "depth_veg.plot.imshow(ax=ax, cmap=depth_cmap, alpha=0.85, vmin=0, vmax=depth_vmax,\n", + " add_colorbar=True,\n", + " cbar_kwargs={'label': 'Flood depth (m)', 'shrink': 0.7})\n", + "ax.set_title('Vegetation-aware: Manning normal depth')\n", + "ax.set_axis_off()\n", + "\n", + "plt.tight_layout()" + ] + }, + { + "cell_type": "markdown", + "id": "fdv-roughness-effect-md", + "metadata": {}, + "source": [ + "### Effect of roughness on flood depth\n", + "\n", + "To isolate the roughness effect, we compare `flood_depth_vegetation` with low uniform n (bare ground) vs. the NLCD-derived n, keeping all other inputs the same." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fdv-roughness-effect", + "metadata": {}, + "outputs": [], + "source": [ + "depth_bare = flood_depth_vegetation(hand_raster, slope_raster,\n", + " mannings_n=0.025, unit_discharge=2.0)\n", + "depth_veg2 = flood_depth_vegetation(hand_raster, slope_raster,\n", + " mannings_n=n_raster, unit_discharge=2.0)\n", + "\n", + "for label, d in [('Bare (n=0.025)', depth_bare),\n", + " ('Vegetation (NLCD n)', depth_veg2)]:\n", + " vals = d.values[~np.isnan(d.values)]\n", + " if len(vals) > 0:\n", + " print(f\"{label:>25s}: {len(vals):>6d} cells flooded, \"\n", + " f\"median={np.median(vals):.3f} m, max={vals.max():.3f} m\")\n", + "\n", + "# Depth difference where both are flooded\n", + "both_flooded = ~np.isnan(depth_bare.values) & ~np.isnan(depth_veg2.values)\n", + "if both_flooded.any():\n", + " diff = depth_veg2.values[both_flooded] - depth_bare.values[both_flooded]\n", + " print(f\"\\nDepth difference (veg - bare) where both flooded:\")\n", + " print(f\" mean: {diff.mean():+.3f} m (positive = vegetation makes it deeper)\")\n", + " print(f\" max: {diff.max():+.3f} m\")" + ] + }, + { + "cell_type": "markdown", + "id": "pipeline-md", + "metadata": {}, + "source": [ + "## Full pipeline comparison\n", + "\n", + "Putting it all together: a side-by-side of the standard workflow (uniform parameters) vs. the vegetation-aware workflow (parameters derived from land cover and soil data).\n", + "\n", + "| Step | Standard | Vegetation-aware |\n", + "|------|----------|------------------|\n", + "| Roughness | n = 0.03 everywhere | `vegetation_roughness(nlcd)` |\n", + "| Curve number | CN = 80 everywhere | `vegetation_curve_number(nlcd, soil)` |\n", + "| Runoff | `curve_number_runoff(rain, 80)` | `curve_number_runoff(rain, cn_raster)` |\n", + "| Travel time | `travel_time(fl, slope, 0.03)` | `travel_time(fl, slope, n_raster)` |\n", + "| Flood depth | `flood_depth(hand, 10)` | `flood_depth_vegetation(hand, slope, n, q)` |" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "pipeline-compute", + "metadata": {}, + "outputs": [], + "source": [ + "# Standard pipeline\n", + "runoff_std = curve_number_runoff(rainfall, curve_number=80.0)\n", + "tt_std = travel_time(fl_raster, slope_raster, mannings_n=0.03)\n", + "depth_std = flood_depth(hand_raster, water_level=10)\n", + "\n", + "# Vegetation-aware pipeline\n", + "n_veg = vegetation_roughness(nlcd_raster, mode='nlcd')\n", + "cn_veg = vegetation_curve_number(nlcd_raster, sg_raster)\n", + "runoff_va = curve_number_runoff(rainfall, curve_number=cn_veg)\n", + "tt_va = travel_time(fl_raster, slope_raster, mannings_n=n_veg)\n", + "depth_va = flood_depth_vegetation(hand_raster, slope_raster, n_veg,\n", + " unit_discharge=2.0)\n", + "\n", + "print(\"Metric Standard Vegetation-aware\")\n", + "print(\"-\" * 60)\n", + "\n", + "ro_std_mean = np.nanmean(runoff_std.values)\n", + "ro_va_mean = np.nanmean(runoff_va.values)\n", + "print(f\"Mean runoff (mm) {ro_std_mean:>8.1f} {ro_va_mean:>8.1f}\")\n", + "\n", + "tc_std = np.nanmax(tt_std.values[np.isfinite(tt_std.values)])\n", + "tc_va = np.nanmax(tt_va.values[np.isfinite(tt_va.values)])\n", + "print(f\"Time of concentration {tc_std:>8.1f} {tc_va:>8.1f}\")\n", + "\n", + "d_std_vals = depth_std.values[~np.isnan(depth_std.values)]\n", + "d_va_vals = depth_va.values[~np.isnan(depth_va.values)]\n", + "print(f\"Flooded cells {len(d_std_vals):>8d} {len(d_va_vals):>8d}\")\n", + "if len(d_std_vals) > 0 and len(d_va_vals) > 0:\n", + " print(f\"Median flood depth (m) {np.median(d_std_vals):>8.2f} {np.median(d_va_vals):>8.2f}\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "pipeline-plot", + "metadata": {}, + "outputs": [], + "source": [ + "fig, axes = plt.subplots(2, 3, figsize=(20, 12))\n", + "\n", + "# Row 1: Standard\n", + "titles_std = ['Standard: runoff', 'Standard: travel time', 'Standard: flood depth']\n", + "data_std = [runoff_std, tt_std, depth_std]\n", + "cmaps = [runoff_cmap, tt_cmap, depth_cmap]\n", + "labels = ['Runoff (mm)', 'Travel time', 'Flood depth (m)']\n", + "\n", + "# Row 2: Vegetation-aware\n", + "titles_va = ['Vegetation: runoff', 'Vegetation: travel time', 'Vegetation: flood depth']\n", + "data_va = [runoff_va, tt_va, depth_va]\n", + "\n", + "for row, (titles, datasets) in enumerate([(titles_std, data_std), (titles_va, data_va)]):\n", + " for col, (title, dat, cmap, label) in enumerate(zip(titles, datasets, cmaps, labels)):\n", + " ax = axes[row, col]\n", + " hillshade.plot.imshow(ax=ax, cmap='gray', add_colorbar=False)\n", + " dat.plot.imshow(ax=ax, cmap=cmap, alpha=0.8, add_colorbar=True,\n", + " cbar_kwargs={'label': label, 'shrink': 0.7})\n", + " ax.set_title(title)\n", + " ax.set_axis_off()\n", + "\n", + "plt.tight_layout()" + ] + }, + { + "cell_type": "markdown", + "id": "takeaways", + "metadata": {}, + "source": [ + "### Key takeaways\n", + "\n", + "- **Spatially varying roughness** from `vegetation_roughness` gives a more realistic travel time pattern than any single uniform value. Forested slopes slow water; developed channels don't.\n", + "- **CN from land cover + soil** via `vegetation_curve_number` replaces manual `xr.where` logic with a single lookup call covering 15 NLCD classes and 4 soil groups. The resulting runoff pattern tracks land use in a way that a 2-class approximation can't.\n", + "- **`flood_depth_vegetation`** uses Manning's normal depth instead of a uniform water level, so flood depth responds to both roughness and slope. Higher roughness produces deeper, slower water. This won't replace a hydraulic model, but it's a step closer to physical realism.\n", + "- All three functions accept custom lookup tables, so you can substitute your own NLCD-to-n or NLCD-to-CN mapping without changing the pipeline." + ] + }, + { + "cell_type": "markdown", + "id": "references", + "metadata": {}, + "source": [ + "### References\n", + "\n", + "- Chow, V.T. (1959). Open Channel Hydraulics. Manning's n tables.\n", + "- Arcement, G.J. & Schneider, V.R. (1989). Guide for selecting Manning's roughness coefficients. USGS Water Supply Paper 2339.\n", + "- USDA-NRCS (1986). Urban Hydrology for Small Watersheds, TR-55. SCS Curve Number tables.\n", + "- Baptist, M.J. et al. (2007). On inducing equations for vegetation resistance. Journal of Hydraulic Research 45(4), 435-450." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "name": "python", + "version": "3.10.0" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} \ No newline at end of file diff --git a/xrspatial/__init__.py b/xrspatial/__init__.py index 3eb227e9..19e66638 100644 --- a/xrspatial/__init__.py +++ b/xrspatial/__init__.py @@ -37,8 +37,11 @@ from xrspatial.fire import rdnbr # noqa from xrspatial.flood import curve_number_runoff # noqa from xrspatial.flood import flood_depth # noqa +from xrspatial.flood import flood_depth_vegetation # noqa from xrspatial.flood import inundation # noqa from xrspatial.flood import travel_time # noqa +from xrspatial.flood import vegetation_curve_number # noqa +from xrspatial.flood import vegetation_roughness # noqa from xrspatial.hydro import flow_accumulation # noqa: unified wrapper from xrspatial.hydro import flow_accumulation_d8, flow_accumulation_dinf, flow_accumulation_mfd # noqa from xrspatial.hydro import flow_direction # noqa: unified wrapper diff --git a/xrspatial/accessor.py b/xrspatial/accessor.py index 82a0ce28..5d582279 100644 --- a/xrspatial/accessor.py +++ b/xrspatial/accessor.py @@ -141,6 +141,20 @@ def travel_time(self, slope_agg, mannings_n, **kwargs): from .flood import travel_time return travel_time(self._obj, slope_agg, mannings_n, **kwargs) + def vegetation_roughness(self, **kwargs): + from .flood import vegetation_roughness + return vegetation_roughness(self._obj, **kwargs) + + def vegetation_curve_number(self, soil_group_agg, **kwargs): + from .flood import vegetation_curve_number + return vegetation_curve_number(self._obj, soil_group_agg, **kwargs) + + def flood_depth_vegetation(self, slope_agg, mannings_n, + unit_discharge, **kwargs): + from .flood import flood_depth_vegetation + return flood_depth_vegetation(self._obj, slope_agg, mannings_n, + unit_discharge, **kwargs) + def viewshed(self, x, y, **kwargs): from .viewshed import viewshed return viewshed(self._obj, x, y, **kwargs) @@ -537,6 +551,20 @@ def travel_time(self, slope_agg, mannings_n, **kwargs): from .flood import travel_time return travel_time(self._obj, slope_agg, mannings_n, **kwargs) + def vegetation_roughness(self, **kwargs): + from .flood import vegetation_roughness + return vegetation_roughness(self._obj, **kwargs) + + def vegetation_curve_number(self, soil_group_agg, **kwargs): + from .flood import vegetation_curve_number + return vegetation_curve_number(self._obj, soil_group_agg, **kwargs) + + def flood_depth_vegetation(self, slope_agg, mannings_n, + unit_discharge, **kwargs): + from .flood import flood_depth_vegetation + return flood_depth_vegetation(self._obj, slope_agg, mannings_n, + unit_discharge, **kwargs) + # ---- Classification ---- def natural_breaks(self, **kwargs): diff --git a/xrspatial/flood.py b/xrspatial/flood.py index be49f920..efc983ab 100644 --- a/xrspatial/flood.py +++ b/xrspatial/flood.py @@ -9,6 +9,11 @@ - ``curve_number_runoff``: SCS/NRCS curve number runoff estimation - ``travel_time``: overland flow travel time via simplified Manning's equation +- ``vegetation_roughness``: Manning's n from land cover or NDVI +- ``vegetation_curve_number``: SCS curve number from land cover and + hydrologic soil group +- ``flood_depth_vegetation``: vegetation-adjusted flood depth via + Manning's normal depth equation """ from __future__ import annotations @@ -33,6 +38,53 @@ # Minimum tan(slope) clamp: tan(0.001 deg), same as TWI _TAN_MIN = np.tan(np.radians(0.001)) +# --------------------------------------------------------------------------- +# NLCD-to-Manning's n lookup (Chow 1959; Arcement & Schneider 1989) +# --------------------------------------------------------------------------- +NLCD_MANNINGS_N = { + 11: 0.025, # Open Water + 21: 0.040, # Developed, Open Space + 22: 0.100, # Developed, Low Intensity + 23: 0.080, # Developed, Medium Intensity + 24: 0.150, # Developed, High Intensity + 31: 0.028, # Barren Land + 41: 0.100, # Deciduous Forest + 42: 0.120, # Evergreen Forest + 43: 0.100, # Mixed Forest + 52: 0.070, # Shrub/Scrub + 71: 0.035, # Grassland/Herbaceous + 81: 0.033, # Pasture/Hay + 82: 0.038, # Cultivated Crops + 90: 0.100, # Woody Wetlands + 95: 0.070, # Emergent Herbaceous Wetlands +} + +# --------------------------------------------------------------------------- +# NLCD + Hydrologic Soil Group -> SCS Curve Number (TR-55, HEC-HMS) +# Soil groups: 1=A, 2=B, 3=C, 4=D. "Good" hydrologic condition assumed. +# --------------------------------------------------------------------------- +NLCD_CURVE_NUMBER = { + (11, 1): 100, (11, 2): 100, (11, 3): 100, (11, 4): 100, # Open Water + (21, 1): 39, (21, 2): 61, (21, 3): 74, (21, 4): 80, # Dev Open + (22, 1): 57, (22, 2): 72, (22, 3): 81, (22, 4): 86, # Dev Low + (23, 1): 77, (23, 2): 85, (23, 3): 90, (23, 4): 92, # Dev Med + (24, 1): 89, (24, 2): 92, (24, 3): 94, (24, 4): 95, # Dev High + (31, 1): 77, (31, 2): 86, (31, 3): 91, (31, 4): 94, # Barren + (41, 1): 30, (41, 2): 55, (41, 3): 70, (41, 4): 77, # Deciduous + (42, 1): 30, (42, 2): 55, (42, 3): 70, (42, 4): 77, # Evergreen + (43, 1): 30, (43, 2): 55, (43, 3): 70, (43, 4): 77, # Mixed Forest + (52, 1): 35, (52, 2): 56, (52, 3): 70, (52, 4): 77, # Shrub + (71, 1): 39, (71, 2): 61, (71, 3): 74, (71, 4): 80, # Grassland + (81, 1): 39, (81, 2): 61, (81, 3): 74, (81, 4): 80, # Pasture + (82, 1): 67, (82, 2): 78, (82, 3): 85, (82, 4): 89, # Crops + (90, 1): 30, (90, 2): 58, (90, 3): 71, (90, 4): 78, # Woody Wetland + (95, 1): 30, (95, 2): 58, (95, 3): 71, (95, 4): 78, # Herb Wetland +} + +# NDVI-to-Manning's n interpolation breakpoints +_NDVI_BREAKPOINTS = np.array([0.0, 0.1, 0.3, 0.6, 1.0]) +_NDVI_MANNINGS_N = np.array([0.02, 0.03, 0.05, 0.10, 0.16]) + # --------------------------------------------------------------------------- # flood_depth @@ -443,3 +495,447 @@ def _travel_time_dask_cupy(fl, sl, n): cp.asarray, dtype=result.dtype, meta=cp.array((), dtype=result.dtype), ) + + +# --------------------------------------------------------------------------- +# vegetation_roughness +# --------------------------------------------------------------------------- + +def vegetation_roughness( + vegetation_agg: xr.DataArray, + mode: str = 'nlcd', + lookup: dict | None = None, + name: str = 'mannings_n', +) -> xr.DataArray: + """Derive Manning's n roughness from land cover or NDVI. + + Parameters + ---------- + vegetation_agg : xarray.DataArray + 2D raster. For ``mode='nlcd'`` this should contain integer NLCD + land cover codes. For ``mode='ndvi'`` this should contain + continuous NDVI values (typically -1 to 1). + mode : {'nlcd', 'ndvi'}, default 'nlcd' + How to convert vegetation data to roughness. + + * ``'nlcd'``: categorical lookup from NLCD codes to Manning's n. + * ``'ndvi'``: piecewise linear interpolation of NDVI to Manning's n. + lookup : dict, optional + Override the default NLCD-to-n mapping. Keys are integer NLCD + codes, values are Manning's n floats. Only used when + ``mode='nlcd'``. + name : str, default 'mannings_n' + Name of the output DataArray. + + Returns + ------- + xarray.DataArray + 2D float64 Manning's n raster. Unrecognized NLCD codes and + NaN inputs map to NaN. + """ + _validate_raster(vegetation_agg, func_name='vegetation_roughness', + name='vegetation_agg') + if mode not in ('nlcd', 'ndvi'): + raise ValueError(f"mode must be 'nlcd' or 'ndvi', got {mode!r}") + + data = vegetation_agg.data + + if mode == 'nlcd': + lut_dict = dict(NLCD_MANNINGS_N) + if lookup is not None: + lut_dict.update(lookup) + max_code = max(lut_dict) + lut = np.full(max_code + 1, np.nan, dtype=np.float64) + for code, val in lut_dict.items(): + lut[code] = val + + if has_cuda_and_cupy() and is_dask_cupy(vegetation_agg): + out = _veg_roughness_nlcd_dask_cupy(data, lut) + elif has_cuda_and_cupy() and is_cupy_array(data): + out = _veg_roughness_nlcd_cupy(data, lut) + elif da is not None and isinstance(data, da.Array): + out = _veg_roughness_nlcd_dask(data, lut) + elif isinstance(data, np.ndarray): + out = _veg_roughness_nlcd_numpy(data, lut) + else: + raise TypeError(f"Unsupported array type: {type(data)}") + else: + if has_cuda_and_cupy() and is_dask_cupy(vegetation_agg): + out = _veg_roughness_ndvi_dask_cupy(data) + elif has_cuda_and_cupy() and is_cupy_array(data): + out = _veg_roughness_ndvi_cupy(data) + elif da is not None and isinstance(data, da.Array): + out = _veg_roughness_ndvi_dask(data) + elif isinstance(data, np.ndarray): + out = _veg_roughness_ndvi_numpy(data) + else: + raise TypeError(f"Unsupported array type: {type(data)}") + + return xr.DataArray(out, name=name, coords=vegetation_agg.coords, + dims=vegetation_agg.dims, + attrs=vegetation_agg.attrs) + + +def _veg_roughness_nlcd_numpy(data, lut): + data = np.asarray(data) + out = np.full(data.shape, np.nan, dtype=np.float64) + not_nan = ~np.isnan(data) if data.dtype.kind == 'f' else True + mask = (data >= 0) & (data < len(lut)) & not_nan + idx = data[mask].astype(np.intp) + out[mask] = lut[idx] + return out + + +def _veg_roughness_nlcd_cupy(data, lut): + import cupy as cp + lut_gpu = cp.asarray(lut) + data = cp.asarray(data) + out = cp.full(data.shape, cp.nan, dtype=cp.float64) + not_nan = ~cp.isnan(data) if data.dtype.kind == 'f' else True + mask = (data >= 0) & (data < len(lut_gpu)) & not_nan + idx = data[mask].astype(cp.intp) + out[mask] = lut_gpu[idx] + return out + + +def _veg_roughness_nlcd_dask(data, lut): + import dask.array as _da + + def _apply_lut(block, lut=lut): + return _veg_roughness_nlcd_numpy(block, lut) + + return _da.map_blocks( + _apply_lut, data, dtype=np.float64, + meta=np.array((), dtype=np.float64), + ) + + +def _veg_roughness_nlcd_dask_cupy(data, lut): + import cupy as cp + data_np = data.map_blocks( + lambda b: b.get(), dtype=data.dtype, + meta=np.array((), dtype=data.dtype), + ) + result = _veg_roughness_nlcd_dask(data_np, lut) + return result.map_blocks( + cp.asarray, dtype=result.dtype, + meta=cp.array((), dtype=result.dtype), + ) + + +def _veg_roughness_ndvi_numpy(data): + data = np.asarray(data, dtype=np.float64) + out = np.interp(data, _NDVI_BREAKPOINTS, _NDVI_MANNINGS_N) + out = np.where(np.isnan(data), np.nan, out) + return out + + +def _veg_roughness_ndvi_cupy(data): + import cupy as cp + data = cp.asarray(data, dtype=cp.float64) + bp = cp.asarray(_NDVI_BREAKPOINTS) + mn = cp.asarray(_NDVI_MANNINGS_N) + out = cp.interp(data.ravel(), bp, mn).reshape(data.shape) + out = cp.where(cp.isnan(data), cp.nan, out) + return out + + +def _veg_roughness_ndvi_dask(data): + import dask.array as _da + + def _apply(block): + return _veg_roughness_ndvi_numpy(block) + + return _da.map_blocks( + _apply, data, dtype=np.float64, + meta=np.array((), dtype=np.float64), + ) + + +def _veg_roughness_ndvi_dask_cupy(data): + import cupy as cp + data_np = data.map_blocks( + lambda b: b.get(), dtype=data.dtype, + meta=np.array((), dtype=data.dtype), + ) + result = _veg_roughness_ndvi_dask(data_np) + return result.map_blocks( + cp.asarray, dtype=result.dtype, + meta=cp.array((), dtype=result.dtype), + ) + + +# --------------------------------------------------------------------------- +# vegetation_curve_number +# --------------------------------------------------------------------------- + +def vegetation_curve_number( + landcover_agg: xr.DataArray, + soil_group_agg: xr.DataArray, + lookup: dict | None = None, + name: str = 'curve_number', +) -> xr.DataArray: + """Derive SCS curve number from NLCD land cover and hydrologic soil group. + + Parameters + ---------- + landcover_agg : xarray.DataArray + 2D integer raster of NLCD land cover codes. + soil_group_agg : xarray.DataArray + 2D integer raster of hydrologic soil groups (1=A, 2=B, 3=C, 4=D). + lookup : dict, optional + Override the default ``(nlcd_code, soil_group) -> CN`` mapping. + name : str, default 'curve_number' + Name of the output DataArray. + + Returns + ------- + xarray.DataArray + 2D float64 curve number raster. Unknown ``(code, group)`` + pairs and NaN inputs map to NaN. + """ + _validate_raster(landcover_agg, func_name='vegetation_curve_number', + name='landcover_agg') + _validate_raster(soil_group_agg, func_name='vegetation_curve_number', + name='soil_group_agg') + + lut_dict = dict(NLCD_CURVE_NUMBER) + if lookup is not None: + lut_dict.update(lookup) + + max_code = max(k[0] for k in lut_dict) + max_sg = max(k[1] for k in lut_dict) + lut = np.full((max_code + 1, max_sg + 1), np.nan, dtype=np.float64) + for (code, sg), cn in lut_dict.items(): + lut[code, sg] = cn + + lc_data = landcover_agg.data + sg_data = soil_group_agg.data + + if has_cuda_and_cupy() and is_dask_cupy(landcover_agg): + out = _veg_cn_dask_cupy(lc_data, sg_data, lut) + elif has_cuda_and_cupy() and is_cupy_array(lc_data): + out = _veg_cn_cupy(lc_data, sg_data, lut) + elif da is not None and isinstance(lc_data, da.Array): + out = _veg_cn_dask(lc_data, sg_data, lut) + elif isinstance(lc_data, np.ndarray): + out = _veg_cn_numpy(lc_data, sg_data, lut) + else: + raise TypeError(f"Unsupported array type: {type(lc_data)}") + + return xr.DataArray(out, name=name, coords=landcover_agg.coords, + dims=landcover_agg.dims, + attrs=landcover_agg.attrs) + + +def _veg_cn_numpy(lc, sg, lut): + lc = np.asarray(lc) + sg = np.asarray(sg) + max_code, max_sg = lut.shape + out = np.full(lc.shape, np.nan, dtype=np.float64) + lc_ok = ~np.isnan(lc) if lc.dtype.kind == 'f' else True + sg_ok = ~np.isnan(sg) if sg.dtype.kind == 'f' else True + mask = ( + (lc >= 0) & (lc < max_code) + & (sg >= 0) & (sg < max_sg) + & lc_ok & sg_ok + ) + lc_idx = lc[mask].astype(np.intp) + sg_idx = sg[mask].astype(np.intp) + out[mask] = lut[lc_idx, sg_idx] + return out + + +def _veg_cn_cupy(lc, sg, lut): + import cupy as cp + lut_gpu = cp.asarray(lut) + lc = cp.asarray(lc) + sg = cp.asarray(sg) + max_code, max_sg = lut_gpu.shape + out = cp.full(lc.shape, cp.nan, dtype=cp.float64) + lc_ok = ~cp.isnan(lc) if lc.dtype.kind == 'f' else True + sg_ok = ~cp.isnan(sg) if sg.dtype.kind == 'f' else True + mask = ( + (lc >= 0) & (lc < max_code) + & (sg >= 0) & (sg < max_sg) + & lc_ok & sg_ok + ) + lc_idx = lc[mask].astype(cp.intp) + sg_idx = sg[mask].astype(cp.intp) + out[mask] = lut_gpu[lc_idx, sg_idx] + return out + + +def _veg_cn_dask(lc, sg, lut): + import dask.array as _da + + def _apply(lc_block, sg_block, lut=lut): + return _veg_cn_numpy(lc_block, sg_block, lut) + + return _da.map_blocks( + _apply, lc, sg, dtype=np.float64, + meta=np.array((), dtype=np.float64), + ) + + +def _veg_cn_dask_cupy(lc, sg, lut): + import cupy as cp + lc_np = lc.map_blocks( + lambda b: b.get(), dtype=lc.dtype, + meta=np.array((), dtype=lc.dtype), + ) + sg_np = sg.map_blocks( + lambda b: b.get(), dtype=sg.dtype, + meta=np.array((), dtype=sg.dtype), + ) + result = _veg_cn_dask(lc_np, sg_np, lut) + return result.map_blocks( + cp.asarray, dtype=result.dtype, + meta=cp.array((), dtype=result.dtype), + ) + + +# --------------------------------------------------------------------------- +# flood_depth_vegetation +# --------------------------------------------------------------------------- + +def flood_depth_vegetation( + hand_agg: xr.DataArray, + slope_agg: xr.DataArray, + mannings_n: Union[float, xr.DataArray], + unit_discharge: float, + name: str = 'flood_depth_vegetation', +) -> xr.DataArray: + """Compute vegetation-adjusted flood depth via Manning's normal depth. + + For each cell, the equilibrium (normal) depth for steady uniform flow + is:: + + h = (q * n / sqrt(tan(slope)))^(3/5) + + where *q* is unit discharge (m^2/s), *n* is Manning's roughness, and + *slope* is terrain slope in degrees. The flood depth is + ``h - HAND`` where ``h > HAND``, and NaN elsewhere. + + Higher roughness (denser vegetation) produces deeper, slower water. + + Parameters + ---------- + hand_agg : xarray.DataArray + 2D HAND raster (output of ``hand()``). + slope_agg : xarray.DataArray + 2D slope raster in degrees. + mannings_n : float or xarray.DataArray + Manning's roughness coefficient. A scalar applies uniformly; + a DataArray (e.g. output of ``vegetation_roughness()``) allows + spatially varying roughness. + unit_discharge : float + Unit discharge in m^2/s (discharge per unit width, > 0). + name : str, default 'flood_depth_vegetation' + Name of the output DataArray. + + Returns + ------- + xarray.DataArray + 2D float64 flood depth grid. NaN where not inundated or + where inputs are NaN. + """ + _validate_raster(hand_agg, func_name='flood_depth_vegetation', + name='hand_agg') + _validate_raster(slope_agg, func_name='flood_depth_vegetation', + name='slope_agg') + if not isinstance(unit_discharge, (int, float)): + raise TypeError("unit_discharge must be numeric") + if unit_discharge <= 0: + raise ValueError("unit_discharge must be > 0") + + if isinstance(mannings_n, xr.DataArray): + n_data = mannings_n.data + elif isinstance(mannings_n, (int, float)): + if mannings_n <= 0: + raise ValueError("mannings_n must be > 0") + n_data = float(mannings_n) + else: + raise TypeError( + "mannings_n must be numeric or an xarray.DataArray") + + h_data = hand_agg.data + sl_data = slope_agg.data + q = float(unit_discharge) + + if has_cuda_and_cupy() and is_dask_cupy(hand_agg): + out = _fdv_dask_cupy(h_data, sl_data, n_data, q) + elif has_cuda_and_cupy() and is_cupy_array(h_data): + out = _fdv_cupy(h_data, sl_data, n_data, q) + elif da is not None and isinstance(h_data, da.Array): + out = _fdv_dask(h_data, sl_data, n_data, q) + elif isinstance(h_data, np.ndarray): + out = _fdv_numpy(h_data, sl_data, n_data, q) + else: + raise TypeError(f"Unsupported array type: {type(h_data)}") + + return xr.DataArray(out, name=name, coords=hand_agg.coords, + dims=hand_agg.dims, attrs=hand_agg.attrs) + + +def _fdv_numpy(hand, sl, n, q): + hand = np.asarray(hand, dtype=np.float64) + sl = np.asarray(sl, dtype=np.float64) + n = np.asarray(n, dtype=np.float64) + tan_s = np.tan(np.radians(sl)) + tan_s = np.where(tan_s < _TAN_MIN, _TAN_MIN, tan_s) + h_normal = (q * n / np.sqrt(tan_s)) ** 0.6 + depth = np.where(h_normal > hand, h_normal - hand, np.nan) + # propagate NaN from any input + any_nan = np.isnan(hand) | np.isnan(sl) | np.isnan(n) + depth = np.where(any_nan, np.nan, depth) + return depth + + +def _fdv_cupy(hand, sl, n, q): + import cupy as cp + hand = cp.asarray(hand, dtype=cp.float64) + sl = cp.asarray(sl, dtype=cp.float64) + n = cp.asarray(n, dtype=cp.float64) + tan_s = cp.tan(cp.radians(sl)) + tan_s = cp.where(tan_s < _TAN_MIN, _TAN_MIN, tan_s) + h_normal = (q * n / cp.sqrt(tan_s)) ** 0.6 + depth = cp.where(h_normal > hand, h_normal - hand, cp.nan) + any_nan = cp.isnan(hand) | cp.isnan(sl) | cp.isnan(n) + depth = cp.where(any_nan, cp.nan, depth) + return depth + + +def _fdv_dask(hand, sl, n, q): + import dask.array as _da + tan_s = _da.tan(_da.radians(sl)) + tan_s = _da.where(tan_s < _TAN_MIN, _TAN_MIN, tan_s) + h_normal = (q * n / _da.sqrt(tan_s)) ** 0.6 + depth = _da.where(h_normal > hand, h_normal - hand, np.nan) + any_nan = _da.isnan(hand) | _da.isnan(sl) | _da.isnan(n) + depth = _da.where(any_nan, np.nan, depth) + return depth + + +def _fdv_dask_cupy(hand, sl, n, q): + import cupy as cp + hand_np = hand.map_blocks( + lambda b: b.get(), dtype=hand.dtype, + meta=np.array((), dtype=hand.dtype), + ) + sl_np = sl.map_blocks( + lambda b: b.get(), dtype=sl.dtype, + meta=np.array((), dtype=sl.dtype), + ) + if hasattr(n, 'map_blocks'): + n_np = n.map_blocks( + lambda b: b.get(), dtype=n.dtype, + meta=np.array((), dtype=n.dtype), + ) + else: + n_np = n + result = _fdv_dask(hand_np, sl_np, n_np, q) + return result.map_blocks( + cp.asarray, dtype=result.dtype, + meta=cp.array((), dtype=result.dtype), + ) diff --git a/xrspatial/tests/test_flood.py b/xrspatial/tests/test_flood.py index 27ebcbc4..e3c7304e 100644 --- a/xrspatial/tests/test_flood.py +++ b/xrspatial/tests/test_flood.py @@ -5,10 +5,15 @@ import xarray as xr from xrspatial.flood import ( + NLCD_CURVE_NUMBER, + NLCD_MANNINGS_N, curve_number_runoff, flood_depth, + flood_depth_vegetation, inundation, travel_time, + vegetation_curve_number, + vegetation_roughness, ) from xrspatial.tests.general_checks import ( create_test_raster, @@ -508,3 +513,445 @@ def test_numpy_equals_dask_cupy(self): general_output_checks(fl_dc, result_dc, expected_results=result_np.data) + + +# =================================================================== +# vegetation_roughness +# =================================================================== + +class TestVegRoughnessNLCDKnownValues: + + def test_known_codes(self): + # Deciduous Forest=41 -> 0.100, Grassland=71 -> 0.035 + data = np.array([[41, 71]], dtype=np.int32) + raster = create_test_raster(data, name='nlcd') + result = vegetation_roughness(raster, mode='nlcd') + np.testing.assert_allclose(result.data[0, 0], 0.100) + np.testing.assert_allclose(result.data[0, 1], 0.035) + + def test_all_nlcd_codes(self): + codes = sorted(NLCD_MANNINGS_N.keys()) + data = np.array([codes], dtype=np.int32) + raster = create_test_raster(data, name='nlcd') + result = vegetation_roughness(raster, mode='nlcd') + for i, code in enumerate(codes): + np.testing.assert_allclose( + result.data[0, i], NLCD_MANNINGS_N[code], + err_msg=f"NLCD code {code}") + + def test_unrecognized_code_nan(self): + data = np.array([[41, 99]], dtype=np.int32) + raster = create_test_raster(data, name='nlcd') + result = vegetation_roughness(raster, mode='nlcd') + np.testing.assert_allclose(result.data[0, 0], 0.100) + assert np.isnan(result.data[0, 1]) + + def test_custom_lookup(self): + data = np.array([[41]], dtype=np.int32) + raster = create_test_raster(data, name='nlcd') + result = vegetation_roughness(raster, mode='nlcd', + lookup={41: 0.200}) + np.testing.assert_allclose(result.data[0, 0], 0.200) + + +class TestVegRoughnessNDVIKnownValues: + + def test_breakpoints(self): + # Test at exact breakpoints + data = np.array([[0.0, 0.1, 0.3, 0.6, 1.0]], dtype=np.float64) + raster = create_test_raster(data, name='ndvi') + result = vegetation_roughness(raster, mode='ndvi') + np.testing.assert_allclose( + result.data[0], [0.02, 0.03, 0.05, 0.10, 0.16]) + + def test_interpolation(self): + # Midpoint between 0.1 (n=0.03) and 0.3 (n=0.05) -> n=0.04 + data = np.array([[0.2]], dtype=np.float64) + raster = create_test_raster(data, name='ndvi') + result = vegetation_roughness(raster, mode='ndvi') + np.testing.assert_allclose(result.data[0, 0], 0.04) + + def test_nan_propagation(self): + data = np.array([[np.nan, 0.5]], dtype=np.float64) + raster = create_test_raster(data, name='ndvi') + result = vegetation_roughness(raster, mode='ndvi') + assert np.isnan(result.data[0, 0]) + assert not np.isnan(result.data[0, 1]) + + +class TestVegRoughnessValidation: + + def test_invalid_mode(self): + raster = create_test_raster(np.ones((2, 2), dtype=np.float64)) + with pytest.raises(ValueError, match="mode must be"): + vegetation_roughness(raster, mode='bad') + + +class TestVegRoughnessOutput: + + def test_coords_preserved(self): + data = np.array([[41, 42], [71, 82]], dtype=np.int32) + raster = create_test_raster(data, name='nlcd') + result = vegetation_roughness(raster, mode='nlcd') + general_output_checks(raster, result, verify_attrs=True) + + +@dask_array_available +class TestVegRoughnessDask: + + def test_nlcd_numpy_equals_dask(self): + data = np.array([[41, 71, 82], + [11, 42, 52], + [24, 31, 95]], dtype=np.int32) + r_np = create_test_raster(data, backend='numpy', name='nlcd') + r_da = create_test_raster(data, backend='dask', name='nlcd', + chunks=(2, 2)) + result_np = vegetation_roughness(r_np, mode='nlcd') + result_da = vegetation_roughness(r_da, mode='nlcd') + general_output_checks(r_da, result_da, + expected_results=result_np.data) + + def test_ndvi_numpy_equals_dask(self): + data = np.array([[0.0, 0.2, 0.5], + [0.1, np.nan, 0.8], + [0.3, 0.6, 1.0]], dtype=np.float64) + r_np = create_test_raster(data, backend='numpy', name='ndvi') + r_da = create_test_raster(data, backend='dask', name='ndvi', + chunks=(2, 2)) + result_np = vegetation_roughness(r_np, mode='ndvi') + result_da = vegetation_roughness(r_da, mode='ndvi') + general_output_checks(r_da, result_da, + expected_results=result_np.data) + + +@cuda_and_cupy_available +class TestVegRoughnessCuPy: + + def test_nlcd_numpy_equals_cupy(self): + data = np.array([[41, 71], [11, 82]], dtype=np.int32) + r_np = create_test_raster(data, backend='numpy', name='nlcd') + r_cu = create_test_raster(data, backend='cupy', name='nlcd') + result_np = vegetation_roughness(r_np, mode='nlcd') + result_cu = vegetation_roughness(r_cu, mode='nlcd') + general_output_checks(r_cu, result_cu, + expected_results=result_np.data) + + def test_ndvi_numpy_equals_cupy(self): + data = np.array([[0.0, 0.5], [np.nan, 1.0]], dtype=np.float64) + r_np = create_test_raster(data, backend='numpy', name='ndvi') + r_cu = create_test_raster(data, backend='cupy', name='ndvi') + result_np = vegetation_roughness(r_np, mode='ndvi') + result_cu = vegetation_roughness(r_cu, mode='ndvi') + general_output_checks(r_cu, result_cu, + expected_results=result_np.data) + + +@cuda_and_cupy_available +@dask_array_available +class TestVegRoughnessDaskCuPy: + + def test_nlcd_numpy_equals_dask_cupy(self): + data = np.array([[41, 71], [11, 82]], dtype=np.int32) + r_np = create_test_raster(data, backend='numpy', name='nlcd') + r_dc = create_test_raster(data, backend='dask+cupy', name='nlcd', + chunks=(2, 2)) + result_np = vegetation_roughness(r_np, mode='nlcd') + result_dc = vegetation_roughness(r_dc, mode='nlcd') + general_output_checks(r_dc, result_dc, + expected_results=result_np.data) + + +# =================================================================== +# vegetation_curve_number +# =================================================================== + +class TestVegCNKnownValues: + + def test_basic(self): + # Deciduous Forest (41) on soil group A (1) -> CN=30 + # Developed High (24) on soil group D (4) -> CN=95 + lc = np.array([[41, 24]], dtype=np.int32) + sg = np.array([[1, 4]], dtype=np.int32) + lc_r = create_test_raster(lc, name='nlcd') + sg_r = create_test_raster(sg, name='soil') + result = vegetation_curve_number(lc_r, sg_r) + np.testing.assert_allclose(result.data[0, 0], 30.0) + np.testing.assert_allclose(result.data[0, 1], 95.0) + + def test_all_entries(self): + # Verify every entry in the default lookup table + for (code, sg_val), expected_cn in NLCD_CURVE_NUMBER.items(): + lc = np.array([[code]], dtype=np.int32) + sg = np.array([[sg_val]], dtype=np.int32) + lc_r = create_test_raster(lc, name='nlcd') + sg_r = create_test_raster(sg, name='soil') + result = vegetation_curve_number(lc_r, sg_r) + np.testing.assert_allclose( + result.data[0, 0], expected_cn, + err_msg=f"NLCD={code}, SG={sg_val}") + + def test_unknown_pair_nan(self): + lc = np.array([[99]], dtype=np.int32) + sg = np.array([[1]], dtype=np.int32) + lc_r = create_test_raster(lc, name='nlcd') + sg_r = create_test_raster(sg, name='soil') + result = vegetation_curve_number(lc_r, sg_r) + assert np.isnan(result.data[0, 0]) + + def test_custom_lookup(self): + lc = np.array([[41]], dtype=np.int32) + sg = np.array([[1]], dtype=np.int32) + lc_r = create_test_raster(lc, name='nlcd') + sg_r = create_test_raster(sg, name='soil') + result = vegetation_curve_number(lc_r, sg_r, + lookup={(41, 1): 50}) + np.testing.assert_allclose(result.data[0, 0], 50.0) + + +class TestVegCNOutput: + + def test_coords_preserved(self): + lc = np.array([[41, 42], [71, 82]], dtype=np.int32) + sg = np.array([[1, 2], [3, 4]], dtype=np.int32) + lc_r = create_test_raster(lc, name='nlcd') + sg_r = create_test_raster(sg, name='soil') + result = vegetation_curve_number(lc_r, sg_r) + general_output_checks(lc_r, result, verify_attrs=True) + + +@dask_array_available +class TestVegCNDask: + + def test_numpy_equals_dask(self): + lc = np.array([[41, 71, 82], + [11, 42, 52], + [24, 31, 95]], dtype=np.int32) + sg = np.array([[1, 2, 3], + [4, 1, 2], + [3, 4, 1]], dtype=np.int32) + lc_np = create_test_raster(lc, backend='numpy', name='nlcd') + sg_np = create_test_raster(sg, backend='numpy', name='soil') + lc_da = create_test_raster(lc, backend='dask', name='nlcd', + chunks=(2, 2)) + sg_da = create_test_raster(sg, backend='dask', name='soil', + chunks=(2, 2)) + result_np = vegetation_curve_number(lc_np, sg_np) + result_da = vegetation_curve_number(lc_da, sg_da) + general_output_checks(lc_da, result_da, + expected_results=result_np.data) + + +@cuda_and_cupy_available +class TestVegCNCuPy: + + def test_numpy_equals_cupy(self): + lc = np.array([[41, 71], [24, 82]], dtype=np.int32) + sg = np.array([[1, 2], [4, 3]], dtype=np.int32) + lc_np = create_test_raster(lc, backend='numpy', name='nlcd') + sg_np = create_test_raster(sg, backend='numpy', name='soil') + lc_cu = create_test_raster(lc, backend='cupy', name='nlcd') + sg_cu = create_test_raster(sg, backend='cupy', name='soil') + result_np = vegetation_curve_number(lc_np, sg_np) + result_cu = vegetation_curve_number(lc_cu, sg_cu) + general_output_checks(lc_cu, result_cu, + expected_results=result_np.data) + + +@cuda_and_cupy_available +@dask_array_available +class TestVegCNDaskCuPy: + + def test_numpy_equals_dask_cupy(self): + lc = np.array([[41, 71], [24, 82]], dtype=np.int32) + sg = np.array([[1, 2], [4, 3]], dtype=np.int32) + lc_np = create_test_raster(lc, backend='numpy', name='nlcd') + sg_np = create_test_raster(sg, backend='numpy', name='soil') + lc_dc = create_test_raster(lc, backend='dask+cupy', name='nlcd', + chunks=(2, 2)) + sg_dc = create_test_raster(sg, backend='dask+cupy', name='soil', + chunks=(2, 2)) + result_np = vegetation_curve_number(lc_np, sg_np) + result_dc = vegetation_curve_number(lc_dc, sg_dc) + general_output_checks(lc_dc, result_dc, + expected_results=result_np.data) + + +# =================================================================== +# flood_depth_vegetation +# =================================================================== + +class TestFDVKnownValues: + + def test_basic(self): + # hand=1, slope=45 deg, n=0.10, q=0.5 + # tan(45)=1.0, h_normal = (0.5*0.10/1.0)^0.6 = 0.05^0.6 + # depth = h_normal - 1.0 (if h_normal > 1.0, else NaN) + h_normal = (0.5 * 0.10 / np.sqrt(1.0)) ** 0.6 + hand_data = np.array([[0.0, 1.0]], dtype=np.float64) + sl_data = np.array([[45.0, 45.0]], dtype=np.float64) + hand = create_test_raster(hand_data, name='hand') + sl = create_test_raster(sl_data, name='slope') + result = flood_depth_vegetation(hand, sl, mannings_n=0.10, + unit_discharge=0.5) + # HAND=0 cell: depth = h_normal - 0 = h_normal + np.testing.assert_allclose(result.data[0, 0], h_normal, rtol=1e-10) + + def test_not_inundated(self): + # Large HAND means no inundation -> NaN + hand_data = np.array([[100.0]], dtype=np.float64) + sl_data = np.array([[10.0]], dtype=np.float64) + hand = create_test_raster(hand_data, name='hand') + sl = create_test_raster(sl_data, name='slope') + result = flood_depth_vegetation(hand, sl, mannings_n=0.05, + unit_discharge=0.5) + assert np.isnan(result.data[0, 0]) + + def test_higher_roughness_deeper(self): + # Physical invariant: higher n -> deeper water + hand_data = np.array([[0.0, 0.0]], dtype=np.float64) + sl_data = np.array([[10.0, 10.0]], dtype=np.float64) + n_data = np.array([[0.03, 0.15]], dtype=np.float64) + hand = create_test_raster(hand_data, name='hand') + sl = create_test_raster(sl_data, name='slope') + n_raster = create_test_raster(n_data, name='n') + result = flood_depth_vegetation(hand, sl, mannings_n=n_raster, + unit_discharge=1.0) + # n=0.15 cell should have deeper water than n=0.03 cell + assert result.data[0, 1] > result.data[0, 0] + + def test_zero_slope_no_inf(self): + hand_data = np.array([[0.0]], dtype=np.float64) + sl_data = np.array([[0.0]], dtype=np.float64) + hand = create_test_raster(hand_data, name='hand') + sl = create_test_raster(sl_data, name='slope') + result = flood_depth_vegetation(hand, sl, mannings_n=0.05, + unit_discharge=0.5) + assert not np.any(np.isinf(result.data)) + + +class TestFDVNaN: + + def test_nan_propagation(self): + hand_data = np.array([[np.nan, 0.0]], dtype=np.float64) + sl_data = np.array([[10.0, np.nan]], dtype=np.float64) + hand = create_test_raster(hand_data, name='hand') + sl = create_test_raster(sl_data, name='slope') + result = flood_depth_vegetation(hand, sl, mannings_n=0.05, + unit_discharge=0.5) + assert np.isnan(result.data[0, 0]) + assert np.isnan(result.data[0, 1]) + + +class TestFDVValidation: + + def test_negative_discharge(self): + hand = create_test_raster(np.ones((2, 2), dtype=np.float64)) + sl = create_test_raster(np.ones((2, 2), dtype=np.float64)) + with pytest.raises(ValueError, match="unit_discharge must be > 0"): + flood_depth_vegetation(hand, sl, mannings_n=0.05, + unit_discharge=-1.0) + + def test_zero_discharge(self): + hand = create_test_raster(np.ones((2, 2), dtype=np.float64)) + sl = create_test_raster(np.ones((2, 2), dtype=np.float64)) + with pytest.raises(ValueError, match="unit_discharge must be > 0"): + flood_depth_vegetation(hand, sl, mannings_n=0.05, + unit_discharge=0.0) + + def test_non_numeric_discharge(self): + hand = create_test_raster(np.ones((2, 2), dtype=np.float64)) + sl = create_test_raster(np.ones((2, 2), dtype=np.float64)) + with pytest.raises(TypeError, match="unit_discharge must be numeric"): + flood_depth_vegetation(hand, sl, mannings_n=0.05, + unit_discharge="high") + + def test_negative_mannings_n(self): + hand = create_test_raster(np.ones((2, 2), dtype=np.float64)) + sl = create_test_raster(np.ones((2, 2), dtype=np.float64)) + with pytest.raises(ValueError, match="mannings_n must be > 0"): + flood_depth_vegetation(hand, sl, mannings_n=-0.01, + unit_discharge=0.5) + + +class TestFDVOutput: + + def test_coords_preserved(self): + hand_data = np.zeros((4, 4), dtype=np.float64) + sl_data = np.ones((4, 4), dtype=np.float64) * 10.0 + hand = create_test_raster(hand_data, name='hand') + sl = create_test_raster(sl_data, name='slope') + result = flood_depth_vegetation(hand, sl, mannings_n=0.05, + unit_discharge=0.5) + general_output_checks(hand, result, verify_attrs=True) + + +@dask_array_available +class TestFDVDask: + + @pytest.mark.parametrize("chunks", [(2, 2), (3, 3)]) + def test_numpy_equals_dask(self, chunks): + np.random.seed(42) + hand_data = np.random.uniform(0, 5, (5, 5)).astype(np.float64) + sl_data = np.random.uniform(0.1, 60, (5, 5)).astype(np.float64) + + h_np = create_test_raster(hand_data, backend='numpy', name='hand') + s_np = create_test_raster(sl_data, backend='numpy', name='slope') + h_da = create_test_raster(hand_data, backend='dask', name='hand', + chunks=chunks) + s_da = create_test_raster(sl_data, backend='dask', name='slope', + chunks=chunks) + + result_np = flood_depth_vegetation(h_np, s_np, mannings_n=0.10, + unit_discharge=1.0) + result_da = flood_depth_vegetation(h_da, s_da, mannings_n=0.10, + unit_discharge=1.0) + + general_output_checks(h_da, result_da, + expected_results=result_np.data) + + +@cuda_and_cupy_available +class TestFDVCuPy: + + def test_numpy_equals_cupy(self): + np.random.seed(42) + hand_data = np.random.uniform(0, 5, (5, 5)).astype(np.float64) + sl_data = np.random.uniform(0.1, 60, (5, 5)).astype(np.float64) + + h_np = create_test_raster(hand_data, backend='numpy', name='hand') + s_np = create_test_raster(sl_data, backend='numpy', name='slope') + h_cu = create_test_raster(hand_data, backend='cupy', name='hand') + s_cu = create_test_raster(sl_data, backend='cupy', name='slope') + + result_np = flood_depth_vegetation(h_np, s_np, mannings_n=0.10, + unit_discharge=1.0) + result_cu = flood_depth_vegetation(h_cu, s_cu, mannings_n=0.10, + unit_discharge=1.0) + + general_output_checks(h_cu, result_cu, + expected_results=result_np.data) + + +@cuda_and_cupy_available +@dask_array_available +class TestFDVDaskCuPy: + + def test_numpy_equals_dask_cupy(self): + np.random.seed(42) + hand_data = np.random.uniform(0, 5, (5, 5)).astype(np.float64) + sl_data = np.random.uniform(0.1, 60, (5, 5)).astype(np.float64) + + h_np = create_test_raster(hand_data, backend='numpy', name='hand') + s_np = create_test_raster(sl_data, backend='numpy', name='slope') + h_dc = create_test_raster(hand_data, backend='dask+cupy', + name='hand', chunks=(3, 3)) + s_dc = create_test_raster(sl_data, backend='dask+cupy', + name='slope', chunks=(3, 3)) + + result_np = flood_depth_vegetation(h_np, s_np, mannings_n=0.10, + unit_discharge=1.0) + result_dc = flood_depth_vegetation(h_dc, s_dc, mannings_n=0.10, + unit_discharge=1.0) + + general_output_checks(h_dc, result_dc, + expected_results=result_np.data)