FIX: support Polygon/MultiPolygon in Voronoi graph builder (#688)

libpysal/graph/_triangulation.py

@@ -431,6 +431,64 @@ def _voronoi(coordinates, coplanar, clip="bounding_box", rook=True):
 #### utilities
 
 
+def _voronoi_polygon(
+    geoms,
+    ids,
+    clip="bounding_box",
+    rook=True,
+    **kwargs,
+):
+    """
+    Compute contiguity weights according to a clipped Voronoi diagram
+    for non-point geometries (Polygon, MultiPolygon).
+
+    Parameters
+    ---------
+    geoms : geopandas.GeoSeries or geopandas.GeoDataFrame
+        Geometries to compute the Voronoi diagram. Accepts Polygon
+        and MultiPolygon geometries. Boundaries are discretised into
+        points and dissolved back using ``voronoi_frames()``.
+    ids : numpy.ndarray
+        ids to use for each sample in geoms.
+    clip : str (default: 'bounding_box')
+        Clipping method passed to ``libpysal.cg.voronoi_frames()``.
+        Options are ``None``, ``'bounding_box'``, ``'convex_hull'``,
+        ``'alpha_shape'``, or a ``shapely.Polygon``.
+    rook : bool, optional
+        Contiguity method. If True, two geometries are considered
+        neighbours if they share at least one edge. If False, they
+        are considered neighbours if they share at least one vertex.
+        By default True.
+    **kwargs
+        Additional keyword arguments passed to ``voronoi_frames()``.
+        Supports ``segment`` (float) to control boundary discretisation
+        and ``shrink`` (float) to shrink geometries before tessellation.
+
+    Returns
+    -------
+    heads : numpy.ndarray
+    tails : numpy.ndarray
+    weights : numpy.ndarray
+    """
+
+    voronoi_kwargs = {k: v for k, v in kwargs.items() if k in ("segment", "shrink")}
+
+    cells = voronoi_frames(
+        geoms, clip=clip, return_input=False, as_gdf=False, **voronoi_kwargs
+    )
+    heads_ix, tails_ix, _ = _vertex_set_intersection(cells, rook=rook)
+
+    valid_ids = set(ids)
+    mask = numpy.isin(heads_ix, list(valid_ids)) & numpy.isin(tails_ix, list(valid_ids))
+    heads_ix = heads_ix[mask]
+    tails_ix = tails_ix[mask]
+
+    heads = heads_ix
+    tails = tails_ix
+    weights = numpy.ones(len(heads_ix), dtype=numpy.int8)
+    return heads, tails, weights
+
+
 @njit
 def _edges_from_simplices(simplices):
     """

libpysal/graph/base.py

@@ -24,7 +24,13 @@
 from ._set_ops import SetOpsMixin
 from ._spatial_lag import _lag_spatial
 from ._summary import GraphSummary
-from ._triangulation import _delaunay, _gabriel, _relative_neighborhood, _voronoi
+from ._triangulation import (
+    _delaunay,
+    _gabriel,
+    _relative_neighborhood,
+    _voronoi,
+    _voronoi_polygon,
+)
 from ._utils import (
     _compute_stats,
     _evaluate_index,
@@ -1377,6 +1383,7 @@ def build_triangulation(
         coplanar="raise",
         taper=True,
         decay=False,
+        **kwargs,
     ):
         """Generate Graph from geometry based on triangulation
 
@@ -1442,6 +1449,9 @@ def build_triangulation(
             or negative) at some very large (possibly infinite) distance.
             Otherwise, kernel functions are treated as proper
             volume-preserving probability distributions.
+        **kwargs: Additional keyword arguments passed to ``voronoi_frames()`` when
+        ``method="voronoi"``. Supports ``segment`` (float) and ``shrink``
+        (float). Ignored for other methods.
 
         Returns
         -------
@@ -1519,15 +1529,24 @@ def build_triangulation(
                 taper=taper,
             )
         elif method == "voronoi":
-            head, tail, weights = _voronoi(
-                data,
-                ids=ids,
-                clip=clip,
-                rook=rook,
-                coplanar=coplanar,
-                decay=decay,
-                taper=taper,
-            )
+            if hasattr(data, "geom_type") and not set(data.geom_type) <= {"Point"}:
+                head, tail, weights = _voronoi_polygon(
+                    data,
+                    ids=ids,
+                    clip=clip,
+                    rook=rook,
+                    **kwargs,
+                )
+            else:
+                head, tail, weights = _voronoi(
+                    data,
+                    ids=ids,
+                    clip=clip,
+                    rook=rook,
+                    coplanar=coplanar,
+                    decay=decay,
+                    taper=taper,
+                )
         else:
             raise ValueError(
                 f"Method '{method}' is not supported. Use one of ['delaunay', "

libpysal/graph/tests/test_triangulation.py

@@ -15,13 +15,15 @@
 import pandas as pd
 import pytest
 import shapely
+from shapely.geometry import MultiPolygon, Polygon
 
 from libpysal.graph._kernel import _kernel_functions
 from libpysal.graph._triangulation import (
     _delaunay,
     _gabriel,
     _relative_neighborhood,
     _voronoi,
+    _voronoi_polygon,
 )
 from libpysal.graph._utils import CoplanarError
 from libpysal.graph.base import Graph
@@ -289,7 +291,7 @@ def test_coplanar_raise_voronoi(stores):
 
 
 @pytest.mark.network
-def test_coplanar_jitter_voronoi(stores, stores_unique):
+def test_coplanar_jitter_voronoi(stores, stores_unique):  # edit
     cp_heads, cp_tails, cp_w = _voronoi(stores, clip=False, coplanar="jitter")
     unique_heads, unique_tails, unique_w = _voronoi(stores_unique, clip=False)
     assert not np.array_equal(cp_heads, unique_heads)
@@ -498,3 +500,84 @@ def test_wrong_resolver(self):
             match="Recieved option coplanar='nonsense'",
         ):
             _delaunay(self.df_int, coplanar="nonsense")
+
+
+def test_voronoi_polygon():
+    polys = [
+        Polygon([(0, 0), (1, 0), (1, 1), (0, 1)]),
+        Polygon([(1, 0), (2, 0), (2, 1), (1, 1)]),
+        Polygon([(0, 1), (1, 1), (1, 2), (0, 2)]),
+        Polygon([(1, 1), (2, 1), (2, 2), (1, 2)]),
+    ]
+    gdf = geopandas.GeoDataFrame(geometry=polys)
+    ids = np.arange(4)
+    heads, tails, weights = _voronoi_polygon(gdf, ids=ids)
+    assert len(np.unique(heads)) == 4
+
+
+def test_voronoi_polygon_kwargs():
+    polys = [
+        Polygon([(0, 0), (1, 0), (1, 1), (0, 1)]),
+        Polygon([(1, 0), (2, 0), (2, 1), (1, 1)]),
+        Polygon([(0, 1), (1, 1), (1, 2), (0, 2)]),
+        Polygon([(1, 1), (2, 1), (2, 2), (1, 2)]),
+    ]
+    gdf = geopandas.GeoDataFrame(geometry=polys)
+    ids = np.arange(4)
+    heads, tails, weights = _voronoi_polygon(gdf, ids=ids, segment=0.5, shrink=0.4)
+    assert len(np.unique(heads)) == 4
+
+
+def test_voronoi_polygon_via_build_triangulation():
+    polys = [
+        Polygon([(0, 0), (1, 0), (1, 1), (0, 1)]),
+        Polygon([(1, 0), (2, 0), (2, 1), (1, 1)]),
+        Polygon([(0, 1), (1, 1), (1, 2), (0, 2)]),
+        Polygon([(1, 1), (2, 1), (2, 2), (1, 2)]),
+    ]
+    gdf = geopandas.GeoDataFrame(geometry=polys)
+    graph = Graph.build_triangulation(gdf, method="voronoi")
+    assert graph.n_nodes == 4
+
+
+def test_voronoi_polygon_multipolygon():
+    polys = [
+        MultiPolygon(
+            [
+                Polygon([(0, 0), (1, 0), (1, 1), (0, 1)]),
+                Polygon([(2, 0), (3, 0), (3, 1), (2, 1)]),
+            ]
+        ),
+        Polygon([(1, 1), (2, 1), (2, 2), (1, 2)]),
+        Polygon([(0, 1), (1, 1), (1, 2), (0, 2)]),
+        Polygon([(1, 0), (2, 0), (2, 1), (1, 1)]),
+    ]
+    gdf = geopandas.GeoDataFrame(geometry=polys)
+    ids = np.arange(4)
+    heads, tails, weights = _voronoi_polygon(gdf, ids=ids)
+    assert len(np.unique(heads)) == 4
+
+
+def test_voronoi_polygon_string_ids():
+    polys = [
+        Polygon([(0, 0), (1, 0), (1, 1), (0, 1)]),
+        Polygon([(1, 0), (2, 0), (2, 1), (1, 1)]),
+        Polygon([(0, 1), (1, 1), (1, 2), (0, 2)]),
+        Polygon([(1, 1), (2, 1), (2, 2), (1, 2)]),
+    ]
+    gdf = geopandas.GeoDataFrame(geometry=polys, index=["a", "b", "c", "d"])
+    ids = np.array(["a", "b", "c", "d"])
+    heads, tails, weights = _voronoi_polygon(gdf, ids=ids)
+    assert set(heads) == {"a", "b", "c", "d"}
+
+
+def test_voronoi_polygon_point_backward_compat():
+    pts = [
+        shapely.Point(0, 0),
+        shapely.Point(1, 0),
+        shapely.Point(0, 1),
+        shapely.Point(1, 1),
+    ]
+    gdf = geopandas.GeoDataFrame(geometry=pts)
+    graph = Graph.build_triangulation(gdf, method="voronoi")
+    assert graph.n_nodes == 4