ENH: Add fast rectangle clipping using .clip_by_rect() (#1928)

* TST: Add test for rectangle clipping

* ENH: Add `GeometryArray().clip_by_rect()`

* ENH: Add `GeoPandasBase().clip_by_rect()`

Closes #1902

* Change signatures to use tuple rectangle

* Change signatures to use tuple rectangle

* Raise exception when shapely < 1.7

Exception is raised as `shapely.ops.clip_by_rect` was introduced in shapely 1.7

* Reorder tests to make them easier to group

* Parameterize all test which use single_rectangle_df to also test with bounds

* Paramterize test with larger single rectangle

* Add bounding box clip capability to clip using .clip_by_rect()

* Add documentation for `.clip_by_rect()` and changes in `clip()`

* Skip clip tests if shapely < 1.7

* Improve test coverage

* Remove shapely >1.7 compatibility checks

Can be removed as #2358 bumped shapely to ">= 1.7"

* Apply suggestions from code review

Co-authored-by: Martin Fleischmann <martin@martinfleischmann.net>

* Only check for empty geometries after `.clip_by_rect()`

* Move changes to `clip.py` to its own PR

* Change docstring of `base.clip_by_rect()` according to PR comments

* add interpshinx links

Co-authored-by: Martin Fleischmann <martin@martinfleischmann.net>

doc/source/docs/reference/geoseries.rst

@@ -69,6 +69,7 @@ Set-theoretic Methods
 .. autosummary::
    :toctree: api/
 
+   GeoSeries.clip_by_rect
    GeoSeries.difference
    GeoSeries.intersection
    GeoSeries.symmetric_difference

geopandas/_vectorized.py

@@ -11,6 +11,7 @@ import pandas as pd
 
 import shapely.geometry
 import shapely.geos
+import shapely.ops
 import shapely.wkb
 import shapely.wkt
 
@@ -728,6 +729,20 @@ def almost_equals(self, other, decimal):
 #
 
 
+def clip_by_rect(data, xmin, ymin, xmax, ymax):
+    if compat.USE_PYGEOS:
+        return pygeos.clip_by_rect(data, xmin, ymin, xmax, ymax)
+    else:
+        clipped_geometries = np.empty(len(data), dtype=object)
+        clipped_geometries[:] = [
+            shapely.ops.clip_by_rect(s, xmin, ymin, xmax, ymax)
+            if s is not None
+            else None
+            for s in data
+        ]
+        return clipped_geometries
+
+
 def difference(data, other):
     if compat.USE_PYGEOS:
         return _binary_method("difference", data, other)

geopandas/array.py

@@ -564,6 +564,11 @@ class GeometryArray(ExtensionArray):
     # Binary operations that return new geometries
     #
 
+    def clip_by_rect(self, xmin, ymin, xmax, ymax):
+        return GeometryArray(
+            vectorized.clip_by_rect(self.data, xmin, ymin, xmax, ymax), crs=self.crs
+        )
+
     def difference(self, other):
         return GeometryArray(
             self._binary_method("difference", self, other), crs=self.crs

geopandas/base.py

@@ -2536,6 +2536,75 @@ GeometryCollection
         """
         return _binary_geo("intersection", self, other, align)
 
+    def clip_by_rect(self, xmin, ymin, xmax, ymax):
+        """Returns a ``GeoSeries`` of the portions of geometry within the given
+        rectangle.
+
+        Note that the results are not exactly equal to
+        :meth:`~GeoSeries.intersection()`. E.g. in edge cases,
+        :meth:`~GeoSeries.clip_by_rect()` will not return a point just touching the
+        rectangle. Check the examples section below for some of these exceptions.
+
+        The geometry is clipped in a fast but possibly dirty way. The output is not
+        guaranteed to be valid. No exceptions will be raised for topological errors.
+
+        Note: empty geometries or geometries that do not overlap with the specified
+        bounds will result in ``GEOMETRYCOLLECTION EMPTY``.
+
+        Parameters
+        ----------
+        xmin: float
+            Minimum x value of the rectangle
+        ymin: float
+            Minimum y value of the rectangle
+        xmax: float
+            Maximum x value of the rectangle
+        ymax: float
+            Maximum y value of the rectangle
+
+        Returns
+        -------
+        GeoSeries
+
+        Examples
+        --------
+        >>> from shapely.geometry import Polygon, LineString, Point
+        >>> s = geopandas.GeoSeries(
+        ...     [
+        ...         Polygon([(0, 0), (2, 2), (0, 2)]),
+        ...         Polygon([(0, 0), (2, 2), (0, 2)]),
+        ...         LineString([(0, 0), (2, 2)]),
+        ...         LineString([(2, 0), (0, 2)]),
+        ...         Point(0, 1),
+        ...     ],
+        ...     crs=3857,
+        ... )
+        >>> bounds = (0, 0, 1, 1)
+        >>> s
+        0    POLYGON ((0.000 0.000, 2.000 2.000, 0.000 2.00...
+        1    POLYGON ((0.000 0.000, 2.000 2.000, 0.000 2.00...
+        2                LINESTRING (0.000 0.000, 2.000 2.000)
+        3                LINESTRING (2.000 0.000, 0.000 2.000)
+        4                                  POINT (0.000 1.000)
+        dtype: geometry
+        >>> s.clip_by_rect(*bounds)
+        0    POLYGON ((0.000 0.000, 0.000 1.000, 1.000 1.00...
+        1    POLYGON ((0.000 0.000, 0.000 1.000, 1.000 1.00...
+        2                LINESTRING (0.000 0.000, 1.000 1.000)
+        3                             GEOMETRYCOLLECTION EMPTY
+        4                             GEOMETRYCOLLECTION EMPTY
+        dtype: geometry
+
+        See also
+        --------
+        GeoSeries.intersection
+        """
+        from .geoseries import GeoSeries
+
+        geometry_array = GeometryArray(self.geometry.values)
+        clipped_geometry = geometry_array.clip_by_rect(xmin, ymin, xmax, ymax)
+        return GeoSeries(clipped_geometry.data, index=self.index, crs=self.crs)
+
     #
     # Other operations
     #

geopandas/tests/test_geom_methods.py

@@ -85,8 +85,12 @@ class TestGeomMethods:
         self.g8 = GeoSeries([self.t1, self.t5])
         self.empty = GeoSeries([])
         self.all_none = GeoSeries([None, None])
+        self.all_geometry_collection_empty = GeoSeries(
+            [GeometryCollection([]), GeometryCollection([])]
+        )
         self.empty_poly = Polygon()
         self.g9 = GeoSeries(self.g0, index=range(1, 8))
+        self.g10 = GeoSeries([self.t1, self.t4])
 
         # Crossed lines
         self.l3 = LineString([(0, 0), (1, 1)])
@@ -254,6 +258,15 @@ class TestGeomMethods:
             assert len(self.g0.intersection(self.g9, align=True) == 8)
         assert len(self.g0.intersection(self.g9, align=False) == 7)
 
+    def test_clip_by_rect(self):
+        self._test_binary_topological(
+            "clip_by_rect", self.g1, self.g10, *self.sq.bounds
+        )
+        # self.g1 and self.t3.bounds do not intersect
+        self._test_binary_topological(
+            "clip_by_rect", self.all_geometry_collection_empty, self.g1, *self.t3.bounds
+        )
+
     def test_union_series(self):
         self._test_binary_topological("union", self.sq, self.g1, self.g2)