Implementing polygon self-intersection in Python

Question

Is there a Python library and function that would allow me to carry out polygon self-intersection on a single polygons layer, such as what is achieved by using the “SAGA Polygon self-intersection” tool in QGIS?

Every other python function I am finding for polygon intersections requires two separate layers to be intersected, rather than simply identifying the intersection between polygons in a single layer. This means that after running a polygons layer through this “polygon self-intersection” python function I am looking for, the geodataframe (if using geopandas) would then have additional rows added to it for additionally added spaces, the now identified intersection zones. There would as well involve the “ID” column including the IDs of all polygons that “make up” each intersection zone. For instance, polygon 1 and polygon 2 non overlapped spaces would have an “ID” of “1” and “2” respectively, and the overlapped region between these two polygons would have an “ID” of “1 | 2”, just like would be shown in QGIS after using the “SAGA Polygon self-intersection” tool. And as well, if we had a triple intersection with a polygon 3, the newly identified overlapped region between these three polygons would have an “ID” of “1 | 2 | 3”, just like would shown in QGIS after using the “SAGA Polygon self-intersection” tool.

I have not been able to find this within GeoPandas, as the GeoPandas polygon overlay functions appear to require two separate layers as inputs.

Here is what my shapefile looks like:

For reference: The polygons I used were turned into their centroids and then buffered to 1/4 mile distance.

Here are the centroid coordinates of my buffers:

0     POINT (-13604211.017 4554804.595)
1     POINT (-13610727.382 4559245.985)
2     POINT (-13614171.974 4565380.476)
3     POINT (-13611036.059 4562595.117)
4     POINT (-13611986.274 4564049.187)
5     POINT (-13612916.586 4563936.448)
6     POINT (-13602234.338 4563015.607)
7     POINT (-13612399.115 4562654.820)
8     POINT (-13609971.081 4559031.424)
9     POINT (-13613293.685 4558725.535)
10    POINT (-13611922.699 4565876.467)
11    POINT (-13613733.483 4556423.327)
12    POINT (-13604485.660 4555089.412)
13    POINT (-13605022.256 4551148.583)
14    POINT (-13603752.634 4550764.525)
15    POINT (-13611500.291 4548909.827)
16    POINT (-13606620.852 4547777.264)
17    POINT (-13604949.018 4551387.605)
18    POINT (-13610794.814 4553121.618)
19    POINT (-13610077.870 4553680.083)
20    POINT (-13615078.994 4566347.140)
21    POINT (-13609346.498 4550470.678)
22    POINT (-13609117.309 4554812.264)
23    POINT (-13609179.017 4551117.796)
24    POINT (-13603320.178 4550288.519)
25    POINT (-13611568.431 4547933.755)
26    POINT (-13601285.925 4550978.405)
27    POINT (-13608327.200 4553650.316)
28    POINT (-13607025.788 4546911.740)
29    POINT (-13625355.222 4550774.185)
30    POINT (-13614118.477 4562602.638)
31    POINT (-13611030.521 4562315.516)
32    POINT (-13604582.190 4547459.957)
33    POINT (-13610192.502 4553936.027)
34    POINT (-13602412.265 4547187.654)
35    POINT (-13605394.544 4549989.133)
36    POINT (-13604610.731 4552881.663)
37    POINT (-13600146.915 4543612.151)
38    POINT (-13600058.162 4543634.806)
39    POINT (-13605164.136 4547391.426)
40    POINT (-13604729.135 4547419.234)
41    POINT (-13602534.183 4547255.310)
42    POINT (-13612207.081 4549163.977)
43    POINT (-13603950.618 4548378.891)
44    POINT (-13607061.377 4549046.392)
45    POINT (-13602507.710 4547248.595)
46    POINT (-13602514.041 4547218.106)
47    POINT (-13609998.343 4553042.678)
dtype: geometry

Answer 1

You can cross join the df to itself, then check for intersections and intersect.

This can be slow, or wont work with large dataframes. For example 1000 rows will become 1 million with the cross join.

import geopandas as gpd
import numpy as np

df['id'] = np.arange(df.shape[0]) #Create an id column

#df.shape
#Out[48]: (48, 3)

cross = pd.merge(left=df, right=df, how='cross') #Cross join the dataframe to itself
#cross.shape
#Out[50]: (2304, 6)

cross = cross.loc[cross['id_x'] < cross['id_y']] #Remove self joins
#cross.shape
#Out[53]: (1128, 6)

cross = cross.loc[cross.geometry_x.intersects(cross.geometry_y)] #Select only polygons intersecting
#cross.shape
#Out[55]: (138, 6)

cross['inter'] = cross.geometry_x.intersection(cross.geometry_y) #Intersect them

cross = cross.set_geometry('inter') #The error was here. I thought set_geometry was inplace, but it is returning the dataframe
#cross.shape
#Out[60]: (138, 7)

cross = cross[[column for column in cross.columns if 'geom' not in column]] #Drop old geometries
cross.to_file(r'/home/bera/Desktop/GIStest/1milebuffers_overlaps.shp')

Input purple, output black