Repository: Toblerity/Shapely
Pull request #881 (2020-03-31 03:14:15)

Add GEOSVoronoiDiagram to shapely
Per #833. 

I thought it made sense to keep Delauney / Voronoi close together in the code since they're similar in use and right next to each other in the CAPI, but if you disagree I'm happy to move it to wherever.

I couldn't nail down for sure exactly when VoronoiDiagram was added to the CAPI. [This](https://github.com/libgeos/geos/blob/master/NEWS#L167) entry was the best I could find indicating release version, so I went with 3.5.0.

coveralls: 
[![Coverage Status](https://coveralls.io/builds/29729351/badge)](https://coveralls.io/builds/29729351)

Coverage remained the same at ?% when pulling **8eecd46307a809ec6aa8f5ac2fa872a883e5e522 on tomplex:voronoi-diagram** into **91c70ba508a4030e20f2ef22e6e6acc6c1a69d82 on Toblerity:maint-1.7**.

jorisvandenbossche: The current return value is a list of geometries. Could it also be considered to return the actual GeometryCollection (as it is returned by GEOS)? The user could then still get the list by doing `result.geoms` themselves.

I assume that for users that pass a single geometry, the current API returning a list is easier (and also more consisitent with the existing `triangulate`, which is a good reason of course ..). But returning a geometry collection would make it easier to have it compatible with the PyGEOS version.
tomplex: @jorisvandenbossche, though I personally agree that a `GeometryCollection` feels like a better option since that's what being returned natively by GEOS, I think one of the best parts of Shapely is the consistency and therefore intuitiveness of the API, so I went with a list of geometries. I'd be happy to discuss more if you feel strongly, especially because it seems like there will be a closer union of Shapely and PyGeos in the future, so it makes sense to take such things into account.
tomplex: @sgillies should I rebase this onto master?
sgillies: @tomplex yes, please. 

Returning a GeometryCollection would be more efficient and I'm in favor :+1: 

Patch:
diff --git a/docs/code/voronoi_diagram.py b/docs/code/voronoi_diagram.py
new file mode 100644
index 0000000..09d2c41
--- /dev/null
+++ b/docs/code/voronoi_diagram.py
@@ -0,0 +1,24 @@
+from shapely.geometry import MultiPoint
+from shapely.ops import voronoi_diagram
+
+from matplotlib import pyplot
+from descartes.patch import PolygonPatch
+from figures import SIZE, BLUE, GRAY, set_limits
+
+points = MultiPoint([(0, 0), (1, 1), (0, 2), (2, 2), (3, 1), (1, 0)])
+regions = voronoi_diagram(points)
+
+fig = pyplot.figure(1, figsize=SIZE, dpi=90)
+fig.set_frameon(True)
+ax = fig.add_subplot(111)
+
+for region in regions:
+    patch = PolygonPatch(region, facecolor=BLUE, edgecolor=BLUE, alpha=0.5, zorder=2)
+    ax.add_patch(patch)
+
+for point in points:
+    pyplot.plot(point.x, point.y, 'o', color=GRAY)
+
+set_limits(ax, -1, 4, -1, 3)
+
+pyplot.show()
diff --git a/docs/manual.rst b/docs/manual.rst
index c83db2d..8c18dd2 100644
--- a/docs/manual.rst
+++ b/docs/manual.rst
@@ -2195,6 +2195,51 @@ Delaunay triangulation from a collection of points.
    'POLYGON ((3 1, 1 1, 1 0, 3 1))',
    'POLYGON ((1 0, 1 1, 0 0, 1 0))']
 
+Voronoi Diagram
+---------------
+
+The :func:`~shapely.ops.voronoi_diagram` function in `shapely.ops` constructs a
+Voronoi diagram from a collection points, or the vertices of any geometry.
+
+.. plot:: code/voronoi_diagram.py
+
+.. function:: shapely.ops.voronoi_diagram(geom, envelope=None, tolerance=0.0, edges=False)
+
+   Constructs a Voronoi diagram from the vertices of the input geometry.
+
+   The source may be any geometry type. All vertices of the geometry will be
+   used as the input points to the diagram.
+
+   The `envelope` keyword argument provides an envelope to use to clip the
+   resulting diagram. If `None`, it will be calculated automatically.
+   The diagram will be clipped to the *larger* of the provided envelope
+   or an envelope surrounding the sites.
+
+   The `tolerance` keyword argument sets the snapping tolerance used to improve
+   the robustness of the computation. A tolerance of 0.0 specifies
+   that no snapping will take place. The tolerance `argument` can be
+   finicky and is known to cause the algorithm to fail in several cases.
+   If you're using `tolerance` and getting a failure, try removing it.
+   The test cases in `tests/test_voronoi_diagram.py` show more details.
+
+   If the `edges` keyword argument is `False` a list of `Polygon`s
+   will be returned. Otherwise a list of `LineString` edges is returned.
+
+
+.. code-block:: pycon
+
+  >>> from shapely.ops import voronoi_diagram
+  >>> points = MultiPoint([(0, 0), (1, 1), (0, 2), (2, 2), (3, 1), (1, 0)])
+  >>> regions = voronoi_diagram(points)
+  >>> pprint([region.wkt for region in regions])
+  ['POLYGON ((2 1, 2 0.5, 0.5 0.5, 0 1, 1 2, 2 1))',
+   'POLYGON ((6 5, 6 -3, 3.75 -3, 2 0.5, 2 1, 6 5))',
+   'POLYGON ((0.5 -3, -3 -3, -3 1, 0 1, 0.5 0.5, 0.5 -3))',
+   'POLYGON ((3.75 -3, 0.5 -3, 0.5 0.5, 2 0.5, 3.75 -3))',
+   'POLYGON ((-3 1, -3 5, 1 5, 1 2, 0 1, -3 1))',
+   'POLYGON ((1 5, 6 5, 2 1, 1 2, 1 5))']
+
+
 Nearest points
 --------------
 
diff --git a/shapely/ctypes_declarations.py b/shapely/ctypes_declarations.py
index 7ecfb3a..eb6c803 100644
--- a/shapely/ctypes_declarations.py
+++ b/shapely/ctypes_declarations.py
@@ -248,6 +248,10 @@ def prototype(lgeos, geos_version):
         lgeos.GEOSDelaunayTriangulation.restype = c_void_p
         lgeos.GEOSDelaunayTriangulation.argtypes = [c_void_p, c_double, c_int]
 
+    if geos_version >= (3, 5, 0):
+        lgeos.GEOSVoronoiDiagram.restype = c_void_p
+        lgeos.GEOSVoronoiDiagram.argtypes = [c_void_p, c_void_p, c_double, c_int]
+
     lgeos.GEOSLineMerge.restype = c_void_p
     lgeos.GEOSLineMerge.argtypes = [c_void_p]
 
diff --git a/shapely/geos.py b/shapely/geos.py
index d16041f..17ebb22 100644
--- a/shapely/geos.py
+++ b/shapely/geos.py
@@ -847,6 +847,7 @@ class LGEOS350(LGEOS340):
     def __init__(self, dll):
         super(LGEOS350, self).__init__(dll)
         self.methods['clip_by_rect'] = self.GEOSClipByRect
+        self.methods['voronoi_diagram'] = self.GEOSVoronoiDiagram
 
 
 if geos_version >= (3, 5, 0):
diff --git a/shapely/ops.py b/shapely/ops.py
index 4601f5a..104089e 100644
--- a/shapely/ops.py
+++ b/shapely/ops.py
@@ -14,7 +14,7 @@ from shapely.algorithms.polylabel import polylabel
 
 __all__ = ['cascaded_union', 'linemerge', 'operator', 'polygonize',
            'polygonize_full', 'transform', 'unary_union', 'triangulate',
-           'split', 'section']
+           'voronoi_diagram', 'split', 'section']
 
 
 class CollectionOperator(object):
@@ -172,6 +172,64 @@ def triangulate(geom, tolerance=0.0, edges=False):
     gc = geom_factory(func(geom._geom, tolerance, int(edges)))
     return [g for g in gc.geoms]
 
+
+def voronoi_diagram(geom, envelope=None, tolerance=0.0, edges=False):
+    """
+    Constructs a Voronoi Diagram [1] from the given geometry.
+    Returns a list of geometries.
+
+    Parameters
+    ----------
+    geom: geometry
+        the input geometry whose vertices will be used to calculate
+        the final diagram.
+    envelope: geometry, None
+        clipping envelope for the returned diagram, automatically
+        determined if None. The diagram will be clipped to the larger
+        of this envelope or an envelope surrounding the sites.
+    tolerance: float, 0.0
+        sets the snapping tolerance used to improve the robustness
+        of the computation. A tolerance of 0.0 specifies that no
+        snapping will take place.
+    edges: bool, False
+        If False, return regions as polygons. Else, return only
+        edges e.g. LineStrings.
+
+    GEOS documentation can be found at [2]
+
+    Returns
+    -------
+    GeometryCollection
+        geometries representing the Voronoi regions.
+
+    Notes
+    -----
+    The tolerance `argument` can be finicky and is known to cause the
+    algorithm to fail in several cases. If you're using `tolerance`
+    and getting a failure, try removing it. The test cases in
+    tests/test_voronoi_diagram.py show more details.
+
+
+    References
+    ----------
+    [1] https://en.wikipedia.org/wiki/Voronoi_diagram
+    [2] https://geos.osgeo.org/doxygen/geos__c_8h_source.html  (line 730)
+    """
+    func = lgeos.methods['voronoi_diagram']
+    envelope = envelope._geom if envelope else None
+    try:
+        result = geom_factory(func(geom._geom, envelope, tolerance, int(edges)))
+    except ValueError:
+        errstr = "Could not create Voronoi Diagram with the specified inputs."
+        if tolerance:
+            errstr += " Try running again with default tolerance value."
+        raise ValueError(errstr)
+
+    if result.type != 'GeometryCollection':
+        return GeometryCollection([result])
+    return result
+
+
 class ValidateOp(object):
     def __call__(self, this):
         return lgeos.GEOSisValidReason(this._geom)


Test patch:
diff --git a/tests/test_voronoi_diagram.py b/tests/test_voronoi_diagram.py
new file mode 100644
index 0000000..d3fd902
--- /dev/null
+++ b/tests/test_voronoi_diagram.py
@@ -0,0 +1,135 @@
+"""
+Test cases for Voronoi Diagram creation.
+
+Overall, I'm trying less to test the correctness of the result
+and more to cover input cases and behavior, making sure
+that we return a sane result without error or raise a useful one.
+"""
+
+import pytest
+
+from shapely.geos import geos_version
+from shapely.geometry import MultiPoint
+from shapely.wkt import loads as load_wkt
+
+from shapely.ops import voronoi_diagram
+
+requires_geos_35 = pytest.mark.skipif(geos_version < (3, 5, 0), reason='GEOS >= 3.5.0 is required.')
+
+
+@requires_geos_35
+def test_no_regions():
+    mp = MultiPoint(points=[(0.5, 0.5)])
+    regions = voronoi_diagram(mp)
+
+    assert len(regions) == 0
+
+
+@requires_geos_35
+def test_two_regions():
+    mp = MultiPoint(points=[(0.5, 0.5), (1.0, 1.0)])
+    regions = voronoi_diagram(mp)
+
+    assert len(regions) == 2
+
+
+@requires_geos_35
+def test_edges():
+    mp = MultiPoint(points=[(0.5, 0.5), (1.0, 1.0)])
+    regions = voronoi_diagram(mp, edges=True)
+
+    assert len(regions) == 1
+    assert all(r.type == 'LineString' for r in regions)
+
+
+@requires_geos_35
+def test_smaller_envelope():
+    mp = MultiPoint(points=[(0.5, 0.5), (1.0, 1.0)])
+    poly = load_wkt('POLYGON ((0 0, 0.5 0, 0.5 0.5, 0 0.5, 0 0))')
+
+    regions = voronoi_diagram(mp, envelope=poly)
+
+    assert len(regions) == 2
+    assert sum(r.area for r in regions) > poly.area
+
+
+@requires_geos_35
+def test_larger_envelope():
+    """When the envelope we specify is larger than the
+    area of the input feature, the created regions should
+    expand to fill that area."""
+    mp = MultiPoint(points=[(0.5, 0.5), (1.0, 1.0)])
+    poly = load_wkt('POLYGON ((0 0, 2 0, 2 2, 0 2, 0 0))')
+
+    regions = voronoi_diagram(mp, envelope=poly)
+
+    assert len(regions) == 2
+    assert sum(r.area for r in regions) == poly.area
+
+
+@requires_geos_35
+def test_from_polygon():
+    poly = load_wkt('POLYGON ((0 0, 1 0, 1 1, 0 1, 0 0))')
+    regions = voronoi_diagram(poly)
+
+    assert len(regions) == 4
+
+
+@requires_geos_35
+def test_from_polygon_with_enough_tolerance():
+    poly = load_wkt('POLYGON ((0 0, 0.5 0, 0.5 0.5, 0 0.5, 0 0))')
+    regions = voronoi_diagram(poly, tolerance=1.0)
+
+    assert len(regions) == 2
+
+
+@requires_geos_35
+def test_from_polygon_without_enough_tolerance():
+    poly = load_wkt('POLYGON ((0 0, 0.5 0, 0.5 0.5, 0 0.5, 0 0))')
+    with pytest.raises(ValueError) as exc:
+        voronoi_diagram(poly, tolerance=0.6)
+
+    assert exc.match("Could not create Voronoi Diagram with the specified inputs. "
+                     "Try running again with default tolerance value.")
+
+
+@requires_geos_35
+def test_from_polygon_without_floating_point_coordinates():
+    poly = load_wkt('POLYGON ((0 0, 1 0, 1 1, 0 1, 0 0))')
+    with pytest.raises(ValueError) as exc:
+        voronoi_diagram(poly, tolerance=0.1)
+
+    assert exc.match("Could not create Voronoi Diagram with the specified inputs. "
+                     "Try running again with default tolerance value")
+
+
+@requires_geos_35
+def test_from_multipoint_without_floating_point_coordinates():
+    """A Multipoint with the same "shape" as the above Polygon raises the same error..."""
+    mp = load_wkt('MULTIPOINT (0 0, 1 0, 1 1, 0 1)')
+
+    with pytest.raises(ValueError) as exc:
+        voronoi_diagram(mp, tolerance=0.1)
+
+    assert exc.match("Could not create Voronoi Diagram with the specified inputs. "
+                     "Try running again with default tolerance value")
+
+
+@requires_geos_35
+def test_from_multipoint_with_tolerace_without_floating_point_coordinates():
+    """This multipoint will not work with a tolerance value."""
+    mp = load_wkt('MULTIPOINT (0 0, 1 0, 1 2, 0 1)')
+    with pytest.raises(ValueError) as exc:
+        voronoi_diagram(mp, tolerance=0.1)
+
+    assert exc.match("Could not create Voronoi Diagram with the specified inputs. "
+                     "Try running again with default tolerance value")
+
+
+@requires_geos_35
+def test_from_multipoint_without_tolerace_without_floating_point_coordinates():
+    """But it's fine without it."""
+    mp = load_wkt('MULTIPOINT (0 0, 1 0, 1 2, 0 1)')
+    regions = voronoi_diagram(mp)
+    assert len(regions) == 4
+