检查一个点是否包含在多边形多角形中--对于许多点来说

Question

我有一张用矩形网格划分的二维地图--大约有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万次的检查，运行时间很容易就到了几小时。

我想一定有更聪明的方法来处理这个问题--如果可能的话，我宁愿继续使用shapely，但任何建议都是非常感激的。

先谢谢你。

Answer 1

如果你想减少比较操作，你可以尝试使用shapely str-tree功能。考虑以下代码。

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 = tree.query(p). # run the query (a single shot). 

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

这段代码的结果是：

['POINT (3 0)', 'POINT (2 0)', 'POINT (1 0)', 'POINT (1 1)', 'POINT (3 1)',
 'POINT (2 1)', 'POINT (2 2)', 'POINT (3 2)', 'POINT (1 2)', 'POINT (2 3)',
 'POINT (3 3)', 'POINT (1 3)', 'POINT (2 4)', 'POINT (1 4)', 'POINT (3 4)']