Add GPU (CuPy) backend for cost_distance() #905
Description
Problem
cost_distance() in xrspatial/cost_distance.py is CPU-only with a 🔄 (CPU fallback) marker for GPU inputs. The implementation uses a numba-JIT multi-source Dijkstra with a binary min-heap — inherently sequential due to the priority queue.
The dask backend uses an iterative tile-based Dijkstra with boundary seeding, which works but is slow due to multiple convergence iterations.
Proposed Approach
The standard sequential Dijkstra does not parallelise well. GPU implementations require a different algorithm:
- Delta-stepping (Meyer & Sanders): well-suited to GPU and proven effective for grid graphs with bounded edge weights. Processes vertices in "buckets" of similar distance, with each bucket processed in parallel.
- Parallel Bellman-Ford variant: simpler but requires more iterations to converge.
- dask+cupy: use the same tile-based iterative approach as dask+numpy but with cupy kernels per tile.
Complexity
This is the most architecturally complex change in the roadmap — it requires implementing a non-trivial parallel shortest-path algorithm on CUDA with careful synchronisation. Consider a separate design spike before implementation.
Impact
Cost-distance analysis is fundamental to environmental planning, transportation analysis, and habitat connectivity modelling. For a 10,000×10,000 DEM with complex friction, CPU Dijkstra takes minutes; a GPU implementation could reduce this to seconds.