pythongeometry

What's the fastest way of checking if a point is inside a polygon in python


I found two main methods to look if a point belongs inside a polygon. One is using the ray tracing method used here, which is the most recommended answer, the other is using matplotlib path.contains_points (which seems a bit obscure to me). I will have to check lots of points continuously. Does anybody know if any of these two is more recommendable than the other or if there are even better third options?

I checked the two methods and matplotlib looks much faster.

from time import time
import numpy as np
import matplotlib.path as mpltPath

# regular polygon for testing
lenpoly = 100
polygon = [[np.sin(x)+0.5,np.cos(x)+0.5] for x in np.linspace(0,2*np.pi,lenpoly)[:-1]]

# random points set of points to test 
N = 10000
points = np.random.rand(N,2)


# Ray tracing
def ray_tracing_method(x,y,poly):

    n = len(poly)
    inside = False

    p1x,p1y = poly[0]
    for i in range(n+1):
        p2x,p2y = poly[i % n]
        if y > min(p1y,p2y):
            if y <= max(p1y,p2y):
                if x <= max(p1x,p2x):
                    if p1y != p2y:
                        xints = (y-p1y)*(p2x-p1x)/(p2y-p1y)+p1x
                    if p1x == p2x or x <= xints:
                        inside = not inside
        p1x,p1y = p2x,p2y

    return inside

start_time = time()
inside1 = [ray_tracing_method(point[0], point[1], polygon) for point in points]
print("Ray Tracing Elapsed time: " + str(time()-start_time))

# Matplotlib mplPath
start_time = time()
path = mpltPath.Path(polygon)
inside2 = path.contains_points(points)
print("Matplotlib contains_points Elapsed time: " + str(time()-start_time))

which gives,

Ray Tracing Elapsed time: 0.441395998001
Matplotlib contains_points Elapsed time: 0.00994491577148

Same relative difference was obtained one using a triangle instead of the 100 sides polygon. I will also check shapely since it looks a package just devoted to these kind of problems.


Solution

  • You can consider shapely:

    from shapely.geometry import Point
    from shapely.geometry.polygon import Polygon
    
    point = Point(0.5, 0.5)
    polygon = Polygon([(0, 0), (0, 1), (1, 1), (1, 0)])
    print(polygon.contains(point))
    

    From the methods you've mentioned I've only used the second, path.contains_points, and it works fine. In any case depending on the precision you need for your test I would suggest creating a numpy bool grid with all nodes inside the polygon to be True (False if not). If you are going to make a test for a lot of points this might be faster (although notice this relies you are making a test within a "pixel" tolerance):

    from matplotlib import path
    import matplotlib.pyplot as plt
    import numpy as np
    
    first = -3
    size  = (3-first)/100
    xv,yv = np.meshgrid(np.linspace(-3,3,100),np.linspace(-3,3,100))
    p = path.Path([(0,0), (0, 1), (1, 1), (1, 0)])  # square with legs length 1 and bottom left corner at the origin
    flags = p.contains_points(np.hstack((xv.flatten()[:,np.newaxis],yv.flatten()[:,np.newaxis])))
    grid = np.zeros((101,101),dtype='bool')
    grid[((xv.flatten()-first)/size).astype('int'),((yv.flatten()-first)/size).astype('int')] = flags
    
    xi,yi = np.random.randint(-300,300,100)/100,np.random.randint(-300,300,100)/100
    vflag = grid[((xi-first)/size).astype('int'),((yi-first)/size).astype('int')]
    plt.imshow(grid.T,origin='lower',interpolation='nearest',cmap='binary')
    plt.scatter(((xi-first)/size).astype('int'),((yi-first)/size).astype('int'),c=vflag,cmap='Greens',s=90)
    plt.show()
    

    The result is this:

    point inside polygon within pixel tolerance