Isolating the intersecting polygon where all polygons intersect in Python with the shapely library

Question

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:

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

Answer 1

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

#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")

Answer 2

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)

Answer 3

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: