11

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.

  • im thinking maybe store the coords in a dictionary and loop through it while using from itertools import combinations. I will post soon – ziggy Apr 17 '17 at 15:43
  • What do you mean by "their intersections"? Do you mean all areas that intersect with at least one other polygon, or the areas that all the inputs intersect? – jpmc26 Apr 17 '17 at 21:29
  • I mean the intersection of all polygons, not at least one. – splinter Apr 17 '17 at 21:30
  • You should clarify this above (perhaps with an example output). I'm fairly certain most of the answers don't behave as you wish. (And the fact several answerers have misunderstood is evidence enough that the question needs clarification.) – jpmc26 Apr 17 '17 at 21:36
  • 1
    @jpmc26 I've just added an update to my answer where rtree is used. The approach is more efficient and scalable now. Hope this helps! – Antonio Falciano Apr 18 '17 at 9:06
6

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
  • I don't believe this does what the OP wants. It gives back the areas that at least 2 polygons cover, whereas the OP is only looking for areas covered by all the polygons in the set. See clarification in comments. – jpmc26 Apr 18 '17 at 0:37
3

Why not use a iteration or recursivity? something like :

from shapely.geometry import Point

def intersection(circle1, circle2):
    return circle1.intersection(circle2)

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)
circles = [circle1, circle2, circle3]
intersectionResult = None

for j, circle  in enumerate(circles[:-1]):

    #first loop is 0 & 1
    if j == 0:
        circleA = circle
        circleB = circles[j+1]
     #use the result if the intersection
    else:
        circleA = intersectionResult
        circleB = circles[j+1]
    intersectionResult = intersection(circleA, circleB)

result= intersectionResult
2

Give this code a shot. its pretty simple in concept and I believe gets you what you are looking for.

from shapely.geometry import Point
from itertools import combinations
dic ={}
dic['coord1']=Point(0,0).buffer(1)
dic['coord2']=Point(1,1).buffer(1)
dic['coord3']=Point(1,0).buffer(1)
inter = {k[0]+v[0]:k[1].intersection(v[1]) for k,v in combinations(dic.items(),2)}
print inter

and if you want the output to be stored as a shapefile use fiona:

from shapely.geometry import Point,mapping
import fiona
from itertools import combinations
schema = {'geometry': 'Polygon', 'properties': {'Place': 'str'}}
dic ={}
dic['coord1']=Point(0,0).buffer(1)
dic['coord2']=Point(1,1).buffer(1)
dic['coord3']=Point(1,0).buffer(1)
inter = {k[0]+v[0]:k[1].intersection(v[1]) for k,v in combinations(dic.items(),2)}
print inter
with fiona.open(r'C:\path\abid', "w", "ESRI Shapefile", schema) as output:
    for x,y in inter.items():
        output.write({'properties':{'Place':x},'geometry':mapping(y)})

this outputs -

enter image description here

  • 3
    I don't believe this does what the OP wants. It gives back the areas that at least 2 polygons cover, whereas the OP is only looking for areas covered by all the polygons in the set. See clarification in comments. Additionally, k and v are poor choices for variable names in your dict comprehensions. Those variables each refer to different elements of dic.items(), not a key-value pair. Something like a, b would be less misleading. – jpmc26 Apr 18 '17 at 0:38
  • 1
    ohh okay yeah I did not understand what he meant – ziggy Apr 18 '17 at 0:46
  • and point well taken about my k,v choices--i just automatically use k,v when looping through a dictionary..didnt give it much thought – ziggy Apr 18 '17 at 0:48

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.