3

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: enter image description here

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
2

1 Answer 1

3

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')

enter image description here Input purple, output black

11
  • Thank you for this very thorough suggestion, this very much looks like what I am going for. And yes, to clarify, I mean polygons in the same layer intersecting, and not one polygon intersecting itself. Though in trying out your code suggestion, I am a bit confused about the output dataframe, since when I took a look at the final geodataframe "cross" at the end of your code, I was left with just 2 rows, but I am not sure why I ended up with only two rows, if when referring to my visual representation of my polygons I just added above, I have what I see is over a dozen intersection zones. Sep 8, 2022 at 22:53
  • No, just polygons, where each row in the geodataframe corresponds to a unique, separate, individual polygon. In my case, each of these unique polygons are just distance buffers. Sep 12, 2022 at 19:11
  • I just added the centroid coordinates to my post at the end. Sep 13, 2022 at 19:55
  • The crs is "EPSG:3857". And here is the full properties of my polygons shapefile: <Derived Projected CRS: EPSG:3857> Name: WGS 84 / Pseudo-Mercator Axis Info [cartesian]: - X[east]: Easting (metre) - Y[north]: Northing (metre) Area of Use: - name: World between 85.06°S and 85.06°N. - bounds: (-180.0, -85.06, 180.0, 85.06) Coordinate Operation: - name: Popular Visualisation Pseudo-Mercator - method: Popular Visualisation Pseudo Mercator Datum: World Geodetic System 1984 ensemble - Ellipsoid: WGS 84 - Prime Meridian: Greenwich Sep 19, 2022 at 17:12
  • 1
    Try groupby by geometry and list on ids
    – BERA
    Sep 21, 2022 at 5:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.