Repository: geopandas/geopandas
Pull request #2148 (2021-09-30 19:17:12)

BUG: Fix multipoint clipping / PERF: query for all geometry types at once
This is based on https://github.com/geopandas/geopandas/pull/1909 (thanks to @bstadlbauer for that!), but only keeping the part about doing the actual clipping for all geometry types together, instead of splitting the original geodataframe in chunks per geometry type, clip those separately, and then concat (+ sort to preserve original order) them again.

This should keep the performance improvement of the original PR (and also the MultiPoint bug), but the diff should be a bit easier to follow since it has less other changes.

While doing this, I noticed that I could do the clipping with preserving the original order by setting the intersections with a mask, which avoids concatting / sorting the point/non-point parts. So that's something added compared to https://github.com/geopandas/geopandas/pull/1909


Patch:
diff --git a/geopandas/geodataframe.py b/geopandas/geodataframe.py
index f4c51790..6b57e037 100644
--- a/geopandas/geodataframe.py
+++ b/geopandas/geodataframe.py
@@ -47,6 +47,20 @@ def _ensure_geometry(data, crs=None):
             return out
 
 
+def _crs_mismatch_warning():
+    # TODO: raise error in 0.9 or 0.10.
+    warnings.warn(
+        "CRS mismatch between CRS of the passed geometries "
+        "and 'crs'. Use 'GeoDataFrame.set_crs(crs, "
+        "allow_override=True)' to overwrite CRS or "
+        "'GeoDataFrame.to_crs(crs)' to reproject geometries. "
+        "CRS mismatch will raise an error in the future versions "
+        "of GeoPandas.",
+        FutureWarning,
+        stacklevel=3,
+    )
+
+
 class GeoDataFrame(GeoPandasBase, DataFrame):
     """
     A GeoDataFrame object is a pandas.DataFrame that has a column
@@ -103,9 +117,9 @@ class GeoDataFrame(GeoPandasBase, DataFrame):
 
     _geometry_column_name = DEFAULT_GEO_COLUMN_NAME
 
-    def __init__(self, *args, geometry=None, crs=None, **kwargs):
+    def __init__(self, data=None, *args, geometry=None, crs=None, **kwargs):
         with compat.ignore_shapely2_warnings():
-            super().__init__(*args, **kwargs)
+            super().__init__(data, *args, **kwargs)
 
         # need to set this before calling self['geometry'], because
         # getitem accesses crs
@@ -117,6 +131,15 @@ class GeoDataFrame(GeoPandasBase, DataFrame):
         # but within a try/except because currently non-geometries are
         # allowed in that case
         # TODO do we want to raise / return normal DataFrame in this case?
+
+        # if gdf passed in and geo_col is set, we use that for geometry
+        if geometry is None and isinstance(data, GeoDataFrame):
+            self._geometry_column_name = data._geometry_column_name
+            if crs is not None and data.crs != crs:
+                _crs_mismatch_warning()
+                # TODO: raise error in 0.9 or 0.10.
+            return
+
         if geometry is None and "geometry" in self.columns:
             # Check for multiple columns with name "geometry". If there are,
             # self["geometry"] is a gdf and constructor gets recursively recalled
@@ -136,16 +159,7 @@ class GeoDataFrame(GeoPandasBase, DataFrame):
                     and crs
                     and not self["geometry"].values.crs == crs
                 ):
-                    warnings.warn(
-                        "CRS mismatch between CRS of the passed geometries "
-                        "and 'crs'. Use 'GeoDataFrame.set_crs(crs, "
-                        "allow_override=True)' to overwrite CRS or "
-                        "'GeoDataFrame.to_crs(crs)' to reproject geometries. "
-                        "CRS mismatch will raise an error in the future versions "
-                        "of GeoPandas.",
-                        FutureWarning,
-                        stacklevel=2,
-                    )
+                    _crs_mismatch_warning()
                     # TODO: raise error in 0.9 or 0.10.
                 self["geometry"] = _ensure_geometry(self["geometry"].values, crs)
             except TypeError:
@@ -165,16 +179,7 @@ class GeoDataFrame(GeoPandasBase, DataFrame):
                 and crs
                 and not geometry.crs == crs
             ):
-                warnings.warn(
-                    "CRS mismatch between CRS of the passed geometries "
-                    "and 'crs'. Use 'GeoDataFrame.set_crs(crs, "
-                    "allow_override=True)' to overwrite CRS or "
-                    "'GeoDataFrame.to_crs(crs)' to reproject geometries. "
-                    "CRS mismatch will raise an error in the future versions "
-                    "of GeoPandas.",
-                    FutureWarning,
-                    stacklevel=2,
-                )
+                _crs_mismatch_warning()
                 # TODO: raise error in 0.9 or 0.10.
             self.set_geometry(geometry, inplace=True)
 
diff --git a/geopandas/tools/clip.py b/geopandas/tools/clip.py
index 12d41c63..c8750239 100644
--- a/geopandas/tools/clip.py
+++ b/geopandas/tools/clip.py
@@ -7,50 +7,22 @@ A module to clip vector data using GeoPandas.
 """
 import warnings
 
-import numpy as np
-import pandas as pd
-
 from shapely.geometry import Polygon, MultiPolygon
 
 from geopandas import GeoDataFrame, GeoSeries
 from geopandas.array import _check_crs, _crs_mismatch_warn
 
 
-def _clip_points(gdf, poly):
-    """Clip point geometry to the polygon extent.
-
-    Clip an input point GeoDataFrame to the polygon extent of the poly
-    parameter. Points that intersect the poly geometry are extracted with
-    associated attributes and returned.
-
-    Parameters
-    ----------
-    gdf : GeoDataFrame, GeoSeries
-        Composed of point geometry that will be clipped to the poly.
-
-    poly : (Multi)Polygon
-        Reference geometry used to spatially clip the data.
-
-    Returns
-    -------
-    GeoDataFrame
-        The returned GeoDataFrame is a subset of gdf that intersects
-        with poly.
-    """
-    return gdf.iloc[gdf.sindex.query(poly, predicate="intersects")]
-
-
-def _clip_line_poly(gdf, poly):
-    """Clip line and polygon geometry to the polygon extent.
+def _clip_gdf_with_polygon(gdf, poly):
+    """Clip geometry to the polygon extent.
 
-    Clip an input line or polygon to the polygon extent of the poly
-    parameter. Parts of Lines or Polygons that intersect the poly geometry are
-    extracted with associated attributes and returned.
+    Clip an input GeoDataFrame to the polygon extent of the poly
+    parameter.
 
     Parameters
     ----------
     gdf : GeoDataFrame, GeoSeries
-        Line or polygon geometry that is clipped to poly.
+        Dataframe to clip.
 
     poly : (Multi)Polygon
         Reference polygon for clipping.
@@ -63,13 +35,23 @@ def _clip_line_poly(gdf, poly):
     """
     gdf_sub = gdf.iloc[gdf.sindex.query(poly, predicate="intersects")]
 
+    # For performance reasons points don't need to be intersected with poly
+    non_point_mask = gdf_sub.geom_type != "Point"
+
+    if not non_point_mask.any():
+        # only points, directly return
+        return gdf_sub
+
     # Clip the data with the polygon
     if isinstance(gdf_sub, GeoDataFrame):
         clipped = gdf_sub.copy()
-        clipped[gdf.geometry.name] = gdf_sub.intersection(poly)
+        clipped.loc[
+            non_point_mask, clipped._geometry_column_name
+        ] = gdf_sub.geometry.values[non_point_mask].intersection(poly)
     else:
         # GeoSeries
-        clipped = gdf_sub.intersection(poly)
+        clipped = gdf_sub.copy()
+        clipped[non_point_mask] = gdf_sub.values[non_point_mask].intersection(poly)
 
     return clipped
 
@@ -149,42 +131,11 @@ def clip(gdf, mask, keep_geom_type=False):
     else:
         poly = mask
 
-    geom_types = gdf.geometry.type
-    poly_idx = np.asarray((geom_types == "Polygon") | (geom_types == "MultiPolygon"))
-    line_idx = np.asarray(
-        (geom_types == "LineString")
-        | (geom_types == "LinearRing")
-        | (geom_types == "MultiLineString")
-    )
-    point_idx = np.asarray((geom_types == "Point") | (geom_types == "MultiPoint"))
-    geomcoll_idx = np.asarray((geom_types == "GeometryCollection"))
-
-    if point_idx.any():
-        point_gdf = _clip_points(gdf[point_idx], poly)
-    else:
-        point_gdf = None
-
-    if poly_idx.any():
-        poly_gdf = _clip_line_poly(gdf[poly_idx], poly)
-    else:
-        poly_gdf = None
-
-    if line_idx.any():
-        line_gdf = _clip_line_poly(gdf[line_idx], poly)
-    else:
-        line_gdf = None
-
-    if geomcoll_idx.any():
-        geomcoll_gdf = _clip_line_poly(gdf[geomcoll_idx], poly)
-    else:
-        geomcoll_gdf = None
-
-    order = pd.Series(range(len(gdf)), index=gdf.index)
-    concat = pd.concat([point_gdf, line_gdf, poly_gdf, geomcoll_gdf])
+    clipped = _clip_gdf_with_polygon(gdf, poly)
 
     if keep_geom_type:
-        geomcoll_concat = (concat.geom_type == "GeometryCollection").any()
-        geomcoll_orig = geomcoll_idx.any()
+        geomcoll_concat = (clipped.geom_type == "GeometryCollection").any()
+        geomcoll_orig = (gdf.geom_type == "GeometryCollection").any()
 
         new_collection = geomcoll_concat and not geomcoll_orig
 
@@ -210,9 +161,9 @@ def clip(gdf, mask, keep_geom_type=False):
             # Check how many geometry types are in the clipped GeoDataFrame
             clip_types_total = sum(
                 [
-                    concat.geom_type.isin(polys).any(),
-                    concat.geom_type.isin(lines).any(),
-                    concat.geom_type.isin(points).any(),
+                    clipped.geom_type.isin(polys).any(),
+                    clipped.geom_type.isin(lines).any(),
+                    clipped.geom_type.isin(points).any(),
                 ]
             )
 
@@ -226,21 +177,10 @@ def clip(gdf, mask, keep_geom_type=False):
             elif new_collection or more_types:
                 orig_type = gdf.geom_type.iloc[0]
                 if new_collection:
-                    concat = concat.explode()
+                    clipped = clipped.explode()
                 if orig_type in polys:
-                    concat = concat.loc[concat.geom_type.isin(polys)]
+                    clipped = clipped.loc[clipped.geom_type.isin(polys)]
                 elif orig_type in lines:
-                    concat = concat.loc[concat.geom_type.isin(lines)]
+                    clipped = clipped.loc[clipped.geom_type.isin(lines)]
 
-    # Return empty GeoDataFrame or GeoSeries if no shapes remain
-    if len(concat) == 0:
-        return gdf.iloc[:0]
-
-    # Preserve the original order of the input
-    if isinstance(concat, GeoDataFrame):
-        concat["_order"] = order
-        return concat.sort_values(by="_order").drop(columns="_order")
-    else:
-        concat = GeoDataFrame(geometry=concat)
-        concat["_order"] = order
-        return concat.sort_values(by="_order").geometry
+    return clipped


Test patch:
diff --git a/geopandas/tests/test_crs.py b/geopandas/tests/test_crs.py
index c8aa321b..ca237b0c 100644
--- a/geopandas/tests/test_crs.py
+++ b/geopandas/tests/test_crs.py
@@ -220,16 +220,16 @@ class TestGeometryArrayCRS:
         assert df.geometry.values.crs == self.osgb
 
         # different passed CRS than array CRS is ignored
-        with pytest.warns(FutureWarning):
+        with pytest.warns(FutureWarning, match="CRS mismatch"):
             df = GeoDataFrame(geometry=s, crs=4326)
         assert df.crs == self.osgb
         assert df.geometry.crs == self.osgb
         assert df.geometry.values.crs == self.osgb
-        with pytest.warns(FutureWarning):
+        with pytest.warns(FutureWarning, match="CRS mismatch"):
             GeoDataFrame(geometry=s, crs=4326)
-        with pytest.warns(FutureWarning):
+        with pytest.warns(FutureWarning, match="CRS mismatch"):
             GeoDataFrame({"data": [1, 2], "geometry": s}, crs=4326)
-        with pytest.warns(FutureWarning):
+        with pytest.warns(FutureWarning, match="CRS mismatch"):
             GeoDataFrame(df, crs=4326).crs
 
         # manually change CRS
diff --git a/geopandas/tests/test_geodataframe.py b/geopandas/tests/test_geodataframe.py
index cbfebbd5..ea98f478 100644
--- a/geopandas/tests/test_geodataframe.py
+++ b/geopandas/tests/test_geodataframe.py
@@ -930,11 +930,9 @@ class TestConstructor:
             res = GeoDataFrame(df, geometry="other_geom")
             check_geodataframe(res, "other_geom")
 
-        # when passing GeoDataFrame with custom geometry name to constructor
-        # an invalid geodataframe is the result TODO is this desired ?
+        # gdf from gdf should preserve active geometry column name
         df = GeoDataFrame(gpdf)
-        with pytest.raises(AttributeError):
-            df.geometry
+        check_geodataframe(df, "other_geom")
 
     def test_only_geometry(self):
         exp = GeoDataFrame(
diff --git a/geopandas/tools/tests/test_clip.py b/geopandas/tools/tests/test_clip.py
index 255ba5ca..f71ed394 100644
--- a/geopandas/tools/tests/test_clip.py
+++ b/geopandas/tools/tests/test_clip.py
@@ -7,7 +7,14 @@ import numpy as np
 import pandas as pd
 
 import shapely
-from shapely.geometry import Polygon, Point, LineString, LinearRing, GeometryCollection
+from shapely.geometry import (
+    Polygon,
+    Point,
+    LineString,
+    LinearRing,
+    GeometryCollection,
+    MultiPoint,
+)
 
 import geopandas
 from geopandas import GeoDataFrame, GeoSeries, clip
@@ -287,6 +294,15 @@ def test_clip_multipoint(single_rectangle_gdf, multi_point):
     assert clipped.geom_type[0] == "MultiPoint"
     assert hasattr(clipped, "attr")
     # All points should intersect the clip geom
+    assert len(clipped) == 2
+    clipped_mutltipoint = MultiPoint(
+        [
+            Point(2, 2),
+            Point(3, 4),
+            Point(9, 8),
+        ]
+    )
+    assert clipped.iloc[0].geometry.wkt == clipped_mutltipoint.wkt
     assert all(clipped.intersects(single_rectangle_gdf.unary_union))