python - 检查点是否包含在多边形/多边形中 - 对于许多点

Question

我有一个 2D map ，分为矩形网格 - 其中大约 45,000 个 - 以及一些从 shapefile 派生的多边形/多多边形(我目前使用 shapefile 库 pyshp 读取它们)。不幸的是，其中一些多边形相当复杂，由大量点组成(其中一个有 640,000 个点)，并且其中可能有孔。

我要做的是检查每个多边形的哪些单元格中心(我的网格的单元格)位于该特定多边形内。然而，拥有大约 45,000 个单元中心和 150 个多边形需要相当长的时间才能使用 shapely 检查所有内容。这就是我正在做的，或多或少:

# nx and ny are the number of cells in the x and y direction respectively
# x, y are the coordinates of the cell centers in the grid - numpy arrays
# shapes is a list of shapely shapes - either Polygon or MultiPolygon
# Point is a shapely.geometry.Point class
for j in range(ny):
    for i in range(nx):
        p = Point([x[i], y[j]])
        for s in shapes:
            if s.contains(p):
                # The shape s contains the point p - remove the cell

根据所讨论的多边形，每次检查需要 200 微秒到 80 毫秒，并且有超过 600 万次检查，运行时间很容易进入数小时。

我想一定有更聪明的方法来处理这个问题 - 如果可能的话，我宁愿保持身材，但任何建议都非常感谢。

提前谢谢你。

Answer 1

在接受的答案中，在 shapely 中使用 rtrees 仅支持多边形的扩展或边界框。正如 shapely 文档中所述:

https://shapely.readthedocs.io/en/stable/manual.html#str-packed-r-tree

从 1.7.1 开始，rtree.query 之后需要检查多边形是否实际包含该点。这意味着额外的检查，但如果多边形不是非常困惑的配置，仍然可以加快算法

也就是说，代码应该是这样的

from shapely.strtree import STRtree
from shapely.geometry import Polygon, Point    

# this is the grid. It includes a point for each rectangle center. 
points = [Point(i, j) for i in range(10) for j in range(10)]
tree = STRtree(points) # create a 'database' of points

# this is one of the polygons. 
p = Polygon([[0.5, 0], [2, 0.2], [3, 4], [0.8, 0.5], [0.5, 0]])

res = [o for o in tree.query(p) if p.contains(o)] # run the query (a single shot) - and test if the found points are actually inside the polygon. 

# do something with the results
print([o.wkt for o in res])

果然，结果就是现在

['POINT (2 1)', 'POINT (2 2)']

这里可以看出是多边形内部的点:

-- 奖金编辑 ---

我观察到通过在以下意义上定位问题可以更好地提高速度:

在该区域上制作了一个方形网格，STRtree 结构由点和多边形组成。在每个单独的网格多边形上查询点和多边形 - 在每个网格多边形内，检查查询的点是否包含在查询的多边形内。从而减少对单个多边形的检查量。

from shapely.strtree import STRtree
from shapely.geometry import Polygon, Point
import random    

# build a grid spread over the area. 
grid = [Polygon([[i, j],[i+1,j],[i,j+1],[i+1,j+1]]) for i in range(10) for j in range(10)]

pointtree = STRtree([Point(random.uniform(1,10),random.uniform(1,10)) for q in range(50)]) # create a 'database' of 50 random points - and build a STRtree structure over it

# take the same polygon as above, but put in an STRtree 
polytree = STRtree([Polygon([[0.5, 0], [2, 0.2], [3, 4], [0.8, 0.5], [0.5, 0]])])
res = []
#do nested queries inside the grid
for g in grid:
    polyquery = polytree.query(g)
    pointquery = pointtree.query(g)
    for point in pointquery:
        for poly in polyquery:
            if poly.contains(point):
                 res.append(point)

# do something with the results
print([o.wkt for o in res])