programing

파이썬에서 포인트가 폴리곤 내부에 있는지 확인하는 가장 빠른 방법은 무엇입니까?

randomtip 2023. 7. 20. 23:11
반응형

파이썬에서 포인트가 폴리곤 내부에 있는지 확인하는 가장 빠른 방법은 무엇입니까?

점이 폴리곤 내부에 속하는지 확인하는 두 가지 주요 방법을 찾았습니다.하나는 여기서 사용되는 광선 추적 방법을 사용하는 것입니다. 이것은 가장 권장되는 답입니다. 다른 하나는 matplotlib을 사용하는 것입니다.path.contains_points(내가 보기엔 좀 애매해 보입니다.)저는 계속해서 많은 점을 확인해야 할 것입니다.이 두 가지 중에 다른 것보다 더 추천할 만한 것이 있는지 또는 더 나은 세 번째 옵션이 있는지 아는 사람이 있습니까?

업데이트:

두 가지 방법을 확인해보니 matplotlib이 훨씬 빨라 보입니다.

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

그러면,

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

100변 다각형 대신 삼각형을 사용하여 동일한 상대적 차이를 얻었습니다.이런 문제들을 전담하는 패키지로 보여 자세히 확인해 보겠습니다.

구체적으로 고려할 수 있습니다.

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

당신이 말한 방법 중에서 두 번째 방법을 사용했을 뿐입니다.path.contains_points그리고 그것은 잘 작동합니다.어떤 경우든 테스트에 필요한 정밀도에 따라 다각형 내부의 모든 노드가 참(그렇지 않으면 거짓)이 되도록 numpy bool 그리드를 생성하는 것이 좋습니다.많은 점에 대해 테스트를 수행하려는 경우 이 작업이 더 빠를있습니다("픽셀" 공차" 내에서 테스트를 수행해야 합니다).

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

결과는 다음과 같습니다.

point inside polygon within pixel tolerance

속도가 필요하고 추가적인 의존성이 문제가 되지 않는다면 매우 유용할 것입니다(이제 모든 플랫폼에 설치하기가 매우 쉽습니다).클래식ray_tracing 제한접방쉽이수있습다니할식으로 이동할 수.numba을 이용하여numba @jit장식가 및 다각형을 numpy 배열에 캐스팅합니다.코드는 다음과 같아야 합니다.

@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

첫 번째 실행은 다음 호출보다 약간 더 오래 걸립니다.

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

컴파일 후에는 다음과 같이 감소합니다.

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

의 첫 컴파일하여 할 수 있습니다.pycc다음과 같은 기능을 src.py 에 저장합니다.

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

으로 구축:python src.py실행:

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

내가 사용한 numba 코드에서: 'b1(f8, f8, f8[:,:])'

로컴하위해기로 nopython=True는 각변는다이선합니다야언어 앞에 .for loop.

프리빌드 src 코드에서 다음 행을 지정합니다.

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

I/O할 때 됩니다. 출력 " " " I/O " " ", " " " " ", " " " 입니다.b1그리고 두 개의 수레f8그리고 2차원의 플로트 배열.f8[:,:]입력으로서

2021년 1월 4일 편집

사용 사례의 경우, 여러 점이 단일 폴리곤 내에 있는지 확인해야 합니다. 이러한 상황에서 numba 병렬 기능을 활용하여 일련의 점을 반복하는 것이 유용합니다.위의 예는 다음과 같이 변경할 수 있습니다.

from numba import jit, njit
import numba
import numpy as np 

@jit(nopython=True)
def pointinpolygon(x,y,poly):
    n = len(poly)
    inside = False
    p2x = 0.0
    p2y = 0.0
    xints = 0.0
    p1x,p1y = poly[0]
    for i in numba.prange(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


@njit(parallel=True)
def parallelpointinpolygon(points, polygon):
    D = np.empty(len(points), dtype=numba.boolean) 
    for i in numba.prange(0, len(D)):
        D[i] = pointinpolygon(points[i,0], points[i,1], polygon)
    return D    

참고: 위의 코드를 미리 컴파일하면 numba의 병렬 기능이 활성화되지 않습니다(병렬 CPU 대상은 에서 지원되지 않음).pycc/AOT컴파일) 참조: https://github.com/numba/numba/issues/3336

테스트:


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]]
polygon = np.array(polygon)
N = 10000
points = np.random.uniform(-1.5, 1.5, size=(N, 2))

위해서N=1000072 코어 시스템에서 다음을 반환합니다.

%%timeit
parallelpointinpolygon(points, polygon)
# 480 µs ± 8.19 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

21년 2월 17일 편집:

  • 를 작할루프에서 하는 것01(고https워 @마di):

for i in numba.prange(0, len(D))

21년 2월 20일 편집:

@mehdi의 비교에 대한 후속 조치로, 아래에 GPU 기반 방법을 추가합니다.사용하는 방법은 다음과 같습니다.cuspatial라이브러리:

    import numpy as np
    import cudf
    import cuspatial

    N = 100000002
    lenpoly = 1000
    polygon = [[np.sin(x)+0.5,np.cos(x)+0.5] for x in 
    np.linspace(0,2*np.pi,lenpoly)]
    polygon = np.array(polygon)
    points = np.random.uniform(-1.5, 1.5, size=(N, 2))


    x_pnt = points[:,0]
    y_pnt = points[:,1]
    x_poly =polygon[:,0]
    y_poly = polygon[:,1]
    result = cuspatial.point_in_polygon(
        x_pnt,
        y_pnt,
        cudf.Series([0], index=['geom']),
        cudf.Series([0], name='r_pos', dtype='int32'), 
        x_poly, 
        y_poly,
    )

@Mehdi ▁@.N=100000002그리고.lenpoly=1000다음과 같은 결과를 얻었습니다.

 time_parallelpointinpolygon:         161.54760098457336 
 time_mpltPath:                       307.1664695739746 
 time_ray_tracing_numpy_numba:        353.07356882095337 
 time_is_inside_sm_parallel:          37.45389246940613 
 time_is_inside_postgis_parallel:     127.13793849945068 
 time_is_inside_rapids:               4.246025562286377

point in poligon methods comparison, #poins: 10e07

하드웨어 사양:

  • CPU Intel E1240
  • GPU Nvidia GTX 1070

주의:

  • cuspatial.point_in_poligon메소드는 매우 강력하고 강력하며, 여러 개의 복잡한 다각형으로 작업할 수 있는 기능을 제공합니다(성능을 희생해야 할 것 같습니다).

  • numba방법은 GPU에서도 '포팅'할 수 있습니다. 포팅을 포함한 비교를 보는 것은 흥미로울 것입니다.cuda@Mehdi가 언급한 가장 빠른 방법 (is_inside_sm).

검정은 양호하지만 일부 특정 상황만 측정합니다. 즉, 정점이 많은 하나의 다각형과 다각형 내에서 이를 확인할 수 있는 긴 점 배열이 있습니다.

게다가, 저는 당신이 matplot lib-inside-polygon-method 대 ray-method가 아니라 matplot lib-some-optimized-reitation 대 simple-list-reitation을 측정하고 있다고 생각합니다.

N개의 독립적인 비교(N개의 점과 다각형 쌍)를 만들어 봅시다.

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

결과:

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

Matplotlib은 여전히 훨씬 낫지만 100배는 아닙니다.이제 훨씬 더 단순한 다각형을 시도해 보겠습니다...

lenpoly = 5
# ... same code

결과:

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

다양한 방법의 비교

점이 폴리곤 내부에 있는지 확인하는 다른 방법을 찾았습니다(여기).저는 그 중 두 가지 방법(is_inside_sm과 is_inside_postgis)만 테스트했고 결과는 다른 방법과 동일합니다.

@epifanio 덕분에 코드를 병렬화하여 @epifanio 및 @user3274748(ray_tracing_numpy) 메서드와 비교했습니다.참고로 두 방법 모두 버그가 있어서 아래 코드와 같이 수정하였습니다.

가지 더 것입니다.np.linspace(0,2*np.pi,lenpoly)[:-1]따라서 위의 GitHub 저장소에 제공된 코드가 제대로 작동하지 않을 수 있습니다.따라서 닫힌 경로를 만드는 것이 좋습니다(첫 번째와 마지막 점은 동일해야 함).

코드

방법 1: 다각형의 평행점

from numba import jit, njit
import numba
import numpy as np 

@jit(nopython=True)
def pointinpolygon(x,y,poly):
    n = len(poly)
    inside = False
    p2x = 0.0
    p2y = 0.0
    xints = 0.0
    p1x,p1y = poly[0]
    for i in numba.prange(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


@njit(parallel=True)
def parallelpointinpolygon(points, polygon):
    D = np.empty(len(points), dtype=numba.boolean) 
    for i in numba.prange(0, len(D)):   #<-- Fixed here, must start from zero
        D[i] = pointinpolygon(points[i,0], points[i,1], polygon)
    return D  

방법 2: ray_tracing_numpy_numba

@jit(nopython=True)
def ray_tracing_numpy_numba(points,poly):
    x,y = points[:,0], points[:,1]
    n = len(poly)
    inside = np.zeros(len(x),np.bool_)
    p2x = 0.0
    p2y = 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 len(idx):    # <-- Fixed here. If idx is null skip comparisons below.
            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 

방법 3: Matplotlib에 _points가 포함되어 있습니다.

path = mpltPath.Path(polygon,closed=True)  # <-- Very important to mention that the path 
                                           #     is closed (default is false)

방법 4: is_inside_sm(여기서 가져왔습니다)

@jit(nopython=True)
def is_inside_sm(polygon, point):
    length = len(polygon)-1
    dy2 = point[1] - polygon[0][1]
    intersections = 0
    ii = 0
    jj = 1

    while ii<length:
        dy  = dy2
        dy2 = point[1] - polygon[jj][1]

        # consider only lines which are not completely above/bellow/right from the point
        if dy*dy2 <= 0.0 and (point[0] >= polygon[ii][0] or point[0] >= polygon[jj][0]):

            # non-horizontal line
            if dy<0 or dy2<0:
                F = dy*(polygon[jj][0] - polygon[ii][0])/(dy-dy2) + polygon[ii][0]

                if point[0] > F: # if line is left from the point - the ray moving towards left, will intersect it
                    intersections += 1
                elif point[0] == F: # point on line
                    return 2

            # point on upper peak (dy2=dx2=0) or horizontal line (dy=dy2=0 and dx*dx2<=0)
            elif dy2==0 and (point[0]==polygon[jj][0] or (dy==0 and (point[0]-polygon[ii][0])*(point[0]-polygon[jj][0])<=0)):
                return 2

        ii = jj
        jj += 1

    #print 'intersections =', intersections
    return intersections & 1  


@njit(parallel=True)
def is_inside_sm_parallel(points, polygon):
    ln = len(points)
    D = np.empty(ln, dtype=numba.boolean) 
    for i in numba.prange(ln):
        D[i] = is_inside_sm(polygon,points[i])
    return D  

방법 5: is_inside_postgis(여기서 가져옵니다)

@jit(nopython=True)
def is_inside_postgis(polygon, point):
    length = len(polygon)
    intersections = 0

    dx2 = point[0] - polygon[0][0]
    dy2 = point[1] - polygon[0][1]
    ii = 0
    jj = 1

    while jj<length:
        dx  = dx2
        dy  = dy2
        dx2 = point[0] - polygon[jj][0]
        dy2 = point[1] - polygon[jj][1]

        F =(dx-dx2)*dy - dx*(dy-dy2);
        if 0.0==F and dx*dx2<=0 and dy*dy2<=0:
            return 2;

        if (dy>=0 and dy2<0) or (dy2>=0 and dy<0):
            if F > 0:
                intersections += 1
            elif F < 0:
                intersections -= 1

        ii = jj
        jj += 1

    #print 'intersections =', intersections
    return intersections != 0  


@njit(parallel=True)
def is_inside_postgis_parallel(points, polygon):
    ln = len(points)
    D = np.empty(ln, dtype=numba.boolean) 
    for i in numba.prange(ln):
        D[i] = is_inside_postgis(polygon,points[i])
    return D  

벤치마크

enter image description here

1,000만 포인트에 대한 타이밍:

parallelpointinpolygon Elapsed time:      4.0122294425964355
Matplotlib contains_points Elapsed time: 14.117807388305664
ray_tracing_numpy_numba Elapsed time:     7.908452272415161
sm_parallel Elapsed time:                 0.7710440158843994
is_inside_postgis_parallel Elapsed time:  2.131121873855591

여기 코드가 있습니다.

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

np.random.seed(2)

time_parallelpointinpolygon=[]
time_mpltPath=[]
time_ray_tracing_numpy_numba=[]
time_is_inside_sm_parallel=[]
time_is_inside_postgis_parallel=[]
n_points=[]

for i in range(1, 10000002, 1000000): 
    n_points.append(i)
    
    lenpoly = 100
    polygon = [[np.sin(x)+0.5,np.cos(x)+0.5] for x in np.linspace(0,2*np.pi,lenpoly)]
    polygon = np.array(polygon)
    N = i
    points = np.random.uniform(-1.5, 1.5, size=(N, 2))
    
    
    #Method 1
    start_time = time()
    inside1=parallelpointinpolygon(points, polygon)
    time_parallelpointinpolygon.append(time()-start_time)

    # Method 2
    start_time = time()
    path = mpltPath.Path(polygon,closed=True)
    inside2 = path.contains_points(points)
    time_mpltPath.append(time()-start_time)

    # Method 3
    start_time = time()
    inside3=ray_tracing_numpy_numba(points,polygon)
    time_ray_tracing_numpy_numba.append(time()-start_time)

    # Method 4
    start_time = time()
    inside4=is_inside_sm_parallel(points,polygon)
    time_is_inside_sm_parallel.append(time()-start_time)

    # Method 5
    start_time = time()
    inside5=is_inside_postgis_parallel(points,polygon)
    time_is_inside_postgis_parallel.append(time()-start_time)


    
plt.plot(n_points,time_parallelpointinpolygon,label='parallelpointinpolygon')
plt.plot(n_points,time_mpltPath,label='mpltPath')
plt.plot(n_points,time_ray_tracing_numpy_numba,label='ray_tracing_numpy_numba')
plt.plot(n_points,time_is_inside_sm_parallel,label='is_inside_sm_parallel')
plt.plot(n_points,time_is_inside_postgis_parallel,label='is_inside_postgis_parallel')
plt.xlabel("N points")
plt.ylabel("time (sec)")
plt.legend(loc = 'best')
plt.show()

결론

가장 빠른 알고리즘은 다음과 같습니다.

1- is_sm_sm_sm_sm

2- is_is_nots_postgis_nots

폴리곤의 3-평행점(@epifanio)

그냥 여기에 두고, 위의 코드를 numpy를 사용하여 다시 작성하면 누군가 유용하게 사용할 수 있습니다.

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 대상으로 포장

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

100000점에서 테스트한 결과:

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

짝수 홀수 규칙의 순수한 무감각한 벡터화된 구현

다른 답변은 느린 파이썬 루프이거나 외부 의존성 또는 사이톤 처리가 필요합니다.

import numpy as np
        
def points_in_polygon(polygon, pts):
    pts = np.asarray(pts,dtype='float32')
    polygon = np.asarray(polygon,dtype='float32')
    contour2 = np.vstack((polygon[1:], polygon[:1]))
    test_diff = contour2-polygon
    mask1 = (pts[:,None] == polygon).all(-1).any(-1)
    m1 = (polygon[:,1] > pts[:,None,1]) != (contour2[:,1] > pts[:,None,1])
    slope = ((pts[:,None,0]-polygon[:,0])*test_diff[:,1])-(test_diff[:,0]*(pts[:,None,1]-polygon[:,1]))
    m2 = slope == 0
    mask2 = (m1 & m2).any(-1)
    m3 = (slope < 0) != (contour2[:,1] < polygon[:,1])
    m4 = m1 & m3
    count = np.count_nonzero(m4,axis=-1)
    mask3 = ~(count%2==0)
    mask = mask1 | mask2 | mask3
    return mask

    
N = 1000000
lenpoly = 1000
polygon = [[np.sin(x)+0.5,np.cos(x)+0.5] for x in np.linspace(0,2*np.pi,lenpoly)]
polygon = np.array(polygon,dtype='float32')
points = np.random.uniform(-1.5, 1.5, size=(N, 2)).astype('float32')
mask = points_in_polygon(polygon, points)

크기가 1000인 폴리곤이 있는 1mill 포인트는 44s를 찍었습니다.

다른 구현보다 속도는 느리지만 파이썬 루프보다 빠르며 numpy만 사용합니다.

inpolyPython에서 폴리곤 검사를 수행하기 위한 골드 표준이며, 다음과 같이 방대한 쿼리를 처리할 수 있습니다.

https://github.com/dengwirda/inpoly-python

단순 사용:

from inpoly import inpoly2
import numpy as np
    
xmin, xmax, ymin, ymax = 0, 1, 0, 1
x0, y0, x1, y1 = 0.5, 0.5, 0, 1

#define any n-sided polygon
p = np.array([[xmin, ymin],
              [xmax, ymin],
              [xmax, ymax],
              [xmin, ymax],
              [xmin, ymin]])

#define some coords
coords = np.array([[x0, y0],
                   [x1, y1]])

#get boolean mask for points if in or on polygon perimeter
isin, ison = inpoly2(coords, p)

백엔드의 C 구현은 매우 빠릅니다.

지오판다를 조인으로 사용할 수 있습니다.아마도 가장 쉬운.

import geopandas as gpd
gpd.sjoin(points, polygon, op = 'within') 

언급URL : https://stackoverflow.com/questions/36399381/whats-the-fastest-way-of-checking-if-a-point-is-inside-a-polygon-in-python

반응형