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

[brendancol]

Closes #3703

Both erosion kernels wrote the moved-particle bounds check as a pair of rejection tests:

if new_x < 1 or new_x >= width - 2 or new_y < 1 or new_y >= height - 2:
    break

Every comparison against NaN is False, so a droplet whose bilinear stencil contained a nodata cell fell straight through to int(new_x). numba evaluates int(nan) as -9223372036854775808, so the four heightmap reads underneath indexed far outside the array. With the default boundscheck=False that is a silent read of arbitrary process memory; under NUMBA_BOUNDSCHECK=1 it raises IndexError.

h_diff was NaN at that point too. Both branch conditions are False for NaN, so control reached the erosion branch and subtracted NaN from every cell under the brush. One nodata cell turned into 3951 of 4096 on numpy and 3969 of 4096 on cupy after 2000 droplets.

Changes

• Rewrite both bounds guards as inside-the-box tests, in _erode_cpu and _erode_gpu_kernel.
• Add a math.isfinite(h_diff) check before the deposit/erode branch in both kernels. Nodata acts as a barrier now: the droplet dies without touching the heightmap, the NaN cells stay NaN, and the finite part of the raster erodes as before.
• Seven tests in xrspatial/tests/test_erosion.py, which had no NaN or Inf coverage on any backend.
• Record the erosion row in .claude/sweep-accuracy-state.csv.

Backends

numpy, cupy, dask+numpy and dask+cupy all go through the same two kernels, and the fix is applied to both. All four are covered by the new confinement test and were run on a CUDA host.

Verification

Rasters with no nodata are unaffected. Same input, iterations=5000, seed=42, sha1 of the output buffer before and after the fix:

before: b2d576aa19335351799afef7c8c81e15ef7bb886
after:  b2d576aa19335351799afef7c8c81e15ef7bb886

The issue reproduction, on both backends after the fix:

numpy: input nan=1 -> output nan=1, nan at (30,30)=True
numpy: input inf=1 -> output non-finite=1
numpy: cells changed outside nodata=4088
cupy:  input nan=1 -> output nan=1, nan at (30,30)=True
cupy:  input inf=1 -> output non-finite=1
cupy:  cells changed outside nodata=4091

And the out-of-bounds access is gone:

$ NUMBA_BOUNDSCHECK=1 python repro.py
no boundscheck error; nan out: 1

Test plan

• pytest xrspatial/tests/test_erosion.py - 35 passed, cupy and dask+cupy included
• pytest xrspatial/tests/test_terrain.py - 57 passed (generate_terrain(erode=True) calls this)
• flake8 on both changed files adds no new findings

Not addressed here

The accuracy sweep also turned up a numpy/cupy divergence that is out of scope for this fix. The CPU kernel steps droplets one at a time, 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 at iterations=5000, seed=42 the two backends differ by up to 261 on a 0-500 elevation range, and cupy is not reproducible run-to-run at a fixed seed. Choosing between documenting that in the erode() docstring and batching the GPU launch is a maintainer call, so I left it alone and wrote it up in the sweep state notes.

[brendancol]

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

Reviewed the full erosion.py and test_erosion.py, not just the diff. The diagnosis holds up: int(nan) is INT64_MIN in numba, both bounds guards were written as rejection tests, and NaN makes every comparison False. Both kernels get the same two fixes, which is the part that matters most here since numpy, cupy, dask+numpy and dask+cupy all funnel into these two functions.

Blockers (must fix before merge)

None.

Suggestions (should fix, not blocking)

• xrspatial/erosion.py:122 -- the direction guard is still NaN-blind:

  dir_len = (dir_x * dir_x + dir_y * dir_y) ** 0.5
  if dir_len < 1e-10:
      break

  A NaN gradient produces a NaN dir_len, nan < 1e-10 is False, and the droplet carries on to divide by NaN. The PR is still correct because the new new_x / new_y guard eight lines down catches the resulting NaN position, but that makes the memory safety of line 140 depend on a guard somewhere else. This is exactly the shape of the bug being fixed. Worth making the guard stand on its own: if dir_len < 1e-10 or not math.isfinite(dir_len). Same at erosion.py:257 in the CUDA kernel. It also saves a few wasted arithmetic steps per dead droplet.

• xrspatial/erosion.py:459 -- erode() now has defined nodata behaviour (barrier, preserved in the output), and the docstring does not say so. Callers running this on a real DEM will want to know that masked cells survive and that the terrain behind them still erodes. A couple of lines under Parameters or a short Notes section covers it.

Nits (optional improvements)

• xrspatial/tests/test_erosion.py:255 -- test_erode_inf_confined runs on numpy only. The CUDA kernel got the same guard change and has no Inf coverage. Reusing the _check_nodata_confined shape for Inf would close that.

• Informational, not a request: the brush loop at erosion.py:181 still subtracts from cells that are NaN. Harmless for the raster (NaN minus a finite number is NaN, so nodata is preserved), but sediment += amount credits the droplet for material it removed from nodata, so a droplet eroding right at the edge of a mask picks up slightly more sediment than it should. Pre-existing and untouched by this PR. Not worth chasing here.

What looks good

• The clean-raster sha1 check in the PR body is the right evidence. Rewriting a guard's polarity is easy to get subtly wrong, and a bit-identical result on a NaN-free input rules that out.
• _check_nodata_confined asserts both halves of the fix: the NaN count stays at 1 and more than 100 finite cells still changed. Asserting only the first would have passed on a kernel that stopped eroding entirely.
• All four backends are covered by the confinement test and were run on a CUDA host, so the CUDA guard change is verified rather than assumed.
• test_erode_nodata_border_leaves_interior_intact is a good addition. A masked border is what real DEMs actually look like, and it exercises many more droplet paths against the mask than a single interior cell does.

Checklist

• Algorithm matches reference/paper -- no algorithm change; clean-input output is bit-identical
• All implemented backends produce consistent results -- for nodata confinement, yes; the wider numpy/cupy divergence the PR body calls out is pre-existing and out of scope
• NaN handling is correct -- this is the fix; verified on numpy and cupy
• Edge cases are covered by tests -- single NaN, all-NaN, Inf, masked border
• Dask chunk boundaries handled correctly -- erode() materializes before running, so there are no chunk boundaries
• No premature materialization or unnecessary copies -- two predicates added inside an existing loop
• Benchmark exists or is not needed -- not needed, no new function
• README feature matrix updated (if applicable) -- not applicable
• Docstrings present and accurate -- see the second suggestion

[brendancol]

PR Review: follow-up pass on 96699ac

Re-reviewed after the review commit. All three actionable findings from the first pass are addressed.

Disposition of the previous findings

• Fixed -- erosion.py:122 / erosion.py:262: both direction guards are now if dir_len < 1e-10 or not math.isfinite(dir_len). The NaN death happens at the site that produces it instead of eight lines downstream, so neither guard depends on the other for memory safety.
• Fixed -- erode() now has a Notes section stating that non-finite cells act as barriers and come back unchanged. That is the contract the tests assert, so it belongs in the docstring.
• Fixed -- test_erode_inf_confined_cupy closes the Inf gap on the CUDA kernel. Folding Inf into _check_nodata_confined via a nodata parameter was the right call; it means the Inf case now also asserts the finite terrain still erodes, which the old standalone test did not.
• Dismissed -- the brush loop crediting sediment for material "removed" from nodata cells. Pre-existing behaviour, untouched by this PR, and fixing it means deciding what a droplet should do when its brush straddles a mask edge. That is a modelling question, not a bug fix.

Blockers

None.

Suggestions

None.

Nits

None.

Notes on the follow-up commit

np.testing.assert_array_equal(result[30, 30], np.float32(nodata)) is a neat way to cover both cases in one assertion, since assert_array_equal treats NaN as equal to NaN. Worth being aware that it would also pass if the kernel replaced a NaN with a different NaN payload, but that is not a distinction anything here cares about.

Re-ran after the commit:

• pytest xrspatial/tests/test_erosion.py xrspatial/tests/test_terrain.py -- 93 passed, cupy and dask+cupy included
• clean-raster sha1 still b2d576aa19335351799afef7c8c81e15ef7bb886, so the extra guard did not perturb NaN-free results either
• flake8 on both files reports only findings that already exist on main