[GIS] Getting intersection of multiple polygons efficiently in Python

intersectionpolygonpolygon-creationpythonshapely

I would like to get the intersection of multiple polygons. Using Python's shapely package, I can find the intersection of two polygons using the intersection function. Is there a similar efficient function for obtaining the intersection of multiple polygons?

Here is a code snippet to understand what I mean:

from shapely.geometry import Point

coord1 = ( 0,0 )
point1 = Point(coord1)
circle1 = point1.buffer(1)

coord2 = ( 1,1 )
point2 = Point(coord2)
circle2 = point2.buffer(1)

coord3 = ( 1,0 )
point3 = Point(coord3)
circle3 = point3.buffer(1)

An intersection of two circles can be found by circle1.intersection(circle2). I can find the intersection of all three circles by circle1.intersection(circle2).intersection(circle3). However, this approach is not salable to a large number of polygons as it requires increasingly more code. I would like a function that takes an arbitrary number of polygons and returns their intersection.

Best Answer

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
Related Question