From 4a650af20223ffd376c7f0f4307d31a243cbb98e Mon Sep 17 00:00:00 2001
From: James McBride <jdmcbr@gmail.com>
Date: Thu, 18 Feb 2021 12:36:56 -0800
Subject: [PATCH] ENH: Add warning if keep_geom_type drops geometries during
 overlay (#1554)

* Add warning if keep_geom_type drops geometries during overlay

* Expand keep_geom_type warning message

* Add test for keep_geom_type dropping geoms

* Set keep_geom_type=None by default

With True, no warning raised. With None, will drop geoms of other types,
but warn the user.

* Update geopandas/tests/test_overlay.py

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

Co-authored-by: Joris Van den Bossche <jorisvandenbossche@gmail.com>
Co-authored-by: Martin Fleischmann <martin@martinfleischmann.net>
---
 geopandas/tests/test_overlay.py | 18 ++++++++++++++++-
 geopandas/tools/overlay.py      | 34 ++++++++++++++++++++++++++-------
 2 files changed, 44 insertions(+), 8 deletions(-)

diff --git a/geopandas/tests/test_overlay.py b/geopandas/tests/test_overlay.py
index 6ba0c39..b828112 100644
--- a/geopandas/tests/test_overlay.py
+++ b/geopandas/tests/test_overlay.py
@@ -382,10 +382,26 @@ def test_correct_index(dfs):
     expected = GeoDataFrame(
         [[1, 1, i1], [3, 2, i2]], columns=["col3", "col2", "geometry"]
     )
-    result = overlay(df3, df2)
+    result = overlay(df3, df2, keep_geom_type=True)
     assert_geodataframe_equal(result, expected)
 
 
+def test_warn_on_keep_geom_type(dfs):
+
+    df1, df2 = dfs
+    polys3 = GeoSeries(
+        [
+            Polygon([(1, 1), (3, 1), (3, 3), (1, 3)]),
+            Polygon([(-1, 1), (1, 1), (1, 3), (-1, 3)]),
+            Polygon([(3, 3), (5, 3), (5, 5), (3, 5)]),
+        ]
+    )
+    df3 = GeoDataFrame({"geometry": polys3})
+
+    with pytest.warns(UserWarning, match="`keep_geom_type=True` in overlay"):
+        overlay(df2, df3, keep_geom_type=None)
+
+
 @pytest.mark.parametrize(
     "geom_types", ["polys", "poly_line", "poly_point", "line_poly", "point_poly"]
 )
diff --git a/geopandas/tools/overlay.py b/geopandas/tools/overlay.py
index d3c0bf8..9aebbb3 100644
--- a/geopandas/tools/overlay.py
+++ b/geopandas/tools/overlay.py
@@ -137,7 +137,7 @@ def _overlay_union(df1, df2):
     return dfunion.reindex(columns=columns)
 
 
-def overlay(df1, df2, how="intersection", keep_geom_type=True):
+def overlay(df1, df2, how="intersection", keep_geom_type=None):
     """Perform spatial overlay between two GeoDataFrames.
 
     Currently only supports data GeoDataFrames with uniform geometry types,
@@ -156,7 +156,9 @@ def overlay(df1, df2, how="intersection", keep_geom_type=True):
         'identity', 'symmetric_difference' or 'difference'.
     keep_geom_type : bool
         If True, return only geometries of the same geometry type as df1 has,
-        if False, return all resulting gemetries.
+        if False, return all resulting geometries. Default is None,
+        which will set keep_geom_type to True but warn upon dropping
+        geometries.
 
     Returns
     -------
@@ -241,6 +243,12 @@ def overlay(df1, df2, how="intersection", keep_geom_type=True):
     if not _check_crs(df1, df2):
         _crs_mismatch_warn(df1, df2, stacklevel=3)
 
+    if keep_geom_type is None:
+        keep_geom_type = True
+        keep_geom_type_warning = True
+    else:
+        keep_geom_type_warning = False
+
     polys = ["Polygon", "MultiPolygon"]
     lines = ["LineString", "MultiLineString", "LinearRing"]
     points = ["Point", "MultiPoint"]
@@ -280,20 +288,32 @@ def overlay(df1, df2, how="intersection", keep_geom_type=True):
         exploded = result.reset_index(drop=True).explode()
         exploded = exploded.reset_index(level=0)
 
-        type = df1.geom_type.iloc[0]
-        if type in polys:
+        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)]
-        elif type in lines:
+        elif geom_type in lines:
             exploded = exploded.loc[exploded.geom_type.isin(lines)]
-        elif type in points:
+        elif geom_type in points:
             exploded = exploded.loc[exploded.geom_type.isin(points)]
         else:
-            raise TypeError("`keep_geom_type` does not support {}.".format(type))
+            raise TypeError("`keep_geom_type` does not support {}.".format(geom_type))
 
         # 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]
+            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),
+                UserWarning,
+                stacklevel=2,
+            )
+
     result.reset_index(drop=True, inplace=True)
     result.drop(["__idx1", "__idx2"], axis=1, inplace=True)
     return result