Repository: geopandas/geopandas

Issue:
ENH: Speedup to overlay/identity
While reading the code of `overlay`, I noticed that for "identity", "union" is calculated (= "intersection" + "symmetric_difference") and then part of the result is removed. This could be replaced by "intersection" + "difference", which should be faster?

Discussion:
@BMarques: Do you mean that doing something like this 

intersections = gpd.overlay(input_gdf, identity_gdf, how='intersection')
difference = gpd.overlay(input_gdf, identity_gdf, how='difference')
result_gdf = intersections.append(difference, ignore_index=True)

will yield the same results, but quicker?@martinfleis: I think that it was directed towards internal code that lives in https://github.com/geopandas/geopandas/blob/e82949886a9937b39d1d5dce087d1ad5ccaf01d4/geopandas/tools/overlay.py#L326-L328

where `_overlay_union` consists of

https://github.com/geopandas/geopandas/blob/e82949886a9937b39d1d5dce087d1ad5ccaf01d4/geopandas/tools/overlay.py#L140-L141

@theroggy I suppose it is worth trying.@BMarques: Sorry, I was aware of what he was talking about, I kinda super imposed because I'm stuck currently with the slow performance of overlay identity (takes more than 3 hours in my use case), and I was wondering if doing what I said manually also did the same thing, based on his comment. Will stop "hijacking" the issue now.@theroggy: > Sorry, I was aware of what he was talking about, I kinda super imposed because I'm stuck currently with the slow performance of overlay identity (takes more than 3 hours in my use case), and I was wondering if doing what I said manually also did the same thing, based on his comment. Will stop "hijacking" the issue now.

@BMarques Yes, doing this "manually" should give the same performance improvement... if it does improve performance. So, if you would try it on your use case... please report back here because if it is faster it would be a confirmation that it is a good idea to change the internal implementation. If performance would be the same, the only improvement would (probably) be decreased memory usage. If it would slower: bad idea to change, or possibly a case to check in geos why `difference` is slower than `symmetric_difference`, as that sounds quite counter-intuitif.

Regarding your slow "identity" use case: if the above wouldn't help (sufficiently): feel free to start a thread about it in the [discussions](https://github.com/geopandas/geopandas/discussions) section.@Bimbam360: Chiming in as @BMarques never replied.
In my example I've just replaced:

```
overlay_start = time.perf_counter()
gdf_layers['union_cadastres'] = gpd.overlay(
	gdf_layers['cadastre'],
	gdf_layers['union'],
	how="identity",
	keep_geom_type=False
).explode(index_parts=False)
gdf_layers['union_cadastres'] = gdf_layers['union_cadastres'][gdf_layers['union_cadastres'].geometry.geom_type == "Polygon"]
print(f"Overlay Identity took {time.perf_counter() - overlay_start:.2f} seconds")
```

with
```
overlay_start = time.perf_counter()
intersections = gpd.overlay(gdf_layers['cadastre'], gdf_layers['union'], how='intersection', keep_geom_type=False)
difference = gpd.overlay(gdf_layers['cadastre'], gdf_layers['union'], how='difference', keep_geom_type=False)
gdf_layers['union_cadastres_test'] = pd.concat([intersections, difference], ignore_index=True).explode(index_parts=False)
gdf_layers['union_cadastres_test'] = gdf_layers['union_cadastres_test'][gdf_layers['union_cadastres_test'].geometry.geom_type == "Polygon"]
print(f"Faster Overlay Identity took {time.perf_counter() - overlay_start:.2f} seconds")
```

Overlay Identity took 1228.24 seconds
Faster Overlay Identity took 119.04 seconds
A 10x speedup! and to confirm:

`gdf_layers['union_cadastres_test'].equals(gdf_layers['union_cadastres']) `
True

```
from geopandas.testing import assert_geodataframe_equal
assert_geodataframe_equal(gdf_layers['union_cadastres'], gdf_layers['union_cadastres_test'], 
                                      True, 'equiv','equiv',True,False,True,True,True,True)
```
No Error.

Shapely = 2.0.1
Geopandas = 0.14.4
@BMarques: Sorry for not coming back in the meantime, I managed to find out another solution that solved my particular issue, but good to know about that.@theroggy: I implemented the speedup in https://github.com/geopandas/geopandas/pull/3504

The 10x increase will be a relatively specific case, probably the other way around for this case will only be a 10% increase, but nonetheless seems like a good improvement...

Solution patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 4224727e..fe41c2be 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -11,6 +11,7 @@ New features and improvements:
 - Add ``grid_size`` parameter to ``union_all`` and ``dissolve`` (#3445).
 - `GeoDataFrame.plot` now supports `pd.Index` as an input for the `column` keyword (#3463).
 - Avoid change of the plot aspect when plotting missing values (#3438).
+- Improve performance of `overlay` with `how=identity` (#3504).
 
 Bug fixes:
 
diff --git a/geopandas/geodataframe.py b/geopandas/geodataframe.py
index f716ea8c..dfe68ba7 100644
--- a/geopandas/geodataframe.py
+++ b/geopandas/geodataframe.py
@@ -2831,11 +2831,11 @@ chicago_w_groceries[chicago_w_groceries["community"] == "UPTOWN"]
 
         >>> df1.overlay(df2, how='identity')
            df1_data  df2_data                                           geometry
-        0       1.0       1.0                POLYGON ((2 2, 2 1, 1 1, 1 2, 2 2))
-        1       2.0       1.0                POLYGON ((2 2, 2 3, 3 3, 3 2, 2 2))
-        2       2.0       2.0                POLYGON ((4 4, 4 3, 3 3, 3 4, 4 4))
-        3       1.0       NaN      POLYGON ((2 0, 0 0, 0 2, 1 2, 1 1, 2 1, 2 0))
-        4       2.0       NaN  MULTIPOLYGON (((3 4, 3 3, 2 3, 2 4, 3 4)), ((4...
+        0         1       1.0                POLYGON ((2 2, 2 1, 1 1, 1 2, 2 2))
+        1         2       1.0                POLYGON ((2 2, 2 3, 3 3, 3 2, 2 2))
+        2         2       2.0                POLYGON ((4 4, 4 3, 3 3, 3 4, 4 4))
+        3         1       NaN      POLYGON ((2 0, 0 0, 0 2, 1 2, 1 1, 2 1, 2 0))
+        4         2       NaN  MULTIPOLYGON (((3 4, 3 3, 2 3, 2 4, 3 4)), ((4...
 
         See also
         --------
diff --git a/geopandas/tools/overlay.py b/geopandas/tools/overlay.py
index 6332d425..66338419 100644
--- a/geopandas/tools/overlay.py
+++ b/geopandas/tools/overlay.py
@@ -110,6 +110,21 @@ def _overlay_difference(df1, df2):
     return dfdiff
 
 
+def _overlay_identity(df1, df2):
+    """
+    Overlay Identity operation used in overlay function.
+    """
+    dfintersection = _overlay_intersection(df1, df2)
+    dfdifference = _overlay_difference(df1, df2)
+    dfdifference = _ensure_geometry_column(dfdifference)
+    result = pd.concat([dfintersection, dfdifference], ignore_index=True, sort=False)
+    # keep geometry column last
+    columns = list(dfintersection.columns)
+    columns.remove("geometry")
+    columns.append("geometry")
+    return result.reindex(columns=columns)
+
+
 def _overlay_symmetric_diff(df1, df2):
     """
     Overlay Symmetric Difference operation used in overlay function
@@ -225,11 +240,11 @@ def overlay(df1, df2, how="intersection", keep_geom_type=None, make_valid=True):
 
     >>> geopandas.overlay(df1, df2, how='identity')
        df1_data  df2_data                                           geometry
-    0       1.0       1.0                POLYGON ((2 2, 2 1, 1 1, 1 2, 2 2))
-    1       2.0       1.0                POLYGON ((2 2, 2 3, 3 3, 3 2, 2 2))
-    2       2.0       2.0                POLYGON ((4 4, 4 3, 3 3, 3 4, 4 4))
-    3       1.0       NaN      POLYGON ((2 0, 0 0, 0 2, 1 2, 1 1, 2 1, 2 0))
-    4       2.0       NaN  MULTIPOLYGON (((3 4, 3 3, 2 3, 2 4, 3 4)), ((4...
+    0         1       1.0                POLYGON ((2 2, 2 1, 1 1, 1 2, 2 2))
+    1         2       1.0                POLYGON ((2 2, 2 3, 3 3, 3 2, 2 2))
+    2         2       2.0                POLYGON ((4 4, 4 3, 3 3, 3 4, 4 4))
+    3         1       NaN      POLYGON ((2 0, 0 0, 0 2, 1 2, 1 1, 2 1, 2 0))
+    4         2       NaN  MULTIPOLYGON (((3 4, 3 3, 2 3, 2 4, 3 4)), ((4...
 
     See also
     --------
@@ -333,8 +348,7 @@ def overlay(df1, df2, how="intersection", keep_geom_type=None, make_valid=True):
         elif how == "union":
             result = _overlay_union(df1, df2)
         elif how == "identity":
-            dfunion = _overlay_union(df1, df2)
-            result = dfunion[dfunion["__idx1"].notnull()].copy()
+            result = _overlay_identity(df1, df2)
 
         if how in ["intersection", "symmetric_difference", "union", "identity"]:
             result.drop(["__idx1", "__idx2"], axis=1, inplace=True)


Test patch:
diff --git a/geopandas/tests/test_overlay.py b/geopandas/tests/test_overlay.py
index c9f8dad0..ee35bcc4 100644
--- a/geopandas/tests/test_overlay.py
+++ b/geopandas/tests/test_overlay.py
@@ -90,7 +90,6 @@ def test_overlay(dfs_index, how):
         expected = pd.concat(
             [expected_intersection, expected_difference], ignore_index=True, sort=False
         )
-        expected["col1"] = expected["col1"].astype(float)
     else:
         expected = _read(how)
 
@@ -174,6 +173,10 @@ def test_overlay_nybb(how, nybb_filename):
 
     if how == "identity":
         expected = expected[expected.BoroCode.notnull()].copy()
+        boro_code_dtype = result["BoroCode"].dtype
+        if boro_code_dtype in ("int32", "int64"):
+            # Depending on the pandas version, the dtype might be int32 or int64
+            expected["BoroCode"] = expected["BoroCode"].astype(boro_code_dtype)
 
     # Order GeoDataFrames
     expected = expected.sort_values(cols).reset_index(drop=True)
@@ -786,7 +789,7 @@ def test_non_overlapping(how):
     elif how == "identity":
         expected = GeoDataFrame(
             {
-                "col1": [1.0],
+                "col1": [1],
                 "col2": [np.nan],
                 "geometry": [p1],
             }