Was ist der schnellste Weg, um zu überprüfen, ob sich ein Punkt in einem Polygon in Python befindet?

Ich habe zwei Hauptmethoden gefunden, um zu prüfen, ob ein Punkt in ein Polygon gehört. Der eine verwendet die hier verwendete Raytracing-Methode , die die am meisten empfohlene Antwort ist, der andere verwendet Matplotlib path.contains_points(was mir etwas dunkel erscheint). Ich werde ständig viele Punkte überprüfen müssen. Weiß jemand, ob eine dieser beiden Optionen empfehlenswerter ist als die andere oder ob es noch bessere dritte Optionen gibt?

AKTUALISIEREN:

Ich habe die beiden Methoden überprüft und matplotlib sieht viel schneller aus.

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)

was gibt,

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

Der gleiche relative Unterschied wurde unter Verwendung eines Dreiecks anstelle des 100-Seiten-Polygons erhalten. Ich werde auch formschön prüfen, da es sich um ein Paket handelt, das nur diesen Problemen gewidmet ist

Sie können formschön betrachten :

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

Von den Methoden, die Sie erwähnt haben, habe ich nur die zweite verwendet path.contains_points, und es funktioniert gut. In jedem Fall würde ich vorschlagen, abhängig von der Genauigkeit, die Sie für Ihren Test benötigen, ein Numpy-Bool-Gitter zu erstellen, bei dem alle Knoten innerhalb des Polygons True sind (False, wenn nicht). Wenn Sie einen Test für viele Punkte durchführen, ist dies möglicherweise schneller ( obwohl dies davon abhängt, dass Sie einen Test innerhalb einer "Pixel" -Toleranz durchführen ):

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

Das Ergebnis ist:

Wenn Sie Geschwindigkeit benötigen und zusätzliche Abhängigkeiten kein Problem darstellen, finden Sie dies möglicherweise numbasehr nützlich (jetzt ist die Installation auf jeder Plattform recht einfach). Der von ray_tracingIhnen vorgeschlagene klassische Ansatz lässt sich leicht portieren, numbaindem Sie den numba @jitDekorator verwenden und das Polygon in ein numpy-Array umwandeln. Der Code sollte folgendermaßen aussehen:

@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

Die erste Ausführung dauert etwas länger als jeder nachfolgende Aufruf:

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

Was nach der Kompilierung auf Folgendes abfällt:

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

Wenn Sie beim ersten Aufruf der Funktion Geschwindigkeit benötigen, können Sie den Code in einem Modul mit vorkompilieren pycc. Speichern Sie die Funktion in einer src.py wie folgt:

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

Erstellen Sie es mit python src.pyund führen Sie es aus:

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

Im numba-Code habe ich verwendet: 'b1 (f8, f8, f8 [:,:])'

Zum Kompilieren nopython=Truemuss jede Variable vor dem deklariert werden for loop.

Im vorgefertigten src-Code lautet die Zeile:

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

Wird verwendet, um den Funktionsnamen und seine E / A-Vari-Typen, eine boolesche Ausgabe b1und zwei Gleitkommazahlen f8sowie ein zweidimensionales Array von Gleitkommazahlen f8[:,:]als Eingabe zu deklarieren .

Ihr Test ist gut, misst aber nur eine bestimmte Situation: Wir haben ein Polygon mit vielen Eckpunkten und einer langen Reihe von Punkten, um sie innerhalb des Polygons zu überprüfen.

Außerdem nehme ich an, dass Sie nicht die Matplotlib-Inside-Polygon-Methode gegen die Ray-Methode messen, sondern die Matplotlib-irgendwie-optimierte Iteration gegen die einfache Listen-Iteration

Lassen Sie uns N unabhängige Vergleiche durchführen (N Paare von Punkt und Polygon)?

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

Ergebnis:

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

Matplotlib ist immer noch viel besser, aber nicht 100-mal besser. Versuchen wir jetzt ein viel einfacheres Polygon ...

lenpoly = 5
# ... same code

Ergebnis:

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

Ich werde es einfach hier lassen, den obigen Code einfach mit numpy umschreiben, vielleicht findet es jemand nützlich:

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    

Wrapped ray_tracing in

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

Getestet an 100000 Punkten, Ergebnisse:

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