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
    Can't you just give the overlay method the same layer as the second layer? Even if you can't, you could do a deep copy of the layer and give that as the second input.
    – Llaves
    Commented Sep 1, 2022 at 23:25
  • Look at Remove the intersecting parts in a set of polygons
    – gene
    Commented Sep 2, 2022 at 7:21

3 Answers 3

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. Commented 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. Commented Sep 12, 2022 at 19:11
  • I just added the centroid coordinates to my post at the end. Commented 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 Commented Sep 19, 2022 at 17:12
  • 1
    Try groupby by geometry and list on ids
    – Bera
    Commented Sep 21, 2022 at 5:13
0

I think this should give the desired result when more than two polygons can intersect each other:

import geopandas as gpd
import pandas as pd
from shapely import unary_union, from_wkt
from shapely.geometry import MultiPolygon

center_wkts = ['POINT (-13604211.017 4554804.595)',
                'POINT (-13610727.382 4559245.985)',
                'POINT (-13614171.974 4565380.476)',
                'POINT (-13611036.059 4562595.117)',
                'POINT (-13611986.274 4564049.187)',
                'POINT (-13612916.586 4563936.448)',
                'POINT (-13602234.338 4563015.607)',
                'POINT (-13612399.115 4562654.820)',
                'POINT (-13609971.081 4559031.424)',
                'POINT (-13613293.685 4558725.535)',
                'POINT (-13611922.699 4565876.467)',
                'POINT (-13613733.483 4556423.327)',
                'POINT (-13604485.660 4555089.412)',
                'POINT (-13605022.256 4551148.583)',
                'POINT (-13603752.634 4550764.525)',
                'POINT (-13611500.291 4548909.827)',
                'POINT (-13606620.852 4547777.264)',
                'POINT (-13604949.018 4551387.605)',
                'POINT (-13610794.814 4553121.618)',
                'POINT (-13610077.870 4553680.083)',
                'POINT (-13615078.994 4566347.140)',
                'POINT (-13609346.498 4550470.678)',
                'POINT (-13609117.309 4554812.264)',
                'POINT (-13609179.017 4551117.796)',
                'POINT (-13603320.178 4550288.519)',
                'POINT (-13611568.431 4547933.755)',
                'POINT (-13601285.925 4550978.405)',
                'POINT (-13608327.200 4553650.316)',
                'POINT (-13607025.788 4546911.740)',
                'POINT (-13625355.222 4550774.185)',
                'POINT (-13614118.477 4562602.638)',
                'POINT (-13611030.521 4562315.516)',
                'POINT (-13604582.190 4547459.957)',
                'POINT (-13610192.502 4553936.027)',
                'POINT (-13602412.265 4547187.654)',
                'POINT (-13605394.544 4549989.133)',
                'POINT (-13604610.731 4552881.663)',
                'POINT (-13600146.915 4543612.151)',
                'POINT (-13600058.162 4543634.806)',
                'POINT (-13605164.136 4547391.426)',
                'POINT (-13604729.135 4547419.234)',
                'POINT (-13602534.183 4547255.310)',
                'POINT (-13612207.081 4549163.977)',
                'POINT (-13603950.618 4548378.891)',
                'POINT (-13607061.377 4549046.392)',
                'POINT (-13602507.710 4547248.595)',
                'POINT (-13602514.041 4547218.106)',
                'POINT (-13609998.343 4553042.678)']

gdf = gpd.GeoDataFrame(geometry=[from_wkt(x).buffer(1609.344/4) for x in center_wkts], crs='epsg:3857')

gdf['parents'] = [{str(x)} for x in gdf.index] # Use index as IDs of polygons

while True:
    
    # Need to copy geometry column because spatial join does not retain it
    # https://gis.stackexchange.com/questions/393387/keeping-geometry-column-from-both-dataframes-when-applying-sjoin-using-geopand
    gdf['geometry_copy'] = gdf.geometry
    
    # Use spatial self join (can leverage spatial index)
    joined_gdf = gpd.sjoin(gdf, gdf, how='inner', predicate='intersects')
    joined_gdf = joined_gdf.drop(columns='geometry_copy_left') # After joining the copy is not needed anymore
    joined_gdf = joined_gdf.rename(columns={'geometry_copy_right': 'geometry_right'})
    joined_gdf = joined_gdf[~(joined_gdf['parents_right'] <= joined_gdf['parents_left'])] # Right parents not subset of left partents and vice versa
    joined_gdf = joined_gdf[~(joined_gdf['parents_left'] <= joined_gdf['parents_right'])] 
    joined_gdf['parents'] = joined_gdf.apply(lambda x: x['parents_left'].union(x['parents_right']), axis=1) # Use set union to determine new unique parent IDs
    
    difference_gdf = joined_gdf.copy() # Copy in order to calc difference on each intersecting polygon later
    
    without_intersection_gdf = gdf[~gdf.index.isin(joined_gdf.index)].drop(columns='geometry_copy') # These polygons don't intersect with any other polygon
    
    # Intersections of polygons
    joined_gdf = joined_gdf.drop(columns=['parents_left', 'parents_right'])
    joined_gdf['geometry'] = joined_gdf.geometry.intersection(joined_gdf['geometry_right'])
    joined_gdf = joined_gdf[joined_gdf.geometry.type.isin(['Polygon', 'MultiPolygon'])] # Drop all other geom types
    
    # Stop if there are no intersecting polygons
    if joined_gdf.empty:
        gdf = gdf.drop(columns='geometry_copy')
        break
    
    joined_gdf = joined_gdf.drop(columns=['index_right', 'geometry_right'])
    
    # For each polygon with intersection(s), subtract the union of all polygons that intersect it
    difference_gdf = difference_gdf.drop(columns=['parents', 'parents_right']).rename(columns={'parents_left': 'parents'})
    difference_gdf = difference_gdf.groupby(difference_gdf.index).agg({'geometry': 'first', 'parents': 'first', 'geometry_right': unary_union})
    difference_gdf.geometry = difference_gdf['geometry'].difference(difference_gdf['geometry_right'])
    difference_gdf.crs = gdf.crs
    difference_gdf = difference_gdf[~difference_gdf.geometry.is_empty] # Drop empty geometries (in case a polygon completely covers another polygon)
    difference_gdf = difference_gdf.drop(columns='geometry_right')
    
    # Finally, concatenate the three intermediate results
    gdf = pd.concat([without_intersection_gdf, joined_gdf, difference_gdf])
    gdf = gdf.drop_duplicates(subset='parents')
    gdf = gdf.reset_index(drop=True)

gdf['ID'] = gdf['parents'].apply(lambda x: '|'.join(sorted(x)))
gdf = gdf.drop(columns='parents')

# If single polygons and multipolygons are mixed, turn all polygons into multipolygons
# (Required by GeoPackage spec)
if set(gdf.geometry.type.unique()) == {'MultiPolygon', 'Polygon'}:
    gdf.geometry = gdf.geometry.apply(lambda x: MultiPolygon([x]) if x.geom_type=='Polygon' else x)
gdf.to_file('out.gpkg')

In simple terms, this works by iteratively performing these 3 steps:

  • Calculate intersections between polygons
  • Detect polygons that do not intersect any other polygons
  • For a polygon that intersects one or more other polygons, calculate the difference between the polygon and the union of intersecting polygons

In each iteration the results of these steps are concatenated while keeping track of parent polygons and removing duplicates.

EDIT

I found that while the above code works for the data provided by the OP, it fails in complicated scenarios where a lot of polygons intersect each other, as it produces huge intermediate dataframes. The following code should perform better in such cases. It works in a similar way, but starts with a single feature and joins one feature in each iteration.

gdf['parents'] = [{str(x)} for x in gdf.index]

gdf_to_join = gdf.iloc[1:]
gdf = gdf.head(1)

for feat in gdf_to_join.iterfeatures():

    feat_to_join_gdf = gpd.GeoDataFrame.from_features([feat], crs=gdf.crs)
    feat_to_join_gdf['geometry_copy'] = feat_to_join_gdf.geometry
    joined_gdf = gpd.sjoin(gdf,
                           feat_to_join_gdf,
                           how='left', predicate='intersects')

    feat_to_join_gdf.drop(columns='geometry_copy', inplace=True)
    
    joined_gdf = joined_gdf.rename(columns={'geometry_copy': 'geometry_right'})

    intersect_gdf = joined_gdf[~joined_gdf['index_right'].isna()].copy()
    
    if intersect_gdf.empty:
        gdf = pd.concat([gdf, feat_to_join_gdf])
    else:
        wo_inter_gdf = joined_gdf[joined_gdf['index_right'].isna()]
        wo_inter_gdf = wo_inter_gdf.copy()
        wo_inter_gdf.rename(columns={'parents_left': 'parents'}, inplace=True)
        wo_inter_gdf.drop(columns=['index_right',
                                   'parents_right',
                                   'geometry_right'], inplace=True)

        diff_gdf = intersect_gdf.copy()
        diff_gdf['geometry'] = diff_gdf.geometry.difference(
            diff_gdf['geometry_right'])

        diff_gdf.rename(columns={'parents_left': 'parents'}, inplace=True)
        diff_gdf.drop(columns=['index_right',
                               'parents_right',
                               'geometry_right'], inplace=True)

        feat_to_join_gdf.geometry = feat_to_join_gdf.difference(
            unary_union(intersect_gdf.geometry))
        
        intersect_gdf['parents'] = intersect_gdf.apply(
            lambda x: x['parents_left'].union(x['parents_right']), axis=1)

        intersect_gdf.geometry = intersect_gdf.geometry.intersection(
            intersect_gdf['geometry_right'])
        
        intersect_gdf.drop(columns=['parents_left',
                                    'index_right',
                                    'parents_right',
                                    'geometry_right'], inplace=True)

        gdf = pd.concat([wo_inter_gdf,
                         diff_gdf,
                         feat_to_join_gdf,
                         intersect_gdf])

        gdf = gdf[gdf.geometry.type.isin(['Polygon', 'MultiPolygon'])]
        gdf = gdf[~gdf.geometry.is_empty]

gdf['ID'] = gdf['parents'].apply(lambda x: '|'.join(sorted(x)))
gdf = gdf.drop(columns='parents')

if set(gdf.geometry.type.unique()) == {'MultiPolygon', 'Polygon'}:
    gdf.geometry = gdf.geometry.apply(lambda x: MultiPolygon([x])
                                      if x.geom_type=='Polygon' else x)      

gdf.to_file('out.gpkg')
0

Recently, support for self-overlays was added to geofileops. The geofilops.union tool gives a somewhat similar output to the saga tool based on your description, but it's not the same. The result is based on "join" operation, so overlapping areas will be split of of the original polygons but will be in the output file twice, like in the QGIS union tool.

Performance should be good, especially for large input files, as processing is parallelized.

The input is expected to be in a Geopackages file.

Sample script:

from pathlib import Path
import geofileops as gfo
import logging

if __name__ == "__main__":
    # Init logging to have progress logging
    logging.basicConfig(level=logging.INFO)

    input_path = Path("C:/Temp/input.gpkg")
    output_path = input_path.parent / "output.gpkg"

    gfo.union(input1_path=input_path, input2_path=None, output_path=output_path)

Disclaimer: I'm the developer of geofileops

Your Answer

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

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