add union fallback to voronoi_frames (#777)

* add union fallback to voronoi_frames

* add test

* fix the test

Author: martinfleis

libpysal/cg/tests/test_voronoi.py

@@ -220,3 +220,75 @@ def test_deprecated_clip(self, clip):
             match=f"The '{clip}' option for the 'clip' parameter is deprecated",
         ):
             voronoi_frames(self.points2, clip=clip)
+
+    # GH776
+    def test_union_fallback(self):
+        data = {
+            "type": "FeatureCollection",
+            "features": [
+                {
+                    "id": "0",
+                    "type": "Feature",
+                    "properties": {},
+                    "geometry": {
+                        "type": "Polygon",
+                        "coordinates": (
+                            (
+                                (575361.23, 115796.88),
+                                (575364.42, 115767.29),
+                                (575343.54, 115765.04),
+                                (575341.61, 115782.98),
+                                (575347.84, 115783.65),
+                                (575348.94, 115773.43),
+                                (575356.99, 115774.3),
+                                (575354.63, 115796.16),
+                                (575361.23, 115796.88),
+                            ),
+                        ),
+                    },
+                    "bbox": (575341.61, 115765.04, 575364.42, 115796.88),
+                },
+                {
+                    "id": "1",
+                    "type": "Feature",
+                    "properties": {},
+                    "geometry": {
+                        "type": "Polygon",
+                        "coordinates": (
+                            (
+                                (575365.47, 115821.98),
+                                (575375.28, 115823.68),
+                                (575379.67, 115798.26),
+                                (575369.86, 115796.56),
+                                (575365.47, 115821.98),
+                            ),
+                        ),
+                    },
+                    "bbox": (575365.47, 115796.56, 575379.67, 115823.68),
+                },
+                {
+                    "id": "2",
+                    "type": "Feature",
+                    "properties": {},
+                    "geometry": {
+                        "type": "Polygon",
+                        "coordinates": (
+                            (
+                                (575344.39, 115825.85),
+                                (575343.47, 115831.29),
+                                (575348.91, 115832.21),
+                                (575349.83, 115826.77),
+                                (575344.39, 115825.85),
+                            ),
+                        ),
+                    },
+                    "bbox": (575343.47, 115825.85, 575349.83, 115832.21),
+                },
+            ],
+            "bbox": (575341.61, 115765.04, 575379.67, 115832.21),
+        }
+        geometry = gpd.GeoDataFrame.from_features(data)
+        vf = voronoi_frames(
+            geometry, shrink=0.4, segment=0.5, return_input=False, as_gdf=False
+        )
+        assert vf.shape == (3,)

libpysal/cg/voronoi.py

@@ -13,6 +13,7 @@
 import shapely
 from packaging.version import Version
 from scipy.spatial import Voronoi
+from shapely.errors import GEOSException
 
 __author__ = "Serge Rey <[email protected]>"
 
@@ -455,7 +456,7 @@ def voronoi_frames(
         polygons = (
             polygons.iloc[ids_polygons]
             .groupby(objects.index.take(ids_objects))
-            .agg(shapely.coverage_union_all)
+            .agg(_union_with_fallback)
         )
         if geometry.crs is not None:
             polygons = polygons.set_crs(geometry.crs)
@@ -497,6 +498,20 @@ def voronoi_frames(
     return polygons
 
 
+def _union_with_fallback(arr):
+    """
+    Coverage union is finnicky with floating point precision and tends to occasionally
+    raise an error from within GEOS we have no control over. It is not a data issue
+    typically. Falling back to unary union if that happens.
+    """
+    try:
+        r = shapely.coverage_union_all(arr)
+    except GEOSException:
+        r = shapely.union_all(arr)
+
+    return r
+
+
 def _get_limit(points, clip):
     if isinstance(clip, shapely.Geometry):
         return clip