Size the surface_distance Dijkstra heap to its real push bound (#3722)

[brendancol]

Fixes #3722.

What was wrong

_dijkstra and _dijkstra_geodesic sized the heap backing arrays as
height * width. Both are lazy-deletion Dijkstra: a pixel is enqueued again
every time its tentative distance improves, so it can be pushed once per settled
neighbour. The bound is height * width * (n_neighbors + 1).

_heap_push writes keys[size] with no bounds check, and numba compiles without
bounds checking by default, so once the heap outgrew the allocation the kernel
wrote past the end of three arrays.

cost_distance.py already carries the corrected sizing and a comment explaining
it (lines 155-161). surface_distance.py was never brought along.

What triggers it

Dense source masks. An 80x80 checkerboard source raster over random elevation
aborts the interpreter:

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

Under NUMBA_BOUNDSCHECK=1 the same input raises IndexError: index is out of bounds from _dijkstra. Reproduced on numpy planar, numpy geodesic, and the
dask+numpy iterative tile path.

Peak concurrent heap size against the old cap, 8-connectivity, random elevation:

Raster	Seeding	Peak heap	Old 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
never caught it. Rasterized road or stream networks, landcover class masks, and
broad target_values= selections all land in the overflowing regime.

The change

Both kernels now allocate height * width * (n_neighbors + 1) slots, matching
cost_distance.

_BYTES_PER_PIXEL goes 80 -> 272 so _check_memory keeps modelling what the
eager numpy pass actually reserves. The heap is 24 bytes per slot and there are
up to 9 slots per pixel at 8-connectivity, so the heap alone is 216 bytes/pixel.
Without this the guard added in #1305 would under-count the footprint by about
3.4x and the module would OOM instead of raising the guard's MemoryError. The
guard now rejects eager rasters above roughly 117 million pixels on a 64 GB host
and still points at max_distance= and dask in its message.

Tests

Four new tests in xrspatial/tests/test_surface_distance.py:

• test_dense_sources_do_not_overflow_heap[4] / [8] — checkerboard sources on
  numpy, both connectivities
• test_dense_sources_do_not_overflow_heap_geodesic — same for _dijkstra_geodesic
• test_dense_sources_do_not_overflow_heap_dask — same through the iterative
  dask tile path
• test_heap_bound_covers_worst_case_push_count — instruments an oversized copy
  of the kernel and asserts peak usage exceeds height * width but stays inside
  height * width * (n_neighbors + 1)

On the unpatched module the first three kill the pytest process outright with a
glibc free(): invalid pointer abort, which is the point.

Verification

• pytest xrspatial/tests/test_surface_distance.py — 42 passed (was 37)
• pytest xrspatial/tests/test_{surface_distance,cost_distance,pathfinding}.py —
  285 passed
• Original crash repro now completes clean, with and without
  NUMBA_BOUNDSCHECK=1
• cupy and dask+cupy exercised on this host (CUDA available):
  compute-sanitizer --tool memcheck over nine shapes (2x2, 2x97, 97x2, 3x101,
  101x3, 61x59, 127x127, 33x512, 512x33) across all three public functions
  reports ERROR SUMMARY: 0 errors, unchanged before and after
• flake8 clean on both files apart from a pre-existing F841 at
  test_surface_distance.py:234, left alone

Not in this PR

xrspatial/pathfinding.py:304 uses the same max_heap = height * width sizing
with the same _heap_push. It is single-source A* and I could not get it to
overflow, so it is recorded in #3722 as a site to re-check rather than changed
here.

---

From the /deep-sweep security audit of surface_distance, 2026-08-16.

[brendancol]

Heads-up for whoever reviews this: PR #3728 (issue #3723) fixes the same defect, found independently by the performance sweep in the same /deep-sweep run over surface_distance. Both PRs identify the same root cause (lazy-deletion Dijkstra pushing more entries than height * width, _heap_push writing past the end with no bounds check) and both reproduce it as an interpreter abort.

The fixes are mutually exclusive and touch the same two kernels:

• This PR uses the static height * width * (n_neighbors + 1) bound, matching what cost_distance.py already carries.
• Grow the surface_distance Dijkstra heap on demand instead of capping it at height*width (#3723) #3728 grows the heap on demand, arguing the static bound raises the per-pixel working set from 80 to roughly 272 bytes, which would make _check_memory refuse rasters about 3.4x smaller than it does today.

Only one should land. Merging both will conflict.

[brendancol]

Follow-up to the cross-reference above, now that both sweeps have reported in full. The tradeoff between #3727 and #3728 is narrower than my first comment implied, and both fixes are internally consistent:

• Size the surface_distance Dijkstra heap to its real push bound (#3722) #3727 takes the static height * width * (n_neighbors + 1) bound and also raises _BYTES_PER_PIXEL from 80 to 272, so the Guard surface_distance() against unbounded eager allocations (#1303) #1305 memory guard keeps modelling the allocation it guards. Cost: _check_memory refuses rasters roughly 3.4x smaller than it does today. Matches the precedent already set in cost_distance.py.
• Grow the surface_distance Dijkstra heap on demand instead of capping it at height*width (#3723) #3728 grows the heap on demand, keeping the documented 80 bytes-per-pixel footprint and today's raster capacity, at the cost of diverging from that precedent and adding a _heap_grow path.

So it is a straight choice: consistency with cost_distance and a provably sufficient static bound (#3727), against preserved capacity and a smaller memory footprint (#3728). Both reproduce the abort and both verify the fix under NUMBA_BOUNDSCHECK=1.

One item worth carrying forward whichever way this lands: #3722 records that pathfinding.py:304 uses the same height * width heap sizing. The security sweep could not make it overflow because that path is single-source A*, but it should be re-checked rather than assumed safe.