36

I'm looking for an algorithm, a high-level solution, or even a library which can help me determine if two polygons intersect, in Python.

I have the vertices of the two polygons (These are single-part polygons without any holes) in two different arrays. The polygons are 2D (i.e. just X and Y coordinates)

I'll like to make a function which will return a boolean indicating whether these two polygons intersect.

Please note that I cannot use arcpy, or any arcgis components in this.

Can you suggest an algorithm or library for doing this?

6 Answers 6

66

You could try shapely.

They describe spatial relationships and it works on Windows

The spatial data model is accompanied by a group of natural language relationships between geometric objects – contains, intersects, overlaps, touches, etc. – and a theoretical framework for understanding them using the 3x3 matrix of the mutual intersections of their component point sets

The following code shows how you can test for intersection:

from shapely.geometry import Polygon

p1 = Polygon([(0,0), (1,1), (1,0)])
p2 = Polygon([(0,1), (1,0), (1,1)])
print(p1.intersects(p2))
26

You can use the GDAL/OGR Python bindings for that.

from osgeo import ogr

wkt1 = "POLYGON ((1208064.271243039 624154.6783778917, 1208064.271243039 601260.9785661874, 1231345.9998651114 601260.9785661874, 1231345.9998651114 624154.6783778917, 1208064.271243039 624154.6783778917))"
wkt2 = "POLYGON ((1199915.6662253144 633079.3410163528, 1199915.6662253144 614453.958118695, 1219317.1067437078 614453.958118695, 1219317.1067437078 633079.3410163528, 1199915.6662253144 633079.3410163528)))"

poly1 = ogr.CreateGeometryFromWkt(wkt1)
poly2 = ogr.CreateGeometryFromWkt(wkt2)

intersection = poly1.Intersection(poly2)

print intersection.ExportToWkt()

It returns None if they don't intersect. If they intersect it returns the geometry where both intersect.

Also, you can find further info in the GDAL/OGR Cookbook.

3
  • I would love to use this, but I'm on windows, and on both the systems i've tried, I can't get the python bindings to work. I run into the problem described in this post: gis.stackexchange.com/questions/44958/… Commented Mar 18, 2014 at 13:13
  • 1
    Just in case someone else stumbles upon this, it is possible to use GDAL/OGR with Python in Windows (and within ArcGIS, no less): gis.stackexchange.com/questions/74524/… Commented Mar 19, 2014 at 12:16
  • you also can write intersection = poly1.Intersect(poly2) --- the value of intersection will be TRUE or FALSE depending on if the polygons intersect
    – Max
    Commented Nov 19, 2015 at 10:01
4

I know this is an old question, but I've written a python library for handling collisions between concave and convex polygons, as well as circles.

It's pretty simple to use, here you go!

Example:

from collision import *
from collision import Vector as v

p0 = Concave_Poly(v(0,0), [v(-80,0), v(-20,20), v(0,80), v(20,20), v(80,0),  v(20,-20), v(0,-80), v(-20,-20)])
p1 = Concave_Poly(v(20,20), [v(-80,0), v(-20,20), v(0,80), v(20,20), v(80,0),  v(20,-20), v(0,-80), v(-20,-20)])

print(collide(p0,p1))

You can also have it generate a reponse, which includes:

overlap (how much they overlap)
overlap vector (when subtracted from second shapes position, the shapes will no longer be colliding)
overlap vector normalized (vector direction of collision)
a in b (whether the first shape is fully inside the second)
b in a (whether the second shape is fully inside the first)

https://github.com/QwekoDev/collision

4

You can also use PyQGIS.

from qgis.core import QgsGeometry

p1 = QgsGeometry.fromWkt('POLYGON ((0 0, 1 1, 1 0, 0 0))')
p2 = QgsGeometry.fromWkt('POLYGON ((0 1, 1 0, 1 1, 0 1))')

intersection = p1.intersection(p2)

print(intersection)

# OUT
# <QgsGeometry: Polygon ((1 1, 1 0, 0.5 0.5, 1 1))>
3

If you want to know the level you can use this. As an argument, you can give a list of polygons. And as a return value, you get a list of levels. In the list of levels, there are the polygons.

from shapely.geometry import Point
from shapely.geometry.polygon import Polygon

def isPolygonInPolygon(poly1,poly2):
    poly2 = Polygon(poly2)
    for poi in poly1:
        poi = Point(poi)
        if(poly2.contains(poi)):
            return True

def polygonTransformHierarchy(polygon_list):
    polygon_list_hierarchy = []
    for polygon1 in polygon_list:
        level = 0
        for polygon2 in polygon_list:
            if(isPolygonInPolygon(polygon1, polygon2)):
                level += 1
        if(level > len(polygon_list_hierarchy)-1):
            dif = (level+1)- len(polygon_list_hierarchy)
            for _ in range(dif):
                polygon_list_hierarchy.append([])   
        polygon_list_hierarchy[level].append(polygon1)
    return polygon_list_hierarchy
1

If you know or are interested in learning R it has some useful spatial packages. http://cran.r-project.org/web/views/Spatial.html There is Python module to interact with R (RPy*)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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