Python Shapely Intersections – Getting Unique Intersections of Multiple Polygons Efficiently in Python

intersectionpythonshapely

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.

Best Answer

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

Related Solutions

Python – Efficient Methods for Identifying Overlapping Polygons

First, welcome to the site!

Without knowing what other GIS software you've got available I'd suggest looking at a spatially enabled database, such as PostGIS to store your data, and accessing the data with GeoAlchemy in Python (I like SpatailLite as a great lightweight alternative, but access in Python is a bit painful as you have to recompile the sqlite3 library from scratch).

The database will automatically build an RTree index which stores extents of objects. It then queries the extents of objects to see if they intersect before querying the actual geometries to significantly increase the speed of the query.

If you don't want to use a database you can build stand-alone RTree indexes in Python with the Rtree module. Construct an RTree from the bounds of your US County polygons and then compare that to the bounds of your incoming polygons to narrow down your query region. Then do a spatial intersection of the counties with your actual geometry to get whether or not they intersect. If you're having trouble constructing geometries I'd also strongly recommend looking at the Shapely library.

[GIS] Getting intersection of multiple polygons efficiently in Python

One possible approach could be considering the combination of pairs of polygons, their intersections and finally the union of all the intersections via a cascaded union (like suggested here):

from shapely.geometry import Point
from shapely.ops import cascaded_union
from itertools import combinations

circles = [
    Point(0,0).buffer(1),
    Point(1,0).buffer(1),
    Point(1,1).buffer(1),
]

intersection = cascaded_union(
    [a.intersection(b) for a, b in combinations(circles, 2)]
)
print intersection

A more efficient approach should use a spatial index, like Rtree, in order to deal with a lot of geometries (not the case of the three circles):

from shapely.geometry import Point
from shapely.ops import cascaded_union
from rtree import index

circles = [
    Point(0,0).buffer(1),
    Point(1,0).buffer(1),
    Point(1,1).buffer(1),
]
intersections = []
idx = index.Index()

for pos, circle in enumerate(circles):
    idx.insert(pos, circle.bounds)

for circle in circles:
    merged_circles = cascaded_union([circles[pos] for pos in idx.intersection(circle.bounds) if circles[pos] != circle])
    intersections.append(circle.intersection(merged_circles))

intersection = cascaded_union(intersections)
print intersection

Best Answer

Related Solutions

Python – Efficient Methods for Identifying Overlapping Polygons

[GIS] Getting intersection of multiple polygons efficiently in Python

Related Question