surface_distance(): the Dijkstra heap is sized H*W, so dense source rasters write out of bounds and abort the process

[brendancol]

Summary

_dijkstra and _dijkstra_geodesic in xrspatial/surface_distance.py size the
binary-heap backing arrays as

max_heap = height * width          # lines 214 and 277
h_keys = np.empty(max_heap, dtype=np.float64)
h_rows = np.empty(max_heap, dtype=np.int64)
h_cols = np.empty(max_heap, dtype=np.int64)

Both kernels are lazy-deletion Dijkstra: a pixel is pushed onto the heap every
time its tentative distance improves, so it can be enqueued once per settled
neighbour rather than once overall. The push count is bounded by
height * width * (n_neighbors + 1), not height * width.

_heap_push (in cost_distance.py) writes keys[pos] at pos = size with no
bounds check. Once the heap outgrows the allocation, the kernel writes past the
end of three numba arrays. numba compiles without bounds checking by default, so
nothing catches it.

cost_distance.py carried the same defect and was already corrected; see the
comment at xrspatial/cost_distance.py:155-161 and the sizing at line 161.
surface_distance.py was never updated.

Reproduction

A checkerboard source raster over random elevation, 8-connectivity, on the
plain numpy backend:

import numpy as np
import xarray as xr
from xrspatial.surface_distance import surface_distance

rng = np.random.default_rng(0)
H = W = 80
elev = rng.random((H, W)) * 5000.0
rr, cc = np.meshgrid(np.arange(H), np.arange(W), indexing='ij')
src = ((rr + cc) % 2 == 0).astype(np.float64)

coords = {'y': np.arange(H, dtype=float), 'x': np.arange(W, dtype=float)}
raster = xr.DataArray(src, dims=('y', 'x'), coords=coords)
elevation = xr.DataArray(elev, dims=('y', 'x'), coords=coords)

surface_distance(raster, elevation, connectivity=8)

Default build (no bounds checking) — the process dies:

$ python repro.py
double free or corruption (!prev)
Aborted (core dumped)

Same script under NUMBA_BOUNDSCHECK=1, which names the defect:

$ NUMBA_BOUNDSCHECK=1 python repro.py
  File "xrspatial/surface_distance.py", line 474, in _surface_distance_numpy
    _dijkstra(elev_data, H, W, max_distance,
              dy, dx, dd, dist, alloc, src_row, src_col)
IndexError: index is out of bounds

Affected backends

Reproduced under NUMBA_BOUNDSCHECK=1 with the same 80x80 input:

Path	Result
numpy, planar (_dijkstra)	IndexError: index is out of bounds
numpy, geodesic (_dijkstra_geodesic)	IndexError: index is out of bounds
dask+numpy, iterative tile path via _run_tile	IndexError: index is out of bounds

The dask+numpy bounded map_overlap branch runs _surface_distance_numpy per
chunk, so it is exposed on any chunk whose seeded pixel density is high enough.

The cupy and dask+cupy backends use iterative parallel relaxation with no heap
and are not affected. compute-sanitizer --tool memcheck over nine adversarial
shapes (2x2, 2x97, 97x2, 3x101, 101x3, 61x59, 127x127, 33x512, 512x33) across all
three public functions reports ERROR SUMMARY: 0 errors, and the
_sd_relax_kernel bounds guard at line 496 is correct.

How often it bites

Peak concurrent heap size against the allocated cap, measured with an
instrumented copy of _dijkstra (oversized heap so it can run to completion),
8-connectivity, random elevation:

Raster	Seeding	Peak heap	Allocated cap	Overflow
60x60	single source	452	3600	no (0.13x)
100x100	20% random sources	9810	10000	no (0.98x)
40x40	checkerboard	1614	1600	yes (1.01x)
80x80	checkerboard	6646	6400	yes (1.04x)
200x200	checkerboard	41711	40000	yes (1.04x)

Sparse single-source rasters stay well clear, which is why the existing tests
pass. Dense source masks are the trigger, and those are ordinary usage:
a rasterized road or stream network, a landcover class mask, or any
target_values= selection that matches a large share of the raster. A 20%
random source density already reaches 0.98x of the cap.

Issue #3709 notes that the module's own benchmark raster makes 99.95% of pixels
sources, so the benchmark suite is already running in the regime that overflows.

Proposed fix

Size the heap to the actual bound in both kernels, matching what
cost_distance.py already does:

max_heap = height * width * (n_neighbors + 1)

The _BYTES_PER_PIXEL = 80 constant used by _check_memory models the heap at
24 bytes/pixel, so it needs to follow the new allocation or the guard added in
#1305 will under-count the eager numpy footprint by roughly 3.4x.

Related, not filed

xrspatial/pathfinding.py:304 uses the same max_heap = height * width sizing
with the same _heap_push. That kernel is single-source A*, and I could not get
it to overflow, so I am recording it here as a site to re-check rather than
claiming a defect.

---

Found by /deep-sweep security audit of surface_distance on 2026-08-16.