Let the GPU surface-distance relaxation run to convergence

[brendancol]

Closes #3721

What changed

_surface_distance_cupy() solves the shortest-path problem by repeated
parallel relaxation and capped itself at H + W passes. One pass moves
distance information roughly one pixel hop, so what the loop actually needs is
the hop count of the longest shortest path. The two numbers agree on open
terrain and come apart as soon as NaN barriers make paths wind; the
Bellman-Ford bound for a grid graph is the pixel count, not the perimeter.

Reaching the cap exited the loop with no warning and nothing in the output to
say so. Pixels the CPU backend reaches came back NaN, which a caller reads as
"no path exists". On a 16x16 serpentine corridor the GPU found 37 of the 136
reachable pixels; at 24x24 it found 53 of 300.

The ceiling is now H * W, and reaching it warns instead of quietly handing
back a truncated answer.

Throughput

The changed flag still ends the loop as soon as the solution settles, so
ordinary terrain runs exactly as many passes as it did before. Timed on this
host against main:

raster	main	this branch
256x256 random elevation	0.136 s	0.135 s
512x512 random elevation	0.326 s	0.313 s

Both inside run-to-run noise.

Backends

cupy, and the bounded dask+cupy path, which calls this function per chunk
against that chunk's H and W (so the old cap bit even harder there).
numpy and dask+numpy run an exact Dijkstra and were already correct.

Test plan

• pytest xrspatial/tests/test_surface_distance.py — 39 passed (37 before)
• flake8 xrspatial/surface_distance.py clean
• New test_cupy_long_path_matches_numpy, parametrized over
  surface_distance and surface_allocation, builds a serpentine
  corridor and asserts cupy matches numpy. It fails against main's
  surface_distance.py and passes here.
• Checked the repro from the issue at N = 8, 12, 16, 24: the GPU now
  reports the same reachable set as numpy at every size, with zero
  difference on the distances.

Found by /sweep-accuracy against xrspatial/surface_distance.py.

[brendancol]

PR Review: Let the GPU surface-distance relaxation run to convergence

Read the whole of _surface_distance_cupy and _sd_relax_kernel rather than
just the diff, and replayed the issue repro at N = 8, 12, 16, 24 on this host.
The diagnosis holds: H + W is the perimeter, not the hop count of the
longest shortest path, and H * W is the right ceiling for a grid graph with
H * W nodes. Since the in-place relaxation converges at least as fast as
Jacobi Bellman-Ford, the new bound is provably sufficient, which makes the
warning unreachable in theory. That is the correct shape for a safety net.

Blockers (must fix before merge)

None.

Suggestions (should fix, not blocking)

• xrspatial/tests/test_surface_distance.py:565 — the long-path test skips
  DIRECTION. It is parametrized over surface_distance and
  surface_allocation, but _surface_distance_cupy runs the same relaxation
  loop for all three modes and DIRECTION reads srow/scol, which the
  truncated loop leaves unset just as it leaves dist at infinity. Adding
  surface_direction to the parametrize list costs one word and covers the
  third mode.

• xrspatial/surface_distance.py:608 — stacklevel=2 points inside the
  library. From a user's surface_distance() call the chain is the
  supports_dataset wrapper, _compute, then _surface_distance_cupy, so
  level 2 lands on _compute. The module's other warning
  (_surface_distance_dask) uses stacklevel=4 for the same call depth.
  Worth matching so the warning points at the caller.

Nits (optional improvements)

• Runtime on a genuinely pathological raster is now unbounded in practice.
  A raster whose shortest paths really do need H * W passes will run H * W
  kernel launches, and on anything large that is a hang rather than an answer.
  The old cap made that case fast and wrong, which is worse, so this is the
  right trade. But nothing tells the caller they are in it. A cheap
  improvement would be a one-time note once the pass count passes, say,
  4 * (H + W), so a long run is at least explicable. I would not hold the PR
  for it.

• int(changed[0]) syncs the device to the host every pass. Pre-existing,
  and the pass count for converging input is unchanged, so this PR does not
  make it worse. Mentioning it only because raising the ceiling makes the
  per-pass cost matter more if the pathological case ever shows up.

• The warning may not surface from a dask worker. The bounded dask+cupy
  path calls this function inside a map_overlap task, where warnings.warn
  can be swallowed depending on the scheduler. Same caveat applies to the
  module's existing warnings, so no action here.

Coverage gap worth recording

There is no dask+cupy long-path test, and I do not think a meaningful one is
writable today: the unbounded dask+cupy path converts to numpy and runs the
CPU Dijkstra, and the bounded path caps path length by max_distance well
below the point where the pass count matters. The per-chunk cap is still real
(a chunk's H + W is smaller than the raster's), so the fix matters there
even though a test cannot reach it directly.

What looks good

• The failure is characterized by counts, not adjectives: 37 of 136 reachable
  pixels at 16x16, 53 of 300 at 24x24, and the distances that both backends
  did report agreed exactly. That separates truncated propagation from
  arithmetic drift, which is the distinction that matters here.
• The throughput claim is measured against main on the same host rather than
  argued from the early-exit check.
• The comment above max_iterations explains why the old bound looked
  plausible, which is more useful to the next reader than only stating the new
  one.
• _serpentine_elevation is a compact, deterministic way to build a corridor
  much longer than the raster's perimeter, and the test asserts the corridor
  really is reachable on the CPU before comparing, so it cannot pass
  vacuously.

Checklist

• Algorithm matches reference — H * W is the Bellman-Ford bound for a
  grid graph of H * W nodes
• All implemented backends produce consistent results — cupy now matches
  numpy on the maze at every size tested
• NaN handling is correct — the truncation was producing spurious NaN;
  that is what the fix removes
• Edge cases covered by tests — long winding path, the case the old cap
  could not reach
• Dask chunk boundaries — not applicable, no change to the tiling
• No premature materialization — no new host transfers
• Benchmark exists — benchmarks/benchmarks/surface_distance.py; timings
  against main are in the PR body
• README feature matrix — not applicable, no new function
• Docstrings — internal function, the reasoning is in the code comment

[brendancol]

PR Review: follow-up pass after 1240e3a

Re-read the changed hunks and re-ran the suite. Disposition of the first pass:

Fixed

• The long-path test now parametrizes over surface_direction as well as
  surface_distance and surface_allocation. All three modes share the
  relaxation loop and DIRECTION reads srow/scol, so it was the mode most
  worth adding. 40 tests pass (37 on main).
• The convergence warning moved from stacklevel=2 to 4, matching
  _surface_distance_dask's choice for the same call depth, so it points at
  the caller rather than at _compute.

Dismissed, with reasons

• Warn partway through a long run. Any threshold I picked (4 * (H + W), or
  anything else) would be a number with nothing behind it, and it would fire
  on input that is slow but correct. A warning nobody can act on is worse than
  no warning. The real answer for a raster that deep is to route it to the CPU
  Dijkstra, which is a design change and does not belong in a bug fix.
• int(changed[0]) syncs the device to the host every pass. Pre-existing,
  and the pass count for converging input is identical to main, so this PR
  does not touch it. Changing it means restructuring the termination check.
• The warning may not surface from a dask worker. Applies equally to the
  module's existing warnings; nothing specific to this change to fix.

Recorded, not fixed

• No dask+cupy long-path test. The unbounded dask+cupy path converts to numpy
  and runs the CPU Dijkstra, and the bounded path caps path length by
  max_distance well below where the pass count bites, so I could not write
  one that actually reaches the code. The per-chunk cap is still real and the
  fix still applies there.

No new findings. Nothing outstanding from my side.