检查点是否包含在多边形/多多边形中-针对多个点

问题描述

我有一个2D地图，它被分成矩形网格(大约有45,000个)和一些从shapefile派生的多边形/多多边形(我目前使用shapefile库pyshp读取它们)。不幸的是，这些多边形中有几个相当复杂，由大量的点组成(其中一个点有64万个点)，可能会有洞。

我尝试做的是为每个多边形检查哪些单元中心(我的网格的单元)落在那个特定的多边形内。然而，拥有大约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多万次检查的运行时间很容易进入小时。

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

提前感谢您。

---

解决方案

如果您希望执行较少的比较操作，可以尝试使用形状字符串树功能。请考虑以下代码：

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