¿Cuál es la forma más rápida de verificar si un punto está dentro de un polígono en Python?

Encontré dos métodos principales para ver si un punto pertenece a un polígono. Uno usa el método de trazado de rayos que se usa aquí , que es la respuesta más recomendada, el otro usa matplotlib path.contains_points(que me parece un poco oscuro). Tendré que comprobar muchos puntos continuamente. ¿Alguien sabe si alguno de estos dos es más recomendable que el otro o si hay terceras opciones aún mejores?

ACTUALIZAR:

Revisé los dos métodos y matplotlib parece mucho más rápido.

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 = zip(np.random.random(N),np.random.random(N))


# 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)

lo que da,

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

Se obtuvo la misma diferencia relativa usando un triángulo en lugar del polígono de 100 lados. También comprobaré bien proporcionado, ya que parece un paquete dedicado a este tipo de problemas.

Puedes considerar bien formado :

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))

Por los métodos que mencionaste, solo usé el segundo path.contains_points, y funciona bien. En cualquier caso, dependiendo de la precisión que necesite para su prueba, sugeriría crear una cuadrícula bool numpy con todos los nodos dentro del polígono para ser Verdaderos (Falso si no). Si va a hacer una prueba para muchos puntos, esto podría ser más rápido ( aunque tenga en cuenta que esto se basa en que está haciendo una prueba dentro de una tolerancia de "píxeles" ):

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()

, el resultado es este:

Si lo que necesita es velocidad y las dependencias adicionales no son un problema, tal vez lo encuentre numbabastante útil (ahora es bastante fácil de instalar, en cualquier plataforma). El ray_tracingenfoque clásico que propuso se puede migrar fácilmente al numbausar el numba @jitdecorador y convertir el polígono en una matriz numpy. El código debería verse así:

@jit(nopython=True)
def ray_tracing(x,y,poly):
    n = len(poly)
    inside = False
    p2x = 0.0
    p2y = 0.0
    xints = 0.0
    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

La primera ejecución tardará un poco más que cualquier llamada posterior:

%%time
polygon=np.array(polygon)
inside1 = [numba_ray_tracing_method(point[0], point[1], polygon) for 
point in points]

CPU times: user 129 ms, sys: 4.08 ms, total: 133 ms
Wall time: 132 ms

Que, después de la compilación, se reducirá a:

CPU times: user 18.7 ms, sys: 320 µs, total: 19.1 ms
Wall time: 18.4 ms

Si necesita velocidad en la primera llamada de la función, puede compilar previamente el código en un módulo usando pycc. Almacene la función en un src.py como:

from numba import jit
from numba.pycc import CC
cc = CC('nbspatial')


@cc.export('ray_tracing',  'b1(f8, f8, f8[:,:])')
@jit(nopython=True)
def ray_tracing(x,y,poly):
    n = len(poly)
    inside = False
    p2x = 0.0
    p2y = 0.0
    xints = 0.0
    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


if __name__ == "__main__":
    cc.compile()

Constrúyalo python src.pyy ejecútelo:

import nbspatial

import numpy as np
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
# making a list instead of a generator to help debug
points = zip(np.random.random(N),np.random.random(N))

polygon = np.array(polygon)

%%time
result = [nbspatial.ray_tracing(point[0], point[1], polygon) for point in points]

CPU times: user 20.7 ms, sys: 64 µs, total: 20.8 ms
Wall time: 19.9 ms

En el código numba usé: 'b1 (f8, f8, f8 [:,:])'

Para compilar con nopython=True, cada var debe declararse antes de for loop.

En el código src previo a la compilación, la línea:

@cc.export('ray_tracing' , 'b1(f8, f8, f8[:,:])')

Se utiliza para declarar el nombre de la función y sus tipos de var de E / S, una salida booleana b1y dos flotantes f8y una matriz bidimensional de flotantes f8[:,:]como entrada.

Su prueba es buena, pero mide solo una situación específica: tenemos un polígono con muchos vértices y una gran variedad de puntos para verificarlos dentro del polígono.

Además, supongo que no está midiendo matplotlib-inside-polygon-method vs ray-method, sino matplotlib-somehow-optimized-iteration vs simple-list-iteration

Hagamos N comparaciones independientes (N pares de punto y polígono).

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

M = 10000
start_time = time()
# Ray tracing
for i in range(M):
    x,y = np.random.random(), np.random.random()
    inside1 = ray_tracing_method(x,y, polygon)
print "Ray Tracing Elapsed time: " + str(time()-start_time)

# Matplotlib mplPath
start_time = time()
for i in range(M):
    x,y = np.random.random(), np.random.random()
    inside2 = path.contains_points([[x,y]])
print "Matplotlib contains_points Elapsed time: " + str(time()-start_time)

Resultado:

Ray Tracing Elapsed time: 0.548588991165
Matplotlib contains_points Elapsed time: 0.103765010834

Matplotlib sigue siendo mucho mejor, pero no 100 veces mejor. Ahora intentemos un polígono mucho más simple ...

lenpoly = 5
# ... same code

resultado:

Ray Tracing Elapsed time: 0.0727779865265
Matplotlib contains_points Elapsed time: 0.105288982391

Lo dejaré aquí, simplemente reescribí el código anterior usando numpy, tal vez alguien lo encuentre útil:

def ray_tracing_numpy(x,y,poly):
    n = len(poly)
    inside = np.zeros(len(x),np.bool_)
    p2x = 0.0
    p2y = 0.0
    xints = 0.0
    p1x,p1y = poly[0]
    for i in range(n+1):
        p2x,p2y = poly[i % n]
        idx = np.nonzero((y > min(p1y,p2y)) & (y <= max(p1y,p2y)) & (x <= max(p1x,p2x)))[0]
        if p1y != p2y:
            xints = (y[idx]-p1y)*(p2x-p1x)/(p2y-p1y)+p1x
        if p1x == p2x:
            inside[idx] = ~inside[idx]
        else:
            idxx = idx[x[idx] <= xints]
            inside[idxx] = ~inside[idxx]    

        p1x,p1y = p2x,p2y
    return inside    

Ray_tracing envuelto en

def ray_tracing_mult(x,y,poly):
    return [ray_tracing(xi, yi, poly[:-1,:]) for xi,yi in zip(x,y)]

Probado en 100000 puntos, resultados:

ray_tracing_mult 0:00:00.850656
ray_tracing_numpy 0:00:00.003769