Saturday 18 November 2017

python - Polygon overlay with Shapely


I'm trying to capture all the non-overlapping polygons indicated below using Shapely (given polygons A,B & C). Moreover, I'm hoping to do so without iteration, testing for intersect etc. The accepted answer to this question expresses the PostGIS method but it would seem that 'union' means different things to different people.



polygon overlay



Answer



You need to iterate at some level. (Update: I've edited to remove all "for" loops, except for one list comprehension)


# imports used throughout this example
from shapely.geometry import Point
from shapely.ops import cascaded_union
from itertools import combinations

# Here are your input shapes (circles A, B, C)
A = Point(3, 6).buffer(4)

B = Point(6, 2).buffer(4)
C = Point(1, 2).buffer(4)

# list the shapes so they are iterable
shapes = [A, B, C]

First you need the union of all intersections (use a cascaded union), using the combination pair of each shape. Then you remove (via difference) the intersections from the union of all shapes.


# All intersections
inter = cascaded_union([pair[0].intersection(pair[1]) for pair in combinations(shapes, 2)])
# Remove from union of all shapes

nonoverlap = cascaded_union(shapes).difference(inter)

Here is what nonoverlap looks like (via JTS Test Builder): nonoverlap


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...