[FIXED] Find overlapping region contained in all polygons within single geodataframe

Issue

I am trying to find a way to display only the overlapping region shared by all polygons in a geodataframe. My geodataframe includes polygon and multipolygon objects from a json file. The implementation I have now involves a loop to iterate over each polygon in the dataframe, however, this is very slow with large number of polygons and also it gives me every single pair of overlaps which is not what I want.

poly1 = Polygon((-119.74583075160622,  35.66006995660339), (-116.30697167993921,35.66006995660339),(-116.30697167993921,34.64381654506288),(-119.74583075160622, 34.64381654506288), (-119.74583075160622,35.66006995660339))
poly2 = Polygon((-118.33337318007834,34.66887933945114),(-117.85950530106018,35.64583276246657),(-117.46804922708861,35.34949851066835),(-116.97357839680875,35.50618853311429),(-116.49971051779052,34.43697449261936),(-116.7881518354538,34.52755087908893),(-118.33337318007834,34.66887933945114))
poly3 = Polygon((-118.73689979855544,34.71741592303824),(-118.0153587441971,35.95953776175309),(-117.45151590111556,34.512388830922866),(-118.27885217295977,34.67014942956098),(-118.73689979855544,34.71741592303824))

above are the contents in my gdf for this example.

overlaps =[]
for index, row in gdf.iterrows():
  for jindex, jrow in gdf.iloc[index+1:].iterrows():
    overlaps.append(jrow["geometry"].intersection(row["geometry"]))

There can be more than 3 overlapping polygons with a shared region, however, It can be assumed each additional polygon will share a common region and not be a outlier.

I would like it to show me the overlapped region within polygons (1 and 2 and 3) only. Not each pair: (1 and 2) (1 and 3) etc...

Is there a way where I can also use sjoin for the comparisons between each polygon since sjoin is faster?

Solution

Without a minimal reproducible example, I can't confirm this, but theoretically you could chain intersections one after the other, using functools.reduce():

import functools
functools.reduce(lambda poly1, poly2: poly1.intersection(poly2), gdf["geometry"])

To make sure I understand the problem correctly, supposing you have polygons A, B, and C -- you want the intersection of all three?

If so, then my suggestion above should work, since intersection is associative, i.e., A ∩ B ∩ C = A ∩ (B ∩ C) = (A ∩ B) ∩ C.

Answered By - ddejohn