Stop NaN input from reading out of bounds and flooding the erosion output (#3703)

.claude/sweep-accuracy-state.csv

@@ -10,6 +10,7 @@ dasymetric,2026-06-20,3403,MEDIUM,2;5,"Cat2/Cat5: disaggregate(limiting_variable
 diffusion,2026-05-01,,LOW,1;2;5,"LOW: no Kahan summation across long iterations (drift over 100k steps, standard for explicit Euler); lap=n+s+w+e-4*val has catastrophic cancellation for nearly-uniform large values; res=0 in attrs causes div-by-zero (no guard); dask+cupy boundary='nan' relies on dask accepting cp.nan as fill. CPU/GPU NaN handling consistent (np.isnan vs val!=val). depth=1 matches stencil radius. Memory guards, CFL check, step cap all in place. No CRIT/HIGH."
 edge_detection,2026-05-01,,,,Thin wrappers around convolve_2d with fixed Sobel/Prewitt/Laplacian kernels; no issues found
 emerging_hotspots,2026-04-30,,MEDIUM,2;3,MEDIUM: threshold_90 uses int() (truncation) instead of ceil() so n_times=11 requires only 9/11 (81.8%) instead of 90%. MEDIUM: NaN time steps produce gi_bin=0 which classifier counts as 'non-significant' rather than missing; threshold_90 uses full n_times not valid count. LOW: 'global_std == 0' check does not catch NaN std for fully/mostly NaN inputs.
+erosion,2026-08-08,3703,HIGH,2;3;5,"HIGH (Cat 2+3, fixed in PR for #3703): NaN/Inf in the input made grad and new_x NaN; the guard `if new_x < 1 or new_x >= width-2` is False for NaN so control fell through to int(nan) = INT64_MIN and indexed heightmap out of bounds (confirmed with NUMBA_BOUNDSCHECK=1 -> IndexError); h_diff NaN then took the erosion branch and painted NaN over the whole brush, so 1 NaN cell became 3951/4096 on numpy and 3969/4096 on cupy after 2000 droplets. Same two guards wrong in _erode_gpu_kernel. Fixed by writing both guards as inside-the-box tests plus a math.isfinite(h_diff) check; clean-input results are bit-identical (sha1 match before/after). MEDIUM (Cat 5, unfixed, needs a maintainer call): numpy and cupy run genuinely different simulations. The CPU kernel steps droplets sequentially so each sees the terrain the previous ones carved; the CUDA kernel launches one thread per droplet so they all see roughly the pristine terrain and race through cuda.atomic.add. On a 64x64 raster, iterations=5000, seed=42: max |numpy-cupy| = 261 on a 0-500 range, total volume change -201k vs -352k, and cupy is not even reproducible run-to-run at a fixed seed (max diff 634). The erode() docstring documents neither; `seed` reads as a reproducibility guarantee. Fix is either a docstring contract or a batched GPU launch, not something to decide inside an accuracy sweep. LOW (Cat 5): erode() always returns float32, so a float64 DEM is silently downcast. LOW (Cat 4): the simulation is cell-index-space only, so it ignores attrs['res'] entirely and erodes isotropically on anisotropic cells; same shape as the sky_view_factor bug #3626, but erode() makes no physical-units claim so this is a feature question. Cat 1 clean: accumulation is float64 with one cast at the end, no unguarded division. Cat 6: no reference tool for droplet erosion; osgeo-unavailable richdem-unavailable, validated by invariants (flat terrain unchanged, nodata confinement, seed determinism on CPU) instead. TEST-COVERAGE GAP: test_erosion.py had zero NaN/Inf tests on any backend before this run (7 added), and still has no numpy-vs-cupy value comparison, which is why the Cat 5 divergence above was never noticed."
 fire,2026-06-19,3394,MEDIUM,5,"Cat5: dask+numpy map_blocks declared float64 (meta default) while ngjit kernels return float32; numpy/cupy/dask+cupy all float32. Fixed 6 wrappers with dtype=np.float32 (PR #3396); bsc already dtype=int8. Cats 1-4 clean: per-pixel ops, no stencil/accumulation/projected-distance; NaN via x!=x; CUDA bounds strict <; rdnbr/ros divisions guarded. cupy+dask+cupy tests run on GPU host."
 flood,2026-06-25,3499,MEDIUM,5,"Cat5 backend dtype divergence (#3499/this PR): flood_depth and curve_number_runoff document float64 output; numpy/cupy cast to float64 but _flood_depth_dask/_cn_runoff_dask skipped the cast, so float32 input leaked float32 on dask + dask+cupy (numpy/cupy=float64, dask/dask+cupy=float32). Confirms the 2026-04-30 note. Fix: cast hand/p to float64 at the top of both dask helpers; dask+cupy wrappers reuse them so all 4 backends now return float64. Other flood fns unaffected: travel_time/flood_depth_vegetation upcast via float64 _TAN_MIN clamp, inundation via 1.0/0.0 literals, vegetation_roughness via np.interp. 8 new dtype tests across all 4 backends. Cats 1-4 clean; pure vectorized numpy/cupy/dask, no numba/cuda kernels, no neighborhood stencil, no geodesic math. CUDA available; cupy + dask+cupy verified (96 tests pass). LOW (not fixed, documented): curve_number/mannings_n DataArray inputs bypass scalar range validation (CN in (0,100], n>0)."
 focal,2026-06-10,3214,MEDIUM,1;5,"mean() dtype divergence: numpy/dask+numpy cast to float64 (astype(float)) while cupy/dask+cupy forced float32, so output dtype was backend-dependent and float64 rasters lost precision on GPU (offset 1e7: GPU error 0.58 > true spread 0.42, same class as fixed #2831). mean() was left out of the #2769 _promote_float contract that apply/focal_stats follow. Fix #3214: _promote_float in mean(), drop hardcoded cupy.float32 in _mean_cupy/_mean_dask_cupy, excludes cast to working dtype for cross-backend match parity. CUDA available; all 4 backends executed (245 focal tests pass incl new 3214 dtype tests). Cats 2-4 clean: GPU kernels two-pass std/var (#2831 fix verified), NaN checks via v!=v, map_overlap depths == kernel radius, Gi* validated against reference test. LOW (documented, not fixed): mean() excludes mask only the center pixel; excluded sentinel values (e.g. -9999) still contribute to neighboring cells' means on all backends -- docstring says 'left unchanged rather than averaged', backend-consistent."

xrspatial/erosion.py

@@ -1,5 +1,7 @@
 from __future__ import annotations
 
+import math
+
 import numpy as np
 import xarray as xr
 from numba import jit
@@ -119,15 +121,25 @@ def _erode_cpu(heightmap, random_pos, boy, box, bw,
             dir_y = dir_y * inertia - grad_y * (1 - inertia)
 
             dir_len = (dir_x * dir_x + dir_y * dir_y) ** 0.5
-            if dir_len < 1e-10:
+            # A nodata cell in the stencil above makes the gradient, and
+            # therefore dir_len, non-finite.  `dir_len < 1e-10` is False for
+            # NaN, so the finite check has to be spelled out or the droplet
+            # runs on with a NaN position.
+            if dir_len < 1e-10 or not math.isfinite(dir_len):
                 break
             dir_x /= dir_len
             dir_y /= dir_len
 
             new_x = pos_x + dir_x
             new_y = pos_y + dir_y
 
-            if new_x < 1 or new_x >= width - 2 or new_y < 1 or new_y >= height - 2:
+            # Every comparison against NaN is False, so this guard has to be
+            # written as "inside the valid box" and break on anything else.
+            # Written the other way round a NaN falls through to int(new_x),
+            # which numba evaluates to INT64_MIN, and the reads below index
+            # outside the array.
+            if not (new_x >= 1 and new_x < width - 2
+                    and new_y >= 1 and new_y < height - 2):
                 break
 
             h_old = h00 * (1 - fx) * (1 - fy) + h10 * fx * (1 - fy) + \
@@ -144,6 +156,14 @@ def _erode_cpu(heightmap, random_pos, boy, box, bw,
 
             h_diff = h_new - h_old
 
+            # Nodata in either stencil leaves h_diff non-finite.  Both branch
+            # conditions below are False for NaN, so control would reach the
+            # erosion branch and subtract NaN from every cell under the brush,
+            # spreading nodata across the grid one droplet at a time.  Kill
+            # the droplet instead: nodata acts as a barrier.
+            if not math.isfinite(h_diff):
+                break
+
             sed_capacity = max(-h_diff, min_slope) * speed * water * capacity
 
             if sediment > sed_capacity or h_diff > 0:
@@ -237,17 +257,21 @@ def _erode_gpu_kernel(
             dir_y = dir_y * inertia - grad_y * (1 - inertia)
 
             dir_len = (dir_x * dir_x + dir_y * dir_y) ** 0.5
-            if dir_len < 1e-10:
+            # See the matching comment in _erode_cpu: `dir_len < 1e-10` is
+            # False for a NaN gradient coming out of a nodata stencil.
+            if dir_len < 1e-10 or not math.isfinite(dir_len):
                 return
             dir_x /= dir_len
             dir_y /= dir_len
 
             new_x = pos_x + dir_x
             new_y = pos_y + dir_y
 
-            if new_x < 1 or new_x >= width - 2:
+            # Written as a pair of rejection tests a NaN position slips
+            # through both and the reads below index outside the array.
+            if not (new_x >= 1 and new_x < width - 2):
                 return
-            if new_y < 1 or new_y >= height - 2:
+            if not (new_y >= 1 and new_y < height - 2):
                 return
 
             h_old = (h00 * (1 - fx) * (1 - fy) + h10 * fx * (1 - fy) +
@@ -264,6 +288,11 @@ def _erode_gpu_kernel(
 
             h_diff = h_new - h_old
 
+            # Nodata anywhere in either stencil: kill the particle rather than
+            # letting the erosion branch atomically add NaN to the brush.
+            if not math.isfinite(h_diff):
+                return
+
             neg_h_diff = -h_diff
             if neg_h_diff < min_slope:
                 neg_h_diff = min_slope
@@ -460,6 +489,13 @@ def erode(agg, iterations=50000, seed=42, params=None):
         is outside the allowed range.
     MemoryError
         If the projected working set exceeds available memory.
+
+    Notes
+    -----
+    Non-finite cells act as barriers. A droplet whose interpolation stencil
+    covers a NaN or an Inf dies there without touching the heightmap, so
+    nodata cells come back unchanged and the finite terrain around them
+    erodes normally.
     """
     _validate_scalar(
         iterations, func_name='erode', name='iterations',

xrspatial/tests/test_erosion.py

@@ -214,6 +214,92 @@ def test_erode_dask_cupy_runs():
     assert np.isfinite(result_np).all()
 
 
+# ---- nodata handling (issue #3703) ----
+
+def _to_numpy(result):
+    data = result.data
+    if hasattr(data, 'compute'):
+        data = data.compute()
+    if hasattr(data, 'get'):
+        data = data.get()
+    return data
+
+
+def _check_nodata_confined(backend, chunks=(32, 32), nodata=np.nan):
+    """A single non-finite cell must stay a single non-finite cell.
+
+    Before #3703 the droplet ran on past a NaN gradient, indexed the
+    heightmap with int(nan) (out of bounds), and painted NaN across the
+    brush footprint, so one nodata cell grew to cover most of the raster.
+    """
+    data = _make_terrain(size=64)
+    data[30, 30] = nodata
+    agg = _input(data, backend, chunks=chunks)
+    result = _to_numpy(erode(agg, iterations=2000, seed=42))
+
+    np.testing.assert_array_equal(
+        result[30, 30], np.float32(nodata),
+        err_msg="the nodata cell was overwritten",
+    )
+    spread = int((~np.isfinite(result)).sum())
+    assert spread == 1, f"nodata spread to {spread} cells"
+
+    # The rest of the raster still erodes normally.
+    finite = np.isfinite(result)
+    changed = int((result[finite] != data[finite]).sum())
+    assert changed > 100, f"only {changed} finite cells changed"
+
+
+def test_erode_nodata_confined_numpy():
+    _check_nodata_confined('numpy')
+
+
+@dask_array_available
+def test_erode_nodata_confined_dask_numpy():
+    _check_nodata_confined('dask+numpy', chunks=(16, 16))
+
+
+@cuda_and_cupy_available
+def test_erode_nodata_confined_cupy():
+    _check_nodata_confined('cupy')
+
+
+@cuda_and_cupy_available
+@dask_array_available
+def test_erode_nodata_confined_dask_cupy():
+    _check_nodata_confined('dask+cupy', chunks=(16, 16))
+
+
+def test_erode_inf_confined_numpy():
+    """An Inf cell behaves the same way a NaN cell does: it stays put."""
+    _check_nodata_confined('numpy', nodata=np.inf)
+
+
+@cuda_and_cupy_available
+def test_erode_inf_confined_cupy():
+    _check_nodata_confined('cupy', nodata=np.inf)
+
+
+def test_erode_all_nodata_raster():
+    """An all-nodata raster comes back all-nodata instead of crashing."""
+    data = np.full((32, 32), np.nan, dtype=np.float32)
+    result = erode(_input(data, 'numpy'), iterations=1000, seed=42)
+    assert np.isnan(result.data).all()
+
+
+def test_erode_nodata_border_leaves_interior_intact():
+    """A nodata border must not eat into the terrain behind it."""
+    data = _make_terrain(size=48)
+    data[:2, :] = np.nan
+    data[-2:, :] = np.nan
+    data[:, :2] = np.nan
+    data[:, -2:] = np.nan
+    expected_nan = int(np.isnan(data).sum())
+
+    result = erode(_input(data, 'numpy'), iterations=3000, seed=42)
+    assert int(np.isnan(result.data).sum()) == expected_nan
+
+
 # ---- parameter validation (issue #1275) ----
 
 def test_erode_iterations_zero_rejected():