proximity: match GPU max_distance precision to the CPU brute-force (#3389)

.claude/sweep-accuracy-state.csv

@@ -27,7 +27,7 @@ normalize,2026-05-01,,,,rescale and standardize across all 4 backends. NaN/inf f
 perlin,2026-04-10T12:00:00Z,,,,Improved Perlin noise implementation correct. Fade/gradient functions verified. Backend-consistent. Continuous at cell boundaries.
 polygon_clip,2026-06-10,3186,HIGH,5,"Cat5 backend inconsistency: dask+cupy clip_polygon rasterizes the mask with a uniform chunk size from the raster's first chunk, then feeds raster+mask to da.map_blocks (positional block pairing). Non-uniform raster chunks gave the mask a different block layout -> IndexError/ValueError (or silent mis-stamp). Repro (8,6) rechunk ((3,5),(6,)) on dask+cupy raised ValueError Shapes do not align; dask+numpy was fine via xarray.where rechunk. Fix #3186/PR: rechunk cond to raster.data.chunks[-2:] before map_blocks; added non-uniform regression tests for dask+numpy and dask+cupy. use_cuda->gpu migration in that branch was already landed by #3089/#3122. CUDA available; cupy+dask+cupy verified, 25 tests pass. Cats 1-4 clean: numpy path uses raster.where, cupy path operates on raw arrays, NaN inputs preserved, no neighborhood ops/curvature. Prior fix #1197/#1200 (crop+all_touched) merged and unrelated."
 polygonize,2026-05-29,2606,HIGH,5,"Cat 5 HIGH: dask connectivity=8 cross-chunk merge filled diagonal notch where same-value regions meet only at a corner across a chunk boundary; total area exceeded raster. Hole ring was dropped because containment tested hole[0] (on exterior at pinch). Fixed via _ring_interior_point in PR for #2606. numpy, dask+numpy, dask+cupy area parity now holds; 4-conn was already correct. cupy + dask+cupy paths validated on GPU host. Other cats clean: NaN masked on numpy/cupy float paths (tested), _is_close handles +/-inf via exact-equality short-circuit, atol/rtol/simplify_tolerance reject NaN/inf, integer GPU CCL matches numpy."
-proximity,2026-06-09,3108,HIGH,4;5,"Cat5/Cat4: bounded GREAT_CIRCLE dask (numpy+cupy) missed targets across the +/-180 antimeridian seam: _halo_depth sized x-halo as linear parallel-arc sum, but haversine is periodic in lon and chords shorten near poles, so array-space adjacency is no lower bound on spherical distance; numpy/cupy (brute force) found the wrap target (~111 km), dask returned NaN. Fixed in #3108 via chord bound 2R asin(cos(lat_max)|sin(dlon/2)|) + x-axis fold when seam/180-deg chord within max_distance (covers over-pole too). CUDA host: cupy + dask+cupy executed, 417+ tests pass. Cat1-3 clean (float32 output documented; NaN via isfinite consistent; bounds guards correct; tie-break unified in #2881). LOW (not fixed): great_circle_distance uses WGS84 equatorial radius 6378137 as sphere radius (~0.1% vs mean-radius convention) but documented and exposed as param."
+proximity,2026-06-18,3389,MEDIUM,1;5,"Cat1/Cat5 (MEDIUM, fixed #3389): CPU brute-force _distance casts to float32 before best_dist<=max_distance, but the CUDA kernel compared the float64 distance; a max_distance in a target's float32 ulp gap qualified it on numpy/dask+numpy yet rejected it on cupy/dask+cupy (GREAT_CIRCLE proximity + EUCLIDEAN/MANHATTAN allocation/direction). Fixed by rounding best_dist to float32 in the kernel before the range test. CUDA host: 518 tests pass incl cupy + dask+cupy. Out of scope (separate pre-existing edge, predates fix): bounded dask EUCLIDEAN allocation can still NaN a target exactly at a float32-ulp max_distance via the halo-sizing path. Prior #3108 fix (antimeridian halo) still verified. LOW (unfixed): great_circle_distance uses WGS84 equatorial radius 6378137 as sphere radius (~0.1% vs mean-radius), but documented and param-exposed."
 rasterize,2026-06-18,3384,HIGH,1;5,dask all_touched polygon-boundary supercover walk re-extracted per tile in tile-local float pixel coords; floor() tie on an on-grid boundary segment diverged from eager (Cat1 precision -> Cat5 backend split). Fixed by extracting boundary float segments in the global grid frame and shifting by integer tile offset (#3384). Verified numpy/cupy/dask+numpy/dask+cupy.
 reproject,2026-06-12,3274,HIGH,1;4,"3 confirmed bugs, all kernel-vs-PROJ parity: #3274 HIGH LAEA inverse spurious /rq (2.6 km err for 3035) + _authalic_apa inverse series wrong (4.8 m in AEA/CEA inverses; PROJ 3-term = 1.6 mm), CPU+CUDA kernels both; #3275 HIGH _is_wgs84_compatible_ellipsoid passes R-defined spheres (MODIS sinusoidal 18.9 km err) and _aea_params/_cea_params lack the guard entirely (23.8 km on spherical aea/cea); #3276 MEDIUM itrf helmert scale 1e-9 but PROJ +s is ppm (1e-6), ~23 mm err. Verified clean: merc/emerc/UTM/tmerc/LCC/polar stere (incl lat_ts akm1) forward+inverse <=1e-5 m vs pyproj; resampling kernels NaN handling and GDAL renorm match across numpy/cupy/dask (CUDA run, gpu-vs-cpu 1.3e-7); dask footprint chunk-skip bbox is a superset in all probed cases (no holes). LOW (documented only): _source_footprint_in_target probe array typo uses x-midpoint mx as a latitude in last 3 ys entries (bbox superset, correctness unaffected)."
 resample,2026-05-29,2610,HIGH,3;5,"dask interp (nearest/bilinear) overlap depth=1 too small on downsample; block-centered source coord landed past chunk, map_coordinates clamped to edge -> wrong seam rows. Fixed PR #2627 via per-axis _downsample_radius. cupy+dask+cupy verified."

xrspatial/proximity.py

@@ -634,6 +634,12 @@ def _proximity_cuda_kernel(target_xs, target_ys, target_vals, n_targets,
             best_dist = d
             best_idx = k
 
+    # Round the winning distance to float32 before the range test so the
+    # in-range decision matches the CPU brute-force path, where _distance
+    # returns np.float32(d). Comparing the float64 distance here let a target
+    # whose distance rounds down across a float32 ulp pass on the CPU but fail
+    # on the GPU (or the reverse) when max_distance sits in that ulp gap.
+    best_dist = np.float32(best_dist)
     if best_idx >= 0 and best_dist <= max_distance:
         if process_mode == PROXIMITY:
             out[iy, ix] = best_dist

xrspatial/tests/test_proximity.py

@@ -2499,3 +2499,75 @@ def test_proximity_numpy_no_scipy_fallback_is_exact():
     expected = _exact_planar_proximity(
         raster.data, raster['x'].data, raster['y'].data, 'EUCLIDEAN')
     np.testing.assert_allclose(result.data, expected, rtol=1e-5, atol=1e-7)
+
+
+# ---------------------------------------------------------------------------
+# Issue #3389: the max_distance comparison must use the same precision on the
+# CPU brute-force and the CUDA kernel. _distance casts to float32 before the
+# CPU compares best_dist <= max_distance, while the kernel used to compare the
+# float64 distance. A max_distance landing inside the float32 ulp gap of a
+# target's distance then qualified the target on numpy but rejected it on cupy
+# (or the reverse).
+# ---------------------------------------------------------------------------
+
+
+def test_max_distance_float32_boundary_gpu_matches_cpu():
+    """GREAT_CIRCLE proximity at a float32-ulp max_distance agrees CPU/GPU."""
+    if not has_cuda_and_cupy():
+        pytest.skip("Requires CUDA and CuPy")
+    import cupy
+
+    lon = np.linspace(0, 39, 40)
+    lat = np.array([0.0])
+    data = np.zeros((1, 40), dtype=np.float64)
+    data[0, 0] = 1.0
+    j = 20
+    d64 = great_circle_distance(lon[0], lon[j], lat[0], lat[0])
+    # max_distance between the float32-rounded and the true float64 distance.
+    md = (float(np.float32(d64)) + d64) / 2.0
+    assert float(np.float32(d64)) < md < d64, "needs a real float32 ulp gap"
+
+    kwargs = dict(x='lon', y='lat', distance_metric='GREAT_CIRCLE',
+                  max_distance=md)
+    numpy_raster = xr.DataArray(
+        data, dims=['lat', 'lon'], coords={'lon': lon, 'lat': lat})
+    numpy_val = proximity(numpy_raster, **kwargs).data[0, j]
+
+    cupy_raster = xr.DataArray(
+        cupy.asarray(data), dims=['lat', 'lon'],
+        coords={'lon': lon, 'lat': lat})
+    cupy_val = proximity(cupy_raster, **kwargs).data.get()[0, j]
+
+    # Both keep the target: its float32 distance is <= max_distance.
+    assert np.isnan(numpy_val) == np.isnan(cupy_val)
+    np.testing.assert_allclose(numpy_val, cupy_val, rtol=1e-6)
+
+
+def test_max_distance_float32_boundary_allocation_gpu_matches_cpu():
+    """EUCLIDEAN allocation (brute-force on both CPU and GPU) agrees at the
+    float32 ulp boundary."""
+    if not has_cuda_and_cupy():
+        pytest.skip("Requires CUDA and CuPy")
+    import cupy
+
+    lon = np.array([0.0, 0.1, 0.7])
+    lat = np.array([0.0])
+    data = np.zeros((1, 3), dtype=np.float64)
+    data[0, 0] = 9.0
+    j = 2
+    d64 = euclidean_distance(lon[0], lon[j], lat[0], lat[0])
+    md = (float(np.float32(d64)) + d64) / 2.0
+    assert float(np.float32(d64)) < md < d64, "needs a real float32 ulp gap"
+
+    kwargs = dict(x='lon', y='lat', max_distance=md)
+    numpy_raster = xr.DataArray(
+        data, dims=['lat', 'lon'], coords={'lon': lon, 'lat': lat})
+    numpy_val = allocation(numpy_raster, **kwargs).data[0, j]
+
+    cupy_raster = xr.DataArray(
+        cupy.asarray(data), dims=['lat', 'lon'],
+        coords={'lon': lon, 'lat': lat})
+    cupy_val = allocation(cupy_raster, **kwargs).data.get()[0, j]
+
+    assert np.isnan(numpy_val) == np.isnan(cupy_val)
+    np.testing.assert_allclose(numpy_val, cupy_val, rtol=1e-6)