From 855bdbb431f882def0ab71d81637b3c5cdcb274d Mon Sep 17 00:00:00 2001 From: bstadlbauer <11799671+bstadlbauer@users.noreply.github.com> Date: Tue, 1 Mar 2022 14:55:13 +0100 Subject: [PATCH] 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 * 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 --- doc/source/docs/reference/geoseries.rst | 1 + geopandas/_vectorized.py | 15 ++++++ geopandas/array.py | 5 ++ geopandas/base.py | 69 +++++++++++++++++++++++++ geopandas/tests/test_geom_methods.py | 13 +++++ 5 files changed, 103 insertions(+) diff --git a/doc/source/docs/reference/geoseries.rst b/doc/source/docs/reference/geoseries.rst index 0173a8c..1965d11 100644 --- a/doc/source/docs/reference/geoseries.rst +++ b/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 diff --git a/geopandas/_vectorized.py b/geopandas/_vectorized.py index 457a81f..159e578 100644 --- a/geopandas/_vectorized.py +++ b/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) diff --git a/geopandas/array.py b/geopandas/array.py index 5f3f24b..9b47d6a 100644 --- a/geopandas/array.py +++ b/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 diff --git a/geopandas/base.py b/geopandas/base.py index 2edf647..dab414b 100644 --- a/geopandas/base.py +++ b/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 # diff --git a/geopandas/tests/test_geom_methods.py b/geopandas/tests/test_geom_methods.py index 05475a9..05e3d5c 100644 --- a/geopandas/tests/test_geom_methods.py +++ b/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)