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..4dfe8f923 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 @@ -119,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 @@ -127,7 +133,13 @@ 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: + # 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,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 @@ -245,9 +267,11 @@ 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: + # 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', diff --git a/xrspatial/tests/test_erosion.py b/xrspatial/tests/test_erosion.py index 9b8ed7a84..4c6524cfc 100644 --- a/xrspatial/tests/test_erosion.py +++ b/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():