Describe the bug
In xrspatial/proximity.py, the CPU brute-force distance helper _distance casts its result to np.float32 before returning (around line 261). The brute-force search then compares that float32-rounded distance against max_distance via best_dist <= max_distance. The CUDA kernel _proximity_cuda_kernel and its device helpers compute and compare in float64.
When max_distance lands between the float32-rounded distance and the true float64 distance of a target, the backends disagree: the CPU path treats the target as in range and returns a finite value, while the GPU path treats it as out of range and returns NaN (or the reverse).
This hits every path that runs the brute-force search: GREAT_CIRCLE proximity/allocation/direction, and EUCLIDEAN/MANHATTAN allocation/direction. EUCLIDEAN/MANHATTAN proximity on numpy goes through the cKDTree path, which bounds differently.
Reproducer
A single-target lat/lon raster, GREAT_CIRCLE proximity, with max_distance set between the float32 and float64 distance of one column:
import numpy as np, xarray as xr, cupy as cp
from xrspatial import proximity
from xrspatial.proximity import great_circle_distance
lon = np.linspace(0, 39, 40); lat = np.array([0.0])
data = np.zeros((1, 40)); data[0, 0] = 1.0
j = 20
d64 = great_circle_distance(lon[0], lon[j], lat[0], lat[0])
md = (float(np.float32(d64)) + d64) / 2.0 # between f32 and f64 distance
r_np = proximity(xr.DataArray(data, dims=['lat','lon'], coords={'lon':lon,'lat':lat}),
x='lon', y='lat', distance_metric='GREAT_CIRCLE', max_distance=md)
r_cp = proximity(xr.DataArray(cp.asarray(data), dims=['lat','lon'], coords={'lon':lon,'lat':lat}),
x='lon', y='lat', distance_metric='GREAT_CIRCLE', max_distance=md)
print(r_np.data[0, j]) # 2226389.75 (finite)
print(r_cp.data.get()[0, j]) # nannumpy and dask+numpy return the finite distance; cupy and dask+cupy return NaN for the same input.
Expected behavior
All four backends should agree on whether a target is within max_distance. The comparison should use the same precision everywhere.
Severity
Medium. The divergence is silent, with no error raised, but it only surfaces when max_distance is set with sub-ulp precision relative to a target's distance. That is an adversarial input rather than typical usage.
Found by /sweep-accuracy on the proximity module.
Describe the bug
In
xrspatial/proximity.py, the CPU brute-force distance helper_distancecasts its result tonp.float32before returning (around line 261). The brute-force search then compares that float32-rounded distance againstmax_distanceviabest_dist <= max_distance. The CUDA kernel_proximity_cuda_kerneland its device helpers compute and compare in float64.When
max_distancelands between the float32-rounded distance and the true float64 distance of a target, the backends disagree: the CPU path treats the target as in range and returns a finite value, while the GPU path treats it as out of range and returns NaN (or the reverse).This hits every path that runs the brute-force search: GREAT_CIRCLE proximity/allocation/direction, and EUCLIDEAN/MANHATTAN allocation/direction. EUCLIDEAN/MANHATTAN proximity on numpy goes through the cKDTree path, which bounds differently.
Reproducer
A single-target lat/lon raster, GREAT_CIRCLE proximity, with
max_distanceset between the float32 and float64 distance of one column:numpy and dask+numpy return the finite distance; cupy and dask+cupy return NaN for the same input.
Expected behavior
All four backends should agree on whether a target is within
max_distance. The comparison should use the same precision everywhere.Severity
Medium. The divergence is silent, with no error raised, but it only surfaces when
max_distanceis set with sub-ulp precision relative to a target's distance. That is an adversarial input rather than typical usage.Found by /sweep-accuracy on the proximity module.