Getting Unique Intersections of Multiple Polygons Efficiently in Python

Question

I would like to get a list of all polygons that are intersections of other polygons. However, that list of intersections cannot self-overlap. For example, if I have a triple Venn diagram, I would want four polygon intersections returned: One for the center where the three circles intersect; and one each for the intersections of pairs of circles, where those three intersections do not include the center intersection.

Based on the top answer to this question, I have this code:

from shapely.geometry import Polygon
from shapely.geometry import MultiPolygon
from rtree import index

#read in my polygons separately, not shown, store to array of Polygons

intersections = []
idx = index.Index()
for pos, polygon in enumerate(shapely_arr_polygon):
    idx.insert(pos, polygon.bounds)

for polygon in shapely_arr_polygon:
    merged_polygons = cascaded_union([shapely_arr_polygon[pos] for pos in idx.intersection(polygon.bounds) if shapely_arr_polygon[pos] != polygon])
    intersections.append(polygon.intersection(merged_polygons))

intersection = cascaded_union(intersections)

From this code block, I get "intersections" which is a list of Polygons that are all pair-wise intersections (not what I want), and I get "intersection" which is a MultiPolygon of the union of all of those intersections (not what I want).

I'm thinking there should be a way to use either/both of the above to do what I want, or perhaps a different way to get what I want in that loop, but I can't figure it out ... help?

Full ideal solution: I would have a list of polygon intersections that best represent the overlaps of the original polygons, such that in the triple Venn diagram case, 4 overlaps with that center being separate because it's 3 overlaps instead of 2.

Maybe good enough solution: Unique overlaps that just don't overlap with each other. So in the triple Venn diagram case, it could be good enough to just return 3 overlaps, where the center is included in ONE of the pair-wise overlaps. I think that could avoid a cominatorial situation.

Answer 1

Look at the How to find the intersection areas of overlapping buffer zones in single shapefile? and the solution is to
1) use the properties of unary_union with the LinearRings/LineString of the polygons (it cuts the lines at each intersection,Planar graph)

listpoly = [a.intersection(b) for a, b in combinations(shapely_arr_polygon, 2)]
print(len(listpoly)
3 

The result is 3 overlapping polygons

Now use the unary_union predicate

rings = [LineString(list(pol.exterior.coords)) for pol in listpoly]
from shapely.ops import unary_union, polygonize
union = unary_union(rings)

2) then use the polygonize function

result = [geom for geom in polygonize(union)]
print(len(result))
4 

The result is 4 polygons that do not overlap