Polygonize: GPU CCL, promotion, dispatcher, GeoJSON, benchmarks (#916)

* Add dask, cupy, and dask+cupy backends for polygonize, update README

Polygonize now supports all four backends. The cupy and dask+cupy
backends use a hybrid GPU/CPU approach: GPU arrays are transferred to
CPU for boundary tracing, which is inherently sequential. The dask
backend chunks independently, merging cross-boundary polygons.

* Promote polygonize out of experimental, add GPU CCL, dispatcher, GeoJSON

- Move polygonize.py from xrspatial/experimental/ to xrspatial/ and update
  all imports (tests, benchmarks, README, __init__, accessor)
- Replace isinstance dispatch chain with ArrayTypeFunctionMapping
- Add GPU-accelerated CCL via cupyx.scipy.ndimage.label with vectorized
  region assignment and raster-scan renumbering for _scan compatibility
- Refactor _scan to accept pre-computed regions array
- Add GeoJSON FeatureCollection return_type ("geojson")
- Add CCL vs tracing phase benchmarks (PolygonizeCCL, PolygonizeTracing)
- Add end-to-end CuPy benchmark (PolygonizeCuPy)

Author: brendancol

README.md

@@ -213,7 +213,7 @@ In the GIS world, rasters are used for representing continuous phenomena (e.g. e
 
 | Name | Description | NumPy xr.DataArray | Dask xr.DataArray | CuPy GPU xr.DataArray | Dask GPU xr.DataArray |
 |:-----|:------------|:------------------:|:-----------------:|:---------------------:|:---------------------:|
-| [Polygonize](xrspatial/experimental/polygonize.py) | Converts contiguous regions of equal value into vector polygons | ✅️ | | | |
+| [Polygonize](xrspatial/polygonize.py) | Converts contiguous regions of equal value into vector polygons | ✅️ | ✅️ | ✅️ | ✅️ |
 
 --------
 

benchmarks/benchmarks/polygonize.py

@@ -1,7 +1,17 @@
 import numpy as np
 import xarray as xr
 
-from xrspatial.experimental import polygonize
+from xrspatial import polygonize
+from xrspatial.polygonize import (
+    _calculate_regions, _scan,
+)
+
+try:
+    import cupy
+    from xrspatial.polygonize import _calculate_regions_cupy
+    _has_cupy = True
+except ImportError:
+    _has_cupy = False
 
 
 class Polygonize:
@@ -40,3 +50,75 @@ def time_polygonize(self, nx, ret):
                 gpd.GeoDataFrame({"DN": values, "geometry": shapes})
         else:
             polygonize(self.raster, mask=self.mask, return_type=ret)
+
+
+class PolygonizeCuPy:
+    params = ([100, 300, 1000],)
+    param_names = ("nx",)
+
+    def setup(self, nx):
+        if not _has_cupy:
+            raise NotImplementedError("CuPy not available")
+        ny = nx // 2
+        rng = np.random.default_rng(9461713)
+        raster = rng.integers(low=0, high=4, size=(ny, nx), dtype=np.int32)
+        mask = rng.uniform(0, 1, size=(ny, nx)) < 0.9
+        self.raster = xr.DataArray(cupy.asarray(raster))
+        self.mask = xr.DataArray(cupy.asarray(mask))
+
+    def time_polygonize(self, nx):
+        polygonize(self.raster, mask=self.mask)
+
+
+class PolygonizeCCL:
+    """Benchmark connected-component labeling phase."""
+    params = ([100, 300, 1000], ["numpy", "cupy"])
+    param_names = ("nx", "backend")
+
+    def setup(self, nx, backend):
+        if backend == "cupy" and not _has_cupy:
+            raise NotImplementedError("CuPy not available")
+        ny = nx // 2
+        rng = np.random.default_rng(9461713)
+        raster = rng.integers(low=0, high=4, size=(ny, nx), dtype=np.int32)
+        mask = (rng.uniform(0, 1, size=(ny, nx)) < 0.9)
+        if backend == "cupy":
+            self.data = cupy.asarray(raster)
+            self.mask = cupy.asarray(mask)
+        else:
+            self.data = raster
+            self.mask = mask
+        self.ny, self.nx = ny, nx
+        self.backend = backend
+
+    def time_calculate_regions(self, nx, backend):
+        if backend == "cupy":
+            _calculate_regions_cupy(self.data, self.mask, False)
+            cupy.cuda.Device().synchronize()
+        else:
+            _calculate_regions(
+                self.data.ravel(), self.mask.ravel(), False, self.nx, self.ny)
+
+
+class PolygonizeTracing:
+    """Benchmark boundary tracing phase (always CPU)."""
+    params = ([100, 300, 1000],)
+    param_names = ("nx",)
+
+    def setup(self, nx):
+        ny = nx // 2
+        rng = np.random.default_rng(9461713)
+        raster = rng.integers(low=0, high=4, size=(ny, nx), dtype=np.int32)
+        mask = (rng.uniform(0, 1, size=(ny, nx)) < 0.9)
+        values_flat = raster.ravel()
+        mask_flat = mask.ravel()
+        self.regions = _calculate_regions(
+            values_flat, mask_flat, False, nx, ny)
+        self.values = values_flat
+        self.mask = mask_flat
+        self.nx = nx
+        self.ny = ny
+
+    def time_scan(self, nx):
+        _scan(self.regions, self.values, self.mask, False, None,
+              self.nx, self.ny)

xrspatial/__init__.py

@@ -33,6 +33,7 @@
 from xrspatial.proximity import proximity  # noqa
 from xrspatial.slope import slope  # noqa
 from xrspatial.terrain import generate_terrain  # noqa
+from xrspatial.polygonize import polygonize  # noqa
 from xrspatial.viewshed import viewshed  # noqa
 from xrspatial.zonal import apply as zonal_apply  # noqa
 from xrspatial.zonal import crop  # noqa

xrspatial/accessor.py

@@ -145,6 +145,12 @@ def mahalanobis(self, other_bands, **kwargs):
         from .mahalanobis import mahalanobis
         return mahalanobis([self._obj] + list(other_bands), **kwargs)
 
+    # ---- Raster to vector ----
+
+    def polygonize(self, **kwargs):
+        from .polygonize import polygonize
+        return polygonize(self._obj, **kwargs)
+
     # ---- Multispectral (self = NIR band) ----
 
     def ndvi(self, red_agg, **kwargs):

xrspatial/experimental/__init__.py

@@ -1,2 +1 @@
 from .min_observable_height import min_observable_height  # noqa
-from .polygonize import polygonize  # noqa