3

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.

3

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

enter image description here

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

enter image description here

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.