检查一个点是否在多边形/多部分多边形中-对于多个点

Question

我有一张2D的地图，被分成一个矩形网格 - 大约有45,000个矩形 - 和一些从shapefiles导出的多边形/多部分多边形（我目前使用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，但任何建议都非常感谢。

提前感谢您。

- Infinity77

---

也许你可以进行一些初步的排除？我建议先确定每个形状的边界矩形（这只需要进行一次），然后进行快速检查（4个操作）以验证一个点是否可能属于边界矩形。如果它不可能属于边界矩形，那么就没有必要验证它是否可能属于底层形状。 - Joris Schellekens

Answer 1

在被接受的答案中，Shapely 中使用 R 树仅支持多边形的扩展或边界框。正如 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])

- Tarje Bargheer

Answer 2

如果您想执行更少的比较操作，可以尝试使用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)']

- Roy2012

---

谢谢，这似乎很有效并且速度非常快。我的数据或STRtree中可能有一些问题，因为对于嵌入在shapefile中的一个形状，STRtree告诉我该形状内有43,621个点，而总单元格数为43,634，这意味着该形状包含几乎所有点。但是，如果我可视化该形状，情况绝对不是这样... - Infinity77

很有趣。我对Shapey不够熟悉，无法就其发表任何有见地的评论。一个想法：如果你在形状内得到了太多点，请使用“慢速”方法进行验证。总体而言，它仍然会运行得更快。 - Roy2012

顺便提一下，如果这个答案有价值的话，您是否介意接受它，以便未来的人可以参考？ - Roy2012