erode(): NaN input causes an out-of-bounds read and floods the output with NaN

[brendancol]

Description

Neither erosion kernel checks for NaN. A single nodata cell in the input causes two problems:

1. An out-of-bounds array index inside _erode_cpu, and the same code in _erode_gpu_kernel. With numba's default boundscheck=False this is a silent read of arbitrary process memory.
2. The output fills with NaN. On a 64x64 raster, one NaN cell becomes 3951 NaN cells after 2000 droplets, 96% of the grid.

Real DEMs carry nodata as NaN, so this is reachable from ordinary use, including generate_terrain(..., erode=True) on a masked template.

Mechanism

In xrspatial/erosion.py, the bilinear stencil at lines 110-116 reads four neighbours and forms a gradient. If any of them is NaN, grad_x / grad_y are NaN, so dir_x and dir_y are NaN and dir_len is NaN.

The guard at line 122 is

if dir_len < 1e-10:
    break

nan < 1e-10 is False, so the droplet keeps going. new_x and new_y become NaN, and the next guard at line 130 has the same problem:

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

Both comparisons are False for NaN, so execution falls through to line 136:

new_node_x = int(new_x)
new_node_y = int(new_y)
h_new = (heightmap[new_node_y, new_node_x] * ...

In numba, int(nan) is -9223372036854775808, so those four reads index far outside the array.

The NaN flood is a separate consequence. h_diff is NaN, so sediment > sed_capacity and h_diff > 0 are both False and control reaches the erosion branch at line 162, which subtracts amount * bw[k] (NaN) from every cell under the brush. Each droplet that touches the NaN region paints a disk of NaN roughly (2r+1)**2 cells wide, and later droplets spread it further.

The CUDA kernel at lines 240 and 248-251 has the identical guards and the identical fallthrough.

Reproduction

import numpy as np, xarray as xr
from xrspatial.erosion import erode

rng = np.random.default_rng(0)
data = (rng.random((64, 64)).astype(np.float32) * 500)
data[30, 30] = np.nan   # single nodata cell
agg = xr.DataArray(data, dims=['y', 'x'],
                   coords={'y': np.arange(64.0), 'x': np.arange(64.0)})
print("input nan count:", np.isnan(data).sum())
res = erode(agg, iterations=2000, seed=42)
print("output nan count:", np.isnan(res.data).sum())

Output on the numpy backend:

input nan count: 1
output nan count: 3951

The same script on the cupy backend gives 3969 NaN cells out of 4096.

Running the numpy version under NUMBA_BOUNDSCHECK=1 shows the out-of-bounds access directly:

$ NUMBA_BOUNDSCHECK=1 python repro.py
IndexError: index is out of bounds

And the int(nan) behaviour on its own:

from numba import jit
import numpy as np

@jit(nopython=True)
def g(a, x):
    return a[int(x)]

print(g(np.arange(10.0), np.nan))   # -> 4.74e-322, an out-of-bounds read

Expected behaviour

A nodata cell should act as a barrier. A droplet whose stencil contains NaN should die without touching the heightmap, so the NaN cells stay NaN and the finite part of the raster erodes normally. And nothing should index outside the array.

Environment

• xarray-spatial main (a873b49)
• numba 0.x with default boundscheck=False
• Linux, CUDA available, reproduced on numpy, dask+numpy, cupy and dask+cupy

Notes

xrspatial/tests/test_erosion.py has no NaN test on any backend, which is why this survived five commits.

Found by the accuracy sweep (Cat 2: NaN propagation; Cat 3: bounds guard).