From c5f21b19792ac412d4bd3b12673b921d44c7170e Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Sat, 8 Aug 2026 00:38:39 -0400 Subject: [PATCH 1/2] Stop NaN input from reading out of bounds and flooding the erosion output (#3703) Both erosion kernels wrote the moved-particle bounds check as a pair of rejection tests. Every comparison against NaN is False, so a droplet whose stencil contained a nodata cell fell through to int(new_x), which numba evaluates to INT64_MIN, and the four heightmap reads below indexed outside the array. h_diff was then NaN, which sent control to the erosion branch and subtracted NaN from every cell under the brush, so a single nodata cell grew to cover 96% of a 64x64 raster after 2000 droplets. Rewrite both guards as inside-the-box tests and add a math.isfinite(h_diff) check before the deposit/erode branch, in the numba kernel and the CUDA kernel. Nodata now acts as a barrier: the droplet dies, the NaN cells stay put, and the finite part of the raster erodes as before. Results on a raster with no nodata are bit-identical. Also records the erosion row in the accuracy sweep state CSV. --- .claude/sweep-accuracy-state.csv | 1 + xrspatial/erosion.py | 31 ++++++++++-- xrspatial/tests/test_erosion.py | 82 ++++++++++++++++++++++++++++++++ 3 files changed, 111 insertions(+), 3 deletions(-) diff --git a/.claude/sweep-accuracy-state.csv b/.claude/sweep-accuracy-state.csv index 47b4bc4f9..aef62f716 100644 --- a/.claude/sweep-accuracy-state.csv +++ b/.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." diff --git a/xrspatial/erosion.py b/xrspatial/erosion.py index 23d664f13..b696d9ebe 100644 --- a/xrspatial/erosion.py +++ b/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 @@ -127,7 +129,14 @@ def _erode_cpu(heightmap, random_pos, boy, box, bw, 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: + # A nodata cell in the stencil above makes the gradient, and + # therefore new_x / new_y, NaN. Every comparison against NaN is + # False, so the 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 +153,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: @@ -245,9 +262,12 @@ def _erode_gpu_kernel( new_x = pos_x + dir_x new_y = pos_y + dir_y - if new_x < 1 or new_x >= width - 2: + # See the matching comment in _erode_cpu: 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 +284,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 diff --git a/xrspatial/tests/test_erosion.py b/xrspatial/tests/test_erosion.py index 9b8ed7a84..f9042bf88 100644 --- a/xrspatial/tests/test_erosion.py +++ b/xrspatial/tests/test_erosion.py @@ -214,6 +214,88 @@ 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)): + """A single NaN cell must stay a single NaN 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] = np.nan + agg = _input(data, backend, chunks=chunks) + result = _to_numpy(erode(agg, iterations=2000, seed=42)) + + assert np.isnan(result[30, 30]), "the nodata cell was overwritten" + assert int(np.isnan(result).sum()) == 1, ( + f"nodata spread to {int(np.isnan(result).sum())} 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(): + """An Inf cell behaves the same way a NaN cell does: it stays put.""" + data = _make_terrain(size=64) + data[30, 30] = np.inf + result = erode(_input(data, 'numpy'), iterations=2000, seed=42) + assert int((~np.isfinite(result.data)).sum()) == 1 + + +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(): From 96699ac4fd62a770750715f7632b3b6b896add48 Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Sat, 8 Aug 2026 00:41:37 -0400 Subject: [PATCH 2/2] Address review: make the dir_len guard NaN-safe on its own (#3703) The direction guard still read `dir_len < 1e-10`, which is False for the NaN that a nodata stencil produces. The fix was still correct because the new_x / new_y guard downstream caught the resulting NaN position, but that left the memory safety of the interpolation reads depending on a guard eight lines away, which is the shape of the original bug. Both kernels now check math.isfinite(dir_len) at the same site. Also documents the nodata contract in the erode() docstring and extends the Inf test to cupy, which had no Inf coverage. --- xrspatial/erosion.py | 33 ++++++++++++++++++++++----------- xrspatial/tests/test_erosion.py | 26 +++++++++++++++----------- 2 files changed, 37 insertions(+), 22 deletions(-) diff --git a/xrspatial/erosion.py b/xrspatial/erosion.py index b696d9ebe..4dfe8f923 100644 --- a/xrspatial/erosion.py +++ b/xrspatial/erosion.py @@ -121,7 +121,11 @@ 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 @@ -129,12 +133,11 @@ def _erode_cpu(heightmap, random_pos, boy, box, bw, new_x = pos_x + dir_x new_y = pos_y + dir_y - # A nodata cell in the stencil above makes the gradient, and - # therefore new_x / new_y, NaN. Every comparison against NaN is - # False, so the 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. + # 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 @@ -254,7 +257,9 @@ 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 @@ -262,9 +267,8 @@ def _erode_gpu_kernel( new_x = pos_x + dir_x new_y = pos_y + dir_y - # See the matching comment in _erode_cpu: written as a pair of - # rejection tests a NaN position slips through both and the reads - # below index outside the array. + # 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 not (new_y >= 1 and new_y < height - 2): @@ -485,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', diff --git a/xrspatial/tests/test_erosion.py b/xrspatial/tests/test_erosion.py index f9042bf88..4c6524cfc 100644 --- a/xrspatial/tests/test_erosion.py +++ b/xrspatial/tests/test_erosion.py @@ -225,22 +225,24 @@ def _to_numpy(result): return data -def _check_nodata_confined(backend, chunks=(32, 32)): - """A single NaN cell must stay a single NaN cell. +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] = np.nan + data[30, 30] = nodata agg = _input(data, backend, chunks=chunks) result = _to_numpy(erode(agg, iterations=2000, seed=42)) - assert np.isnan(result[30, 30]), "the nodata cell was overwritten" - assert int(np.isnan(result).sum()) == 1, ( - f"nodata spread to {int(np.isnan(result).sum())} cells" + 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) @@ -268,12 +270,14 @@ def test_erode_nodata_confined_dask_cupy(): _check_nodata_confined('dask+cupy', chunks=(16, 16)) -def test_erode_inf_confined(): +def test_erode_inf_confined_numpy(): """An Inf cell behaves the same way a NaN cell does: it stays put.""" - data = _make_terrain(size=64) - data[30, 30] = np.inf - result = erode(_input(data, 'numpy'), iterations=2000, seed=42) - assert int((~np.isfinite(result.data)).sum()) == 1 + _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():