Repository: geopandas/geopandas

Issue:
REGR: performance degrading after updating from version 0.8.2 to 0.9, using the "overlay"-method
- [x] I have searched the [geopandas] tag on [StackOverflow](https://stackoverflow.com/questions/tagged/geopandas) and [GIS StackExchange](https://gis.stackexchange.com/questions/tagged/geopandas) for similar questions.

- [x] I have asked my usage related question on [StackOverflow](https://stackoverflow.com) or [GIS StackExhange](https://gis.stackexchange.com).

---

#### Question about geopandas

Hello! The OS I'm working on is Microsoft Windows 10 LTSC. I often need to split GeoDataFrames with a minimum number of 350.000 Points into approximately 100 Polygons with:
```
cell_points = gpd.overlay(gdf_points, gdf_polygons, how='intersection')
```
So I was very happy when I found this bug report on the geopandas GitHub site: https://github.com/geopandas/geopandas/issues/1793. A bit further down in the thread, I found that geopandas 0.8.2 was released which works with pygeos. After I had updated from 0.8.1 to 0.8.2 my code executed much much faster. But when I recently updated to version 0.9, I quickly noticed that my code executed much much slower. Luckily, switching back to version 0.8.2 brought back the execution improvement.
What I did was to set up a virtual environment. Then I installed first of all geopandas from conda-forge, then I installed pygeos from conda-forge. This is the code that I used to measure the time:
```
import geopandas as gpd
import numpy as np
import pandas as pd
from shapely.geometry import Polygon
from time import perf_counter


print(gpd.__version__)


def create_Polygon(row):
   return Polygon([(i[0], i[1]) for i in row.reshape((4, 2))])


# #### Generate random Points inside (0, 10) ######
df_points = pd.DataFrame(data=np.random.rand(100000, 2) * 10, columns=['X', 'Y'])
gdf_points = gpd.GeoDataFrame(df_points, geometry=gpd.points_from_xy(df_points['X'], df_points['Y']))

# #### Generate 100 Polygons (Squares with edge length of 1) inside (0, 10) #####
step = 1
arr = np.zeros(shape=(100, 8))
counter = 0
for i in range(10):
   for j in range(10):
       arr[counter] = [i, j, i+step, j, i+step, j+step, i, j+step]
       counter += 1
gdf_polygons = gpd.GeoDataFrame(data=list(range(100)),
                               geometry=np.apply_along_axis(create_Polygon, axis=1, arr=arr))

# #### Intersect Polygons with Points #####
time_start = perf_counter()
cell_points = gpd.overlay(gdf_points, gdf_polygons, how='intersection')
print(f'Time: {(perf_counter() - time_start):.2f} seconds')
```
The execution time of version 0.9 was approximately 40 seconds. The execution time of version 0.8.2 was approximately 0.6 seconds! (Both with Python version 3.9.6, executed in PyCharm). My question is: Can someone recreate this? Is there maybe something else that I should have installed, when updating to geopandas version 0.9 to keep the performance improvements from version 0.8.2?


Discussion:
jorisvandenbossche: We should also add more tests for this. Currently we only have `test_keep_geom_type_geometry_collection` that uses the code path that I am changing here, and that's a simple case with 1 geometry in the result (eg we should at least have a test where the result gives multiple rows with a mixture of collections and normal geometries).
jdmcbr: @jorisvandenbossche Thanks for taking this on! I'd been planning to take a shot at this issue myself last weekend, but you got here well before I did (and then I wound up with the flu or something like the flu [though not COVID] the week before last so wouldn't have made it here anyway). 

> 3) add a benchmark for this.

Are you running the `asv` benchmarks regularly yourself? I've used them on occasion, but the fact that we don't have them running on every PR makes it easy to miss regressions. So far as I can see, there isn't a standard way to add `asv` to the suite of CI checks (though there's some fresh-ish activity since I had last looked, e.g. https://github.com/airspeed-velocity/asv/issues/955). I confirmed our current overlay benchmarks wouldn't have caught this particular issue, and I'd be happy to add a new benchmark that does run into the performance regression I introduced, but I'm just generally curious about how others are using `asv`. 

> We should also add more tests for this.

I can also update the test fixtures (let me know if you want me to add that here, or put in a separate PR). 
jorisvandenbossche: > > add a benchmark for this.
>
> Are you running the asv benchmarks regularly yourself?

To be honest, the second I posted that comment, I directly thought: 1) hmm, I should have checked if the current asv benchmark didn't already would have catched this, and 2) we should actually run those benchmarks ... :-)

So no, I am not regularly running those benchmarks. We should make it a habit to at least run them when prepping for a new release (like right now https://github.com/geopandas/geopandas/issues/2076). Ideally we would also regularly run them in an automated way, but that's a bit a harder problem.

And indeed, the existing benchmarks would actually have catched this regression:

```
$ asv continuous -f 1.1 v0.8.2 v0.9.0 -b overlay
...
       before           after         ratio
     [928fe6c2]       [ec4c6805]
     <v0.8.2^0>       <v0.9.0^0>
+        238±30ms         792±60ms     3.32  overlay.Countries.time_overlay('intersection')
+      10.9±0.5ms         25.5±1ms     2.35  overlay.Small.time_overlay('intersection')
```

> I can also update the test fixtures (let me know if you want me to add that here, or put in a separate PR).

That's certainly welcome if you would have time to look at this!

jdmcbr: >To be honest, the second I posted that comment, I directly thought: 1) hmm, I should have checked if the current asv benchmark didn't already would have catched this, and 2) we should actually run those benchmarks ... :-)
> And indeed, the existing benchmarks would actually have catched this regression:

Yeah, I realized that in the course of putting in #2118. It had been long enough since I'd run the benchmarks (but I _thought_ I remembered how to run them) that my first attempt didn't correctly run the benchmark against your branch, leading me to the mistaken conclusion that the existing benchmarks didn't catch this regression. 

> That's certainly welcome if you would have time to look at this!

Do you prefer committing here, or in a separate PR? 
jorisvandenbossche: Separate PR is fine
jorisvandenbossche: @martinfleis if you have some time to take a look here? (this one is important for 0.10 I think)
martinfleis: Can you check and update the docstring example? I think it is just the matter of order of coordinates.
jorisvandenbossche: So I added already one extra test (where the result has both normal polygons as geometrycollection), but we can certainly still further expand the coverage (eg with other geometry types as polygons)

Patch:
diff --git a/geopandas/tools/overlay.py b/geopandas/tools/overlay.py
index ef47fa1e..ca516f81 100644
--- a/geopandas/tools/overlay.py
+++ b/geopandas/tools/overlay.py
@@ -185,8 +185,8 @@ def overlay(df1, df2, how="intersection", keep_geom_type=None, make_valid=True):
     1       2.0       1.0  POLYGON ((2.00000 2.00000, 2.00000 3.00000, 3....
     2       2.0       2.0  POLYGON ((4.00000 4.00000, 4.00000 3.00000, 3....
     3       1.0       NaN  POLYGON ((2.00000 0.00000, 0.00000 0.00000, 0....
-    4       2.0       NaN  MULTIPOLYGON (((4.00000 3.00000, 4.00000 2.000...
-    5       NaN       1.0  MULTIPOLYGON (((3.00000 2.00000, 3.00000 1.000...
+    4       2.0       NaN  MULTIPOLYGON (((3.00000 3.00000, 4.00000 3.000...
+    5       NaN       1.0  MULTIPOLYGON (((2.00000 2.00000, 3.00000 2.000...
     6       NaN       2.0  POLYGON ((3.00000 5.00000, 5.00000 5.00000, 5....
 
     >>> geopandas.overlay(df1, df2, how='intersection')
@@ -198,8 +198,8 @@ def overlay(df1, df2, how="intersection", keep_geom_type=None, make_valid=True):
     >>> geopandas.overlay(df1, df2, how='symmetric_difference')
        df1_data  df2_data                                           geometry
     0       1.0       NaN  POLYGON ((2.00000 0.00000, 0.00000 0.00000, 0....
-    1       2.0       NaN  MULTIPOLYGON (((4.00000 3.00000, 4.00000 2.000...
-    2       NaN       1.0  MULTIPOLYGON (((3.00000 2.00000, 3.00000 1.000...
+    1       2.0       NaN  MULTIPOLYGON (((3.00000 3.00000, 4.00000 3.000...
+    2       NaN       1.0  MULTIPOLYGON (((2.00000 2.00000, 3.00000 2.000...
     3       NaN       2.0  POLYGON ((3.00000 5.00000, 5.00000 5.00000, 5....
 
     >>> geopandas.overlay(df1, df2, how='difference')
@@ -213,7 +213,7 @@ def overlay(df1, df2, how="intersection", keep_geom_type=None, make_valid=True):
     1       2.0       1.0  POLYGON ((2.00000 2.00000, 2.00000 3.00000, 3....
     2       2.0       2.0  POLYGON ((4.00000 4.00000, 4.00000 3.00000, 3....
     3       1.0       NaN  POLYGON ((2.00000 0.00000, 0.00000 0.00000, 0....
-    4       2.0       NaN  MULTIPOLYGON (((4.00000 3.00000, 4.00000 2.000...
+    4       2.0       NaN  MULTIPOLYGON (((3.00000 3.00000, 4.00000 3.000...
 
     See also
     --------
@@ -317,32 +317,64 @@ def overlay(df1, df2, how="intersection", keep_geom_type=None, make_valid=True):
             result = dfunion[dfunion["__idx1"].notnull()].copy()
 
     if keep_geom_type:
-        key_order = result.keys()
-        exploded = result.reset_index(drop=True).explode()
-        exploded = exploded.reset_index(level=0)
+        geom_type = df1.geom_type.iloc[0]
+
+        # First we filter the geometry types inside GeometryCollections objects
+        # (e.g. GeometryCollection([polygon, point]) -> polygon)
+        # we do this separately on only the relevant rows, as this is an expensive
+        # operation (an expensive no-op for geometry types other than collections)
+        is_collection = result.geom_type == "GeometryCollection"
+        if is_collection.any():
+            geom_col = result._geometry_column_name
+            collections = result[[geom_col]][is_collection]
+
+            exploded = collections.reset_index(drop=True).explode()
+            exploded = exploded.reset_index(level=0)
+
+            orig_num_geoms_exploded = exploded.shape[0]
+            if geom_type in polys:
+                exploded = exploded.loc[exploded.geom_type.isin(polys)]
+            elif geom_type in lines:
+                exploded = exploded.loc[exploded.geom_type.isin(lines)]
+            elif geom_type in points:
+                exploded = exploded.loc[exploded.geom_type.isin(points)]
+            else:
+                raise TypeError(
+                    "`keep_geom_type` does not support {}.".format(geom_type)
+                )
+            num_dropped_collection = orig_num_geoms_exploded - exploded.shape[0]
+
+            # level_0 created with above reset_index operation
+            # and represents the original geometry collections
+            # TODO avoiding dissolve to call unary_union in this case could further
+            # improve performance (we only need to collect geometries in their
+            # respective Multi version)
+            dissolved = exploded.dissolve(by="level_0")
+            result.loc[is_collection, geom_col] = dissolved[geom_col].values
+        else:
+            num_dropped_collection = 0
 
+        # Now we filter all geometries (in theory we don't need to do this
+        # again for the rows handled above for GeometryCollections, but filtering
+        # them out is probably more expensive as simply including them when this
+        # is typically about only a few rows)
         orig_num_geoms = result.shape[0]
-        geom_type = df1.geom_type.iloc[0]
         if geom_type in polys:
-            exploded = exploded.loc[exploded.geom_type.isin(polys)]
+            result = result.loc[result.geom_type.isin(polys)]
         elif geom_type in lines:
-            exploded = exploded.loc[exploded.geom_type.isin(lines)]
+            result = result.loc[result.geom_type.isin(lines)]
         elif geom_type in points:
-            exploded = exploded.loc[exploded.geom_type.isin(points)]
+            result = result.loc[result.geom_type.isin(points)]
         else:
             raise TypeError("`keep_geom_type` does not support {}.".format(geom_type))
+        num_dropped = orig_num_geoms - result.shape[0]
 
-        # level_0 created with above reset_index operation
-        # and represents the original geometry collections
-        result = exploded.dissolve(by="level_0")[key_order]
-
-        if (result.shape[0] != orig_num_geoms) and keep_geom_type_warning:
-            num_dropped = orig_num_geoms - result.shape[0]
+        if (num_dropped > 0 or num_dropped_collection > 0) and keep_geom_type_warning:
             warnings.warn(
                 "`keep_geom_type=True` in overlay resulted in {} dropped "
                 "geometries of different geometry types than df1 has. "
                 "Set `keep_geom_type=False` to retain all "
-                "geometries".format(num_dropped),
+                "geometries".format(num_dropped + num_dropped_collection),
                 UserWarning,
                 stacklevel=2,
             )


Test patch:
diff --git a/geopandas/tests/test_overlay.py b/geopandas/tests/test_overlay.py
index a07bf62b..ec9b6820 100644
--- a/geopandas/tests/test_overlay.py
+++ b/geopandas/tests/test_overlay.py
@@ -575,6 +575,11 @@ def test_keep_geom_type_geometry_collection():
     df1 = read_file(os.path.join(DATA, "geom_type", "df1.geojson"))
     df2 = read_file(os.path.join(DATA, "geom_type", "df2.geojson"))
 
+    with pytest.warns(UserWarning, match="`keep_geom_type=True` in overlay"):
+        intersection = overlay(df1, df2, keep_geom_type=None)
+    assert len(intersection) == 1
+    assert (intersection.geom_type == "Polygon").all()
+
     intersection = overlay(df1, df2, keep_geom_type=True)
     assert len(intersection) == 1
     assert (intersection.geom_type == "Polygon").all()
@@ -584,6 +589,44 @@ def test_keep_geom_type_geometry_collection():
     assert (intersection.geom_type == "GeometryCollection").all()
 
 
+def test_keep_geom_type_geometry_collection2():
+    polys1 = [
+        box(0, 0, 1, 1),
+        box(1, 1, 3, 3).union(box(1, 3, 5, 5)),
+    ]
+
+    polys2 = [
+        box(0, 0, 1, 1),
+        box(3, 1, 4, 2).union(box(4, 1, 5, 4)),
+    ]
+    df1 = GeoDataFrame({"left": [0, 1], "geometry": polys1})
+    df2 = GeoDataFrame({"right": [0, 1], "geometry": polys2})
+
+    result1 = overlay(df1, df2, keep_geom_type=True)
+    expected1 = GeoDataFrame(
+        {
+            "left": [0, 1],
+            "right": [0, 1],
+            "geometry": [box(0, 0, 1, 1), box(4, 3, 5, 4)],
+        }
+    )
+    assert_geodataframe_equal(result1, expected1)
+
+    result1 = overlay(df1, df2, keep_geom_type=False)
+    expected1 = GeoDataFrame(
+        {
+            "left": [0, 1, 1],
+            "right": [0, 0, 1],
+            "geometry": [
+                box(0, 0, 1, 1),
+                Point(1, 1),
+                GeometryCollection([box(4, 3, 5, 4), LineString([(3, 1), (3, 2)])]),
+            ],
+        }
+    )
+    assert_geodataframe_equal(result1, expected1)
+
+
 @pytest.mark.parametrize("make_valid", [True, False])
 def test_overlap_make_valid(make_valid):
     bowtie = Polygon([(1, 1), (9, 9), (9, 1), (1, 9), (1, 1)])