Repository: Toblerity/Shapely (python)

Issue:
Add GEOSMakeValid to Shapely
It would be nice to have this functionality available directly in Shapely. It requires [GEOS 3.8](https://github.com/libgeos/geos/blob/312c085bdc9606896281c7d0610c1a1b3931e014/NEWS#L21).

Should it live in `ops`, `validation`, or somewhere else?

Pull request:
Add implementation for make_valid
Closes #883.

Patch:
diff --git a/CHANGES.txt b/CHANGES.txt
index ae11305..39112a3 100644
--- a/CHANGES.txt
+++ b/CHANGES.txt
@@ -7,6 +7,7 @@ Changes
 New features:
 
 - shapely.ops.voronoi_diagram() computes the Voronoi Diagram of a geometry or geometry collection (#833, #851).
+- shapely.validation.make_valid() fixes invalid geometries (#883)
 
 1.7.1 (TBD)
 -----------
diff --git a/docs/code/make_valid_geometrycollection.py b/docs/code/make_valid_geometrycollection.py
new file mode 100644
index 0000000..17c2f1b
--- /dev/null
+++ b/docs/code/make_valid_geometrycollection.py
@@ -0,0 +1,33 @@
+from shapely.geometry import Polygon
+from shapely.validation import make_valid
+
+from matplotlib import pyplot
+from descartes.patch import PolygonPatch
+from figures import SIZE, BLUE, RED, set_limits, plot_line
+
+
+invalid_poly = Polygon([(0, 2), (0, 1), (2, 0), (0, 0), (0, 2)])
+valid_poly = make_valid(invalid_poly)
+
+fig = pyplot.figure(1, figsize=SIZE, dpi=90)
+fig.set_frameon(True)
+
+
+invalid_ax = fig.add_subplot(121)
+
+patch = PolygonPatch(invalid_poly, facecolor=BLUE, edgecolor=BLUE, alpha=0.5, zorder=2)
+invalid_ax.add_patch(patch)
+
+set_limits(invalid_ax, -1, 3, -1, 3)
+
+
+valid_ax = fig.add_subplot(122)
+
+patch = PolygonPatch(valid_poly[0], facecolor=BLUE, edgecolor=BLUE, alpha=0.5, zorder=2)
+valid_ax.add_patch(patch)
+
+plot_line(valid_ax, valid_poly[1], color=RED, linewidth=1)
+
+set_limits(valid_ax, -1, 3, -1, 3)
+
+pyplot.show()
\ No newline at end of file
diff --git a/docs/code/make_valid_multipolygon.py b/docs/code/make_valid_multipolygon.py
new file mode 100644
index 0000000..b156e46
--- /dev/null
+++ b/docs/code/make_valid_multipolygon.py
@@ -0,0 +1,33 @@
+from shapely.geometry import Polygon
+from shapely.validation import make_valid
+
+from matplotlib import pyplot
+from descartes.patch import PolygonPatch
+from figures import SIZE, BLUE, RED, set_limits
+
+invalid_poly = Polygon([(0, 0), (0, 2), (1, 1), (2, 2), (2, 0), (1, 1), (0, 0)])
+valid_poly = make_valid(invalid_poly)
+
+fig = pyplot.figure(1, figsize=SIZE, dpi=90)
+fig.set_frameon(True)
+
+
+invalid_ax = fig.add_subplot(121)
+
+patch = PolygonPatch(invalid_poly, facecolor=BLUE, edgecolor=BLUE, alpha=0.5, zorder=2)
+invalid_ax.add_patch(patch)
+
+set_limits(invalid_ax, -1, 3, -1, 3)
+
+
+valid_ax = fig.add_subplot(122)
+
+patch = PolygonPatch(valid_poly[0], facecolor=BLUE, edgecolor=BLUE, alpha=0.5, zorder=2)
+valid_ax.add_patch(patch)
+
+patch = PolygonPatch(valid_poly[1], facecolor=RED, edgecolor=RED, alpha=0.5, zorder=2)
+valid_ax.add_patch(patch)
+
+set_limits(valid_ax, -1, 3, -1, 3)
+
+pyplot.show()
\ No newline at end of file
diff --git a/docs/manual.rst b/docs/manual.rst
index f13c66f..8baf5cc 100644
--- a/docs/manual.rst
+++ b/docs/manual.rst
@@ -2460,6 +2460,47 @@ be parsed out.
   >>> explain_validity(p)
   'Ring Self-intersection[1 1]'
 
+.. function:: validation.make_valid(ob)
+
+  Returns a valid representation of the geometry, if it is invalid.
+  If it is valid, the input geometry will be returned.
+
+  In many cases, in order to create a valid geometry, the input geometry
+  must be split into multiple parts or multiple geometries. If the geometry
+  must be split into multiple parts of the same geometry type, then a multi-part
+  geometry (e.g. a MultiPolygon) will be returned. if the geometry must be split
+  into multiple parts of different types, then a GeometryCollection will be returned.
+
+  For example, this operation on a geometry with a bow-tie structure:
+
+.. code-block:: pycon
+
+  >>> from shapely.validation import make_valid
+  >>> coords = [(0, 0), (0, 2), (1, 1), (2, 2), (2, 0), (1, 1), (0, 0)]
+  >>> p = Polygon(coords)
+  >>> str(make_valid(p))
+  'MULTIPOLYGON (((0 0, 0 2, 1 1, 0 0)), ((1 1, 2 2, 2 0, 1 1)))'
+
+  Yields a MultiPolygon with two parts:
+
+.. plot:: code/make_valid_multipolygon.py
+
+  While this operation:
+
+.. code-block:: pycon
+
+  >>> from shapely.validation import make_valid
+  >>> coords = [(0, 2), (0, 1), (2, 0), (0, 0), (0, 2)]
+  >>> p = Polygon(coords)
+  >>> str(make_valid(p))
+
+  Yields a GeometryCollection with a Polygon and a LineString:
+
+.. plot:: code/make_valid_geometrycollection.py
+
+  `New in version 1.8`
+  `Requires GEOS > 3.8`
+
 The Shapely version, GEOS library version, and GEOS C API version are
 accessible via :attr:`shapely.__version__`,
 :attr:`shapely.geos.geos_version_string`, and
diff --git a/requirements-dev.txt b/requirements-dev.txt
index 0b686d5..4ee040a 100644
--- a/requirements-dev.txt
+++ b/requirements-dev.txt
@@ -1,5 +1,6 @@
 cython==0.29.14
 descartes==1.0.1
+matplotlib
 numpy>=1.4.1
 pytest
 pytest-cov
diff --git a/shapely/ctypes_declarations.py b/shapely/ctypes_declarations.py
index 9606412..f9d199a 100644
--- a/shapely/ctypes_declarations.py
+++ b/shapely/ctypes_declarations.py
@@ -545,3 +545,8 @@ def prototype(lgeos, geos_version):
 
         lgeos.GEOSMinimumClearance.argtypes = [c_void_p]
         lgeos.GEOSMinimumClearance.restype = c_double
+
+    if geos_version > (3, 8, 0):
+        lgeos.GEOSMakeValid.restype = c_void_p
+        lgeos.GEOSMakeValid.argtypes = [c_void_p]
+
diff --git a/shapely/geos.py b/shapely/geos.py
index 41e1b33..791fda8 100644
--- a/shapely/geos.py
+++ b/shapely/geos.py
@@ -871,7 +871,20 @@ class LGEOS360(LGEOS350):
         self.methods['minimum_clearance'] = self.GEOSMinimumClearance
 
 
-if geos_version >= (3, 6, 0):
+class LGEOS380(LGEOS360):
+    """Proxy for GEOS 3.8.0-CAPI-1.13.0"""
+
+    geos_version = (3, 8, 0)
+    geos_capi_version = (1, 13, 0)
+
+    def __init__(self, dll):
+        super(LGEOS380, self).__init__(dll)
+        self.methods['make_valid'] = self.GEOSMakeValid
+
+
+if geos_version >= (3, 8, 0):
+    L = LGEOS380
+elif geos_version >= (3, 6, 0):
     L = LGEOS360
 elif geos_version >= (3, 5, 0):
     L = LGEOS350
diff --git a/shapely/validation.py b/shapely/validation.py
index 541a779..6ef1f7a 100644
--- a/shapely/validation.py
+++ b/shapely/validation.py
@@ -1,6 +1,55 @@
 # TODO: allow for implementations using other than GEOS
-
+from shapely.geometry.base import geom_factory
 from shapely.geos import lgeos
 
+__all__ = ['explain_validity', 'make_valid']
+
+
 def explain_validity(ob):
+    """
+    Explain the validity of the input geometry, if it is invalid.
+    This will describe why the geometry is invalid, and might
+    include a location if there is a self-intersection or a
+    ring self-intersection.
+
+    Parameters
+    ----------
+    ob: Geometry
+        A shapely geometry object
+
+    Returns
+    -------
+    str
+        A string describing the reason the geometry is invalid.
+
+    """
     return lgeos.GEOSisValidReason(ob._geom)
+
+
+def make_valid(ob):
+    """
+    Make the input geometry valid according to the GEOS MakeValid algorithm.
+
+    If the input geometry is already valid, then it will be returned.
+
+    If the geometry must be split into multiple parts of the same type to be made valid,
+    then a multi-part geometry will be returned.
+
+    If the geometry must be split into multiple parts of different types to be made valid,
+    then a GeometryCollection will be returned.
+
+    Parameters
+    ----------
+    ob : Geometry
+        A shapely geometry object which should be made valid. If the object is already valid,
+        it will be returned as-is.
+
+    Returns
+    -------
+    Geometry
+        The input geometry, made valid according to the GEOS MakeValid algorithm.
+
+    """
+    if ob.is_valid:
+        return ob
+    return geom_factory(lgeos.GEOSMakeValid(ob._geom))


Test patch:
diff --git a/tests/conftest.py b/tests/conftest.py
index 5a6094d..6319788 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -5,6 +5,7 @@ import pytest
 from shapely.geos import geos_version
 
 
+requires_geos_38 = pytest.mark.skipif(geos_version < (3, 8, 0), reason='GEOS >= 3.8.0 is required.')
 requires_geos_342 = pytest.mark.skipif(geos_version < (3, 4, 2), reason="GEOS > 3.4.2 is required.")
 
 
diff --git a/tests/test_make_valid.py b/tests/test_make_valid.py
new file mode 100644
index 0000000..3d1108c
--- /dev/null
+++ b/tests/test_make_valid.py
@@ -0,0 +1,20 @@
+from shapely.geometry import Polygon
+from shapely.validation import make_valid
+
+from tests.conftest import requires_geos_38
+
+
+@requires_geos_38
+def test_make_valid_invalid_input():
+    geom = Polygon([(0, 0), (0, 2), (1, 1), (2, 2), (2, 0), (1, 1), (0, 0)])
+    valid = make_valid(geom)
+    assert len(valid) == 2
+    assert all(geom.type == 'Polygon' for geom in valid)
+
+
+@requires_geos_38
+def test_make_valid_input():
+    geom = Polygon([(0, 0), (1, 0), (1, 1), (0, 1), (0, 0)])
+    valid = make_valid(geom)
+    assert id(valid) == id(geom)
+