polygonize: cupy backend honours float tolerance from numba path (#2151) (#2153)

* polygonize: cupy backend honours numba _is_close tolerance for floats (#2151)

The cupy backend binned float pixels into regions by exact value while
the numba CPU path (used by the numpy, dask and dask+cupy backends)
compares neighbouring pixels with `_is_close` (atol=1e-8, rtol=1e-5).
On float rasters with near-equal values this produced different polygon
sets across backends.

For float dtypes, group sorted unique values in the cupy backend using
the same tolerance and label each group as a single region.  Integer
dtypes keep the existing exact-equality path.

Adds a regression test asserting that the cupy backend yields the same
polygon count and per-value areas as the numpy backend on a float
raster with values 1.0 and 1.000001.

Also records this finding in the accuracy sweep state CSV.

* polygonize: route cupy float path through CPU CCL for true _is_close parity

The first pass at #2151 added value-only greedy chaining of unique float
values in _calculate_regions_cupy.  That fixed the reported pattern
(adjacent near-equal pixels) but still diverged from the CPU path when
three or more transitively-close values were present in the raster
without intermediate adjacent pixels.  Example:

    [1.0,        1.000018]
    [1.000009,   3.0      ]

CPU spatial CCL yields three regions: 1.0 + 1.000009 (close, adjacent),
1.000018 alone (not close to either neighbour), and 3.0 alone.  The
value-grouping heuristic chained 1.0 -> 1.000009 -> 1.000018 and labeled
the connected mask as a single region, returning two regions on cupy.

Route float dtypes through _polygonize_numpy in _polygonize_cupy so the
numba _is_close predicate is applied to spatially adjacent pixels, the
only place CPU and GPU semantics can agree without re-implementing CCL.
The data is already transferred to host for boundary tracing, so the
fall-back cost is one extra _calculate_regions call.

Drop _group_float_values_by_tolerance and the float branch of
_calculate_regions_cupy as dead code.  Extend the regression test to
parametrize over the original adjacent pattern and the transitivity
edge case.

Author: brendancol

.claude/sweep-accuracy-state.csv

@@ -23,7 +23,7 @@ multispectral,2026-03-30T14:00:00Z,1094,,,
 normalize,2026-05-01,,,,rescale and standardize across all 4 backends. NaN/inf filtered via isfinite mask before min/max/mean/std. Constant input handled (range=0 -> new_min; std=0 -> 0.0). Output dtype float64 consistently. Backend parity covered by test_matches_numpy. No accuracy issues found.
 perlin,2026-04-10T12:00:00Z,,,,Improved Perlin noise implementation correct. Fade/gradient functions verified. Backend-consistent. Continuous at cell boundaries.
 polygon_clip,2026-04-13T12:00:00Z,1197,,,crop=True + all_touched=True drops boundary pixels. Fix in PR #1200.
-polygonize,2026-04-13T12:00:00Z,1190,,,NaN pixels not masked in numpy/dask backends. Fix in PR #1194.
+polygonize,2026-05-19,2151,HIGH,5,cupy backend used exact-equality vs numba _is_close tolerance for float dtypes; fixed via tolerance-aware value grouping
 proximity,2026-03-30T15:00:00Z,,,,Direction >= boundary fragile but works due to truncated constant. Float32 truncation is design choice. No wrong-results bugs found.
 reproject,2026-05-17,2025,HIGH,5,"HIGH: _apply_vertical_shift crashed on dask (no fancy-index setitem), cupy (Numba JIT can't take cupy), and 3D multi-band (broadcast mismatch). Fixed by routing dask through map_blocks, round-tripping cupy via host, and looping per band. Previously documented MEDIUM (GPU dispatcher missing NAD27/Helmert wrapper) and MEDIUM (legacy _resample_cupy cubic NaN threshold 1e-6 diverges from CPU bilinear-fallback) still stand. Identity-CRS cupy reproject returns zeros - separate pre-existing bug, not in scope."
 resample,2026-04-14T12:00:00Z,1202,,,Interpolation used edge-aligned coords instead of block-centered. Fix in PR #1204.

xrspatial/polygonize.py

@@ -594,10 +594,16 @@ def _polygonize_numpy(
 
 
 def _calculate_regions_cupy(data, mask_data, connectivity_8):
-    """CuPy GPU backend for connected-component labeling.
+    """CuPy GPU backend for connected-component labeling on integer data.
 
     Uses cupyx.scipy.ndimage.label per unique value to produce a regions
     array compatible with _scan.  Returns a cupy uint32 2D array.
+
+    Only used for integer dtypes; float dtypes route through the CPU
+    ``_calculate_regions`` so that connected components honour the numba
+    ``_is_close`` predicate (``atol=1e-8``, ``rtol=1e-5``) exactly.
+    A purely value-based grouping on the GPU cannot reproduce CPU spatial
+    CCL when transitively-close values are not spatially adjacent.
     """
     import cupy as cp
     from cupyx.scipy.ndimage import label as cp_label
@@ -609,14 +615,10 @@ def _calculate_regions_cupy(data, mask_data, connectivity_8):
 
     regions = cp.zeros(data.shape, dtype=cp.uint32)
 
-    # Build valid mask (unmask + handle float NaN).
-    is_float = cp.issubdtype(data.dtype, cp.floating)
     if mask_data is not None:
         valid = cp.asarray(mask_data, dtype=bool)
-        if is_float:
-            valid &= ~cp.isnan(data)
     else:
-        valid = ~cp.isnan(data) if is_float else None
+        valid = None
 
     unique_vals = data[valid] if valid is not None else data.ravel()
     unique_vals = cp.unique(unique_vals)
@@ -669,9 +671,19 @@ def _renumber_regions(regions, nx, ny):
 
 
 def _polygonize_cupy(data, mask_data, connectivity_8, transform):
-    """Hybrid GPU/CPU: GPU CCL, CPU boundary tracing."""
+    """Hybrid GPU/CPU: GPU CCL for integer data, CPU CCL for float data,
+    CPU boundary tracing in either case.
+
+    Float dtypes route through ``_polygonize_numpy`` so that connected
+    components honour the numba ``_is_close`` tolerance (``atol=1e-8``,
+    ``rtol=1e-5``) for spatially adjacent pixels (#2151).  GPU CCL is a
+    per-value labeling, which cannot reproduce spatial tolerance-aware
+    CCL when transitively-close values are not spatially adjacent.
+    """
     np_data = cupy.asnumpy(data)
     np_mask = cupy.asnumpy(mask_data) if mask_data is not None else None
+    if np.issubdtype(np_data.dtype, np.floating):
+        return _polygonize_numpy(np_data, np_mask, connectivity_8, transform)
     ny, nx = np_data.shape
     if nx == 1:
         # Edge case: fall back to full numpy path (pads array).

xrspatial/tests/test_polygonize.py

@@ -459,6 +459,62 @@ def test_polygonize_dask_cupy_matches_numpy(chunks):
                         err_msg=f"Area mismatch for value {val}")
 
 
+@cuda_and_cupy_available
+@pytest.mark.parametrize(
+    "data",
+    [
+        # Adjacent near-equal values, the originally reported pattern.
+        np.array([
+            [1.0, 1.000001, 2.0],
+            [1.000001, 1.0, 2.0],
+            [3.0, 3.0, 2.0],
+        ], dtype=np.float64),
+        # Transitivity edge case: three values 1.0, 1.000009, 1.000018 where
+        # 1.0 and 1.000018 are NOT pairwise close by ``_is_close`` but a
+        # value-only grouping heuristic would chain them transitively.  The
+        # CPU spatial CCL keeps them in separate regions when they are not
+        # bridged by an adjacent intermediate pixel; the cupy backend must
+        # match.
+        np.array([
+            [1.0, 1.000018],
+            [1.000009, 3.0],
+        ], dtype=np.float64),
+    ],
+    ids=["adjacent_near_equal", "transitive_non_adjacent"],
+)
+def test_polygonize_cupy_float_tolerance_matches_numpy_2151(data):
+    """CuPy backend must use the same float tolerance as the numpy/numba
+    path for grouping pixels into regions (#2151).
+
+    Before the fix, the cupy backend binned pixels by exact value while
+    the numpy/numba backend used ``_is_close`` (atol=1e-8, rtol=1e-5),
+    yielding different polygon counts on float rasters with near-equal
+    values.
+    """
+    import cupy
+
+    raster_np = xr.DataArray(data)
+    vals_np, polys_np = polygonize(raster_np, connectivity=4)
+
+    raster_cp = xr.DataArray(cupy.asarray(data))
+    vals_cp, polys_cp = polygonize(raster_cp, connectivity=4)
+
+    # Both backends should return the same number of polygons.
+    assert len(polys_np) == len(polys_cp), (
+        f"polygon count differs: numpy={len(polys_np)} cupy={len(polys_cp)}; "
+        f"numpy values={vals_np}, cupy values={vals_cp}"
+    )
+    # Per-value total area must agree across backends.
+    areas_np = _area_by_value(vals_np, polys_np)
+    areas_cp = _area_by_value(vals_cp, polys_cp)
+    assert set(areas_np.keys()) == set(areas_cp.keys()), (
+        f"value sets differ: numpy={set(areas_np.keys())} "
+        f"cupy={set(areas_cp.keys())}"
+    )
+    for v, a in areas_np.items():
+        assert_allclose(areas_cp[v], a, atol=1e-10)
+
+
 # --- Performance-related regression tests (#1008) ---
 
 def test_polygonize_1008_jit_merge_helpers():