Repository: geopandas/geopandas (python)

Issue:
ENH: Topological simplification
`GeoSeries.simplify` does not have a notion of topology in relation to other elements, so simplification often introduces gaps between polygons and similar mess:

```py
df = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
df = df[df.continent == 'Africa']
df.simplify(10).plot()
```
![image](https://user-images.githubusercontent.com/36797143/79693907-38888400-8265-11ea-89f9-35f811c1ba21.png)

However, that is likely not what user wants in that case. I just came across new `topojson` package (https://github.com/mattijn/topojson) by @mattijn which support GeoDataFrames directly and can do topological simplification:

```
import topojson as tp
topo = tp.Topology(df, prequantize=False)
topo.toposimplify(5).to_gdf().plot()
```
![image](https://user-images.githubusercontent.com/36797143/79693967-8604f100-8265-11ea-9386-28b4b6d6a6ad.png)

I believe there's a lot of people who would like to see this. It would be nice to expose this functionality as `GeoSeries.toposimplify()`. Should be straightforward to implement.

Apart from that we could use the package to dump Topojson (#645, #610) to file.

Pull request description:
ENH: add GeoSeries.simplify_coverage for topological simplification
Closes #1387 

Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index a797e06d..0cbdee19 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -21,6 +21,8 @@ New features and improvements:
 - Added `maximum_inscribed_circle` method from shapely to GeoSeries/GeoDataframe (#3544).
 - Added `minimum_clearance_line` method from shapely to GeoSeries/GeoDataframe (#3543).
 - Added ``method`` and ``keep_collapsed`` argument to ``make_valid`` (#3548).
+- Added `simplify_coverage` method for topological simplification of polygonal coverages
+  to GeoSeries/GeoDataframe (#3541).
 
 Bug fixes:
 
diff --git a/doc/source/docs/reference/geoseries.rst b/doc/source/docs/reference/geoseries.rst
index 48ebea22..06ac4e89 100644
--- a/doc/source/docs/reference/geoseries.rst
+++ b/doc/source/docs/reference/geoseries.rst
@@ -123,6 +123,7 @@ Constructive methods and attributes
    GeoSeries.segmentize
    GeoSeries.shortest_line
    GeoSeries.simplify
+   GeoSeries.simplify_coverage
    GeoSeries.snap
    GeoSeries.transform
 
diff --git a/geopandas/array.py b/geopandas/array.py
index 44011c88..698c5ecf 100644
--- a/geopandas/array.py
+++ b/geopandas/array.py
@@ -871,6 +871,18 @@ class GeometryArray(ExtensionArray):
             crs=self.crs,
         )
 
+    def simplify_coverage(self, tolerance, simplify_boundary=True):
+        if not (SHAPELY_GE_21 and GEOS_GE_312):
+            raise ImportError(
+                "'simplify_coverage' requires shapely>=2.1 and GEOS>=3.12."
+            )
+        return GeometryArray(
+            shapely.coverage_simplify(
+                self._data, tolerance, simplify_boundary=simplify_boundary
+            ),
+            crs=self.crs,
+        )
+
     def project(self, other, normalized=False):
         if isinstance(other, GeometryArray):
             other = other._data
diff --git a/geopandas/base.py b/geopandas/base.py
index 554e0027..b71227a8 100644
--- a/geopandas/base.py
+++ b/geopandas/base.py
@@ -5289,9 +5289,14 @@ GeometryCollection
         to the straight line is smaller than `tolerance`. It does not
         move any points and it always preserves endpoints of
         the original line or polygon.
-        See http://shapely.readthedocs.io/en/latest/manual.html#object.simplify
+        See https://shapely.readthedocs.io/en/latest/manual.html#object.simplify
         for details
 
+        Simplifies individual geometries independently, without considering
+        the topology of a potential polygonal coverage. If you would like to treat
+        the ``GeoSeries`` as a coverage and simplify its edges, while preserving the
+        coverage topology, see :meth:`simplify_coverage`.
+
         Parameters
         ----------
         tolerance : float
@@ -5311,6 +5316,10 @@ GeometryCollection
         coordinates: two geometries differing only in order of coordinates may be
         simplified differently.
 
+        See also
+        --------
+        simplify_coverage : simplify geometries using coverage simplification
+
         Examples
         --------
         >>> from shapely.geometry import Point, LineString
@@ -5331,6 +5340,80 @@ GeometryCollection
             "simplify", self, tolerance=tolerance, preserve_topology=preserve_topology
         )
 
+    def simplify_coverage(self, tolerance, simplify_boundary=True):
+        """Returns a ``GeoSeries`` containing a simplified representation of
+        polygonal coverage.
+
+        Assumes that the ``GeoSeries`` forms a polygonal coverage. Under this
+        assumption, the method simplifies the edges using the Visvalingam-Whyatt
+        algorithm, while preserving a valid coverage. In the most simplified case,
+        polygons are reduced to triangles.
+
+        A ``GeoSeries`` of valid polygons is considered a coverage if the polygons are:
+
+        * **Non-overlapping** - polygons do not overlap (their interiors do not
+          intersect)
+        * **Edge-Matched** - vertices along shared edges are identical
+
+        The method allows simplification of all edges including the outer boundaries of
+        the coverage or simplification of only the inner (shared) edges.
+
+        If there are other geometry types than Polygons or MultiPolygons present, the
+        method will raise an error.
+
+        If the geometry is polygonal but does not form a valid coverage due to overlaps,
+        it will be simplified but it may result in invalid coverage topology.
+
+        .. versionadded:: 1.1.0
+
+        Parameters
+        ----------
+        tolerance : float
+            The degree of simplification roughly equal to the square root of the area
+            of triangles that will be removed. It has the same units
+            as the coordinate reference system of the GeoSeries.
+            For example, using `tolerance=100` in a projected CRS with meters
+            as units means a distance of 100 meters in reality.
+        simplify_boundary: bool (default True)
+            By default (True), simplifies both internal edges of the coverage as well
+            as its boundary. If set to False, only simplifies internal edges.
+
+
+        See also
+        --------
+        simplify : simplification of individual geometries
+
+        Examples
+        --------
+        >>> from shapely.geometry import Polygon
+        >>> s = geopandas.GeoSeries(
+        ...     [
+        ...         Polygon([(0, 0), (1, 1.1), (2, 0), (1.5, 1), (2, 2), (0, 2)]),
+        ...         Polygon([(2, 0), (4, 0), (4, 2), (2, 2), (1.5, 1)]),
+        ...     ]
+        ... )
+        >>> s
+        0    POLYGON ((0 0, 1 1.1, 2 0, 1.5 1, 2 2, 0 2, 0 0))
+        1           POLYGON ((2 0, 4 0, 4 2, 2 2, 1.5 1, 2 0))
+        dtype: geometry
+
+        >>> s.simplify_coverage(1)
+        0         POLYGON ((2 0, 2 2, 0 2, 2 0))
+        1    POLYGON ((2 0, 4 0, 4 2, 2 2, 2 0))
+        dtype: geometry
+
+        >>> s.simplify_coverage(1, simplify_boundary=False)
+        0    POLYGON ((2 0, 2 2, 0 2, 0 0, 1 1.1, 2 0))
+        1           POLYGON ((2 0, 4 0, 4 2, 2 2, 2 0))
+        dtype: geometry
+        """
+        return _delegate_geo_method(
+            "simplify_coverage",
+            self,
+            tolerance=tolerance,
+            simplify_boundary=simplify_boundary,
+        )
+
     def relate(self, other, align=None):
         """
         Returns the DE-9IM intersection matrices for the geometries


Test patch:
diff --git a/geopandas/tests/test_geom_methods.py b/geopandas/tests/test_geom_methods.py
index 772ad6ab..0a5c6df4 100644
--- a/geopandas/tests/test_geom_methods.py
+++ b/geopandas/tests/test_geom_methods.py
@@ -1483,6 +1483,38 @@ class TestGeomMethods:
         ):
             self.g1.simplify(Series([0.1], index=[99]))
 
+    @pytest.mark.skipif(
+        not (GEOS_GE_312 and SHAPELY_GE_21), reason="GEOS 3.12 and shapely 2.1 needed."
+    )
+    def test_simplify_coverage(self):
+        s = GeoSeries(
+            [
+                shapely.Polygon(
+                    [(0, 0), (10, 1), (20, 0), (20, 10), (10, 5), (0, 10), (0, 0)]
+                ),
+                shapely.Polygon(
+                    [(0, 10), (10, 5), (20, 10), (20, 20), (0, 20), (0, 10)]
+                ),
+            ]
+        )
+        e = GeoSeries(
+            [
+                shapely.Polygon([(0, 0), (20, 0), (20, 10), (0, 10)]),
+                shapely.Polygon([(0, 10), (20, 10), (20, 20), (0, 20)]),
+            ]
+        )
+        assert_geoseries_equal(s.simplify_coverage(8), e.normalize())
+
+        e_boundary = GeoSeries(
+            [
+                shapely.Polygon([(0, 0), (10, 1), (20, 0), (20, 10), (0, 10)]),
+                shapely.Polygon([(0, 10), (20, 10), (20, 20), (0, 20)]),
+            ]
+        )
+        assert_geoseries_equal(
+            s.simplify_coverage(8, simplify_boundary=False), e_boundary.normalize()
+        )
+
     def test_envelope(self):
         e = self.g3.envelope
         assert np.all(e.geom_equals(self.sq))