Intersect a set polygons/multi-polygons at once and get the intersected polygon

Question

I have a Geopandas data frame with five multi-polygons geometries. Now, I want to intersect all the geometries at once and get the intersected polygon out of it. I tried using unary union and polygonize but this is giving me a list of polygons but I want only one polygon that has the intersection of the set of multi-polygon polygons. How can we intersect a set of multi-polygon or polygons together and get the intersected polygon?

df=
location    geometry    
1          MULTIPOLYGON (((-0.304766 51.425882, -0.304904...    
2          MULTIPOLYGON (((-0.305968 51.427425, -0.30608 ...    
3          MULTIPOLYGON (((-0.358358 51.423471, -0.3581 5...    
4          MULTIPOLYGON (((-0.357654 51.413925, -0.357604...    
    rows=[]
    listpoly = [a.intersection(b) for a, b in combinations(df['geometry'], 2)]
    rings = []
    for poly in listpoly:
      if type(poly) == MultiPolygon:
         exterior_lines = LineString()
         for p in poly:
            exterior_lines = exterior_lines.union(p.exterior)
         rings.append(exterior_lines)
      elif type(poly) == Polygon:
         rings.append(LineString(list(poly.exterior.coords)))
    union = unary_union(rings)
    result = [geom for geom in polygonize(union)]
    print(result)
MULTILINESTRING((-0.0345 54.900...))
MULTILINESTRING((-0.045 54.200...))
MULTILINESTRING((-0.05 54.650...))
MULTILINESTRING((-0.04 54.750...))

Can you add a screenshot showing your data and your expected result — BERA, May 11 '23 at 05:48

score 1 · Accepted Answer · answered May 10 '23 at 21:52

This should work for a set of polygons:

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

Intersect a set polygons/multi-polygons at once and get the intersected polygon

1 Answers1

Set initials

Loop unitil only one polygon left

Linked