1

New to working with shapely and I'm trying to find the intersection of multiple polygons where all (or as many as possible) intersect at a single point to produce a single polygon as below:

Text

The second answer on the below link has pointed me in the right direction however it only produces an empty polygon when I convert the geojson geometry to a shapely geometry.

Getting intersection of multiple polygons efficiently in Python

Unsure where to go from here.

Code for finding intersections with polygons as below:

from shapely.geometry import Point

def intersection(circle1, circle2):
    return circle1.intersection(circle2)

coord1 = ( 0,0 )
point1 = Point(coord1)
circle1 = point1.buffer(1)

coord2 = ( 1,1 )
point2 = Point(coord2)    
circle2 = point2.buffer(1)


coord3 = ( 1,0 )
point3 = Point(coord3)
circle3 = point3.buffer(1)
circles = [circle1, circle2, circle3]
intersectionResult = None

for j, circle  in enumerate(circles[:-1]):

    #first loop is 0 & 1
    if j == 0:
        circleA = circle
        circleB = circles[j+1]
     #use the result if the intersection
    else:
        circleA = intersectionResult
        circleB = circles[j+1]
    intersectionResult = intersection(circleA, circleB)
result= intersectionResult
4
  • Welcome to Geographic Information Systems! Welcome to GIS SE! We're a little different from other sites; this isn't a discussion forum but a Q&A site. Your questions should as much as possible describe not just what you want to do, but precisely what you have tried and where you are stuck trying that. Please check out our short tour for more about how the site works
    – Ian Turton
    Commented Jun 1, 2023 at 12:14
  • can you explain more what doesn't work with the linked solution
    – Ian Turton
    Commented Jun 1, 2023 at 12:16
  • 1
    Does this answer your question? Intersect a set polygons/multi-polygons at once and get the intersected polygon
    – Binx
    Commented Jun 1, 2023 at 15:15
  • Thank you Ian. I've edited accordingly to highlight the specific problem I'm having with the original code
    – Tom
    Commented Jun 1, 2023 at 19:13

3 Answers 3

2

You can use GeoPandas and Shapely to:

1 Extract all polygon boundaries

2 Union them

3 Polygonize

4 Join attributes from the original polygons

import geopandas as gpd
import shapely
import matplotlib.pyplot as plt

df = gpd.read_file(r"/home/bera/Desktop/GIStest/overlaps.geojson")
#df.shape
#(100, 2)

boundaries = df.geometry.boundary.tolist() #A list of all polygon boundaries (as multilinestrings)
boundaries_noded = shapely.unary_union(boundaries) #Union to one multilinestring so vertices are created at each line crossing
singleparts = list(boundaries_noded.geoms) #List all linestrings in the multilinestring to be able to polygonize it
new_polys = shapely.get_parts(shapely.polygonize(singleparts)) #Polygonize creates a geometrycollection, extract all polygons in it

new_polys = gpd.GeoDataFrame(geometry=new_polys, crs=df.crs)

#Create a point df by placing a point in each new polygon
new_points = gpd.GeoDataFrame(geometry=new_polys.geometry.representative_point(), crs=df.crs)
#new_points.shape
#(185, 1)

#Join all attributes from the original df to the points
#new_points.columns
#Index(['geometry'], dtype='object')

new_points_with_attrs = gpd.sjoin(left_df=new_points, right_df=df, how="left")
#new_points_with_attrs.shape
#(295, 3) #From 185 to 295: Points are duplicated where polygons overlap

#new_points_with_attrs.columns
#Out[146]: Index(['geometry', 'index_right', 'id'], dtype='object') #The id column from df was joined

#Then join the point data to the new polygons
new_points_with_attrs = new_points_with_attrs[[col for col in new_points_with_attrs.columns if "index" not in col.lower()]] #Drop the index_right colunm
new_polys_with_attrs = gpd.sjoin(left_df=new_polys, right_df=new_points_with_attrs, how="left") #Or the join will fail

#Calculate number of identical geometries (=overlaps)
new_polys_with_attrs["wkt"] = new_polys_with_attrs.apply(lambda x: x.geometry.wkt, axis=1)
new_polys_with_attrs["overlaps"] = new_polys_with_attrs.groupby('wkt')['wkt'].transform('count')

#There are 102 polygons which do not overlap, 120 which overlap another, 57 two others, and 16 where four polygons overlap
# new_polys_with_attrs["overlaps"].value_counts().sort_index()
# 1    102
# 2    120
# 3     57
# 4     16

fig, ax = plt.subplots(figsize=(15, 15))
new_polys_with_attrs.plot(ax=ax, column="overlaps", legend=True, cmap='YlOrRd')

enter image description here

#You can extract the polygons with different number of overlaps with
overlaps_4 = new_polys_with_attrs.loc[new_polys_with_attrs["overlaps"]==4]
overlaps_4.plot(color="red")

enter image description here

0

Here is my answer from this question. Should work the same for you.

import os
os.environ['USE_PYGEOS'] = '0'  # this is required by shapely when using python 3.10
import geopandas


file = geopandas.read_file("C:/Users/user/Desktop/all_multi.shp")
the_crs = file.crs

# Set initials
intersect_list = []
uniq = []
any_intersect = []

# Loop unitil only one polygon left
while True:
    intersect_list = []
    for i, geom in enumerate(file.geometry):
        for j, next_geom in enumerate(file.geometry):
            if i != j:
                if geom.intersects(next_geom):
                    intersect_list.append(geom.intersection(next_geom))
                                      
    # Set last unique list
    last = uniq

    # Remove duplicate geometries
    uniq = []
    for poly in intersect_list:
        if not any(p.equals(poly) for p in uniq):
            uniq.append(poly)    
    
    # Check if there are anymore intersections
    if len(uniq) == 0:
        break

    # Update GeoDataFrame
    file = geopandas.GeoDataFrame(geometry=uniq)

    
geopandas.GeoDataFrame(geometry=last).to_file("C:/Users/user/Desktop/intersection.shp", crs=the_crs)
14
  • Thank you Binx. I've run your code using my geojson files which have been converted to geopandas, however when I run the code the Jupyter Notebook hangs with a [*] next to the cell. Is this code quite resource intensive?
    – Tom
    Commented Jun 1, 2023 at 19:16
  • I have not tested the efficiency of it. But at first glance, if you had a ton of polygons, yes it would most likely take a while due to those two for loops. How many polygons do you have?
    – Binx
    Commented Jun 1, 2023 at 19:21
  • I have a total of 63 polygons spread across 6 multi-polygons
    – Tom
    Commented Jun 1, 2023 at 19:44
  • hmmm I wouldn't think it would take that long. Would you mind sharing your dataset?
    – Binx
    Commented Jun 1, 2023 at 21:01
  • 1
    Hi Binx - no apology needed, very grateful for help. Thanks for the update, I've not got a specific deadline for this. I've reviewed a few other SO pages however had little to no success. I'd be very interested in your findings!
    – Tom
    Commented Jun 12, 2023 at 15:12
0

I think the following code snippet does what you are looking for, with file "intersection_all.shp.zip" containing the shapefile you shared:

from pathlib import Path

import geopandas as gpd
from matplotlib import pyplot as plt
import shapely
import shapely.plotting

script_dir = Path(__file__).resolve().parent
data_gdf = gpd.read_file(script_dir / "intersection_all.shp.zip")
result = shapely.intersection_all(data_gdf.geometry)

data_gdf.plot(alpha=0.50)
shapely.plotting.plot_polygon(result, color="red", add_points=False)
plt.show()

Result, with your original input in blue and the areas where everything intersects in red:

enter image description here

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.