鉴于所有多边形都是一个分区，我如何快速查找包含一个点的多边形？

Question

我有一个函数

get_polygon(polygon_collection, point):
    for polygon in polygon_collection:
        if polygon.intersects(point):
            return polygon
    return None

这种方法有效，但它在 O(n) * O（单多边形检查）中。如果构建树数据结构，这肯定可以减少到 O(log(n)) * O(单多边形检查)。

shapely 是否直接支持这一点？

多边形集合的结构

• 该集合包含 n 个 MultiPolygons - 因此一个对象可以包含许多多边形。
• 每个 MultiPolygon 都可以有飞地/飞地
• 通常 (> 95%)，它们是没有孔的简单多边形
• 通常 (> 50%)，形状相当简单（例如正方形）

用例示例

多边形列表可以是德国的邮政编码区域。那将是数千。然后我有我和一些朋友的 GPS 位置，也有几千个。我想说的是我们在哪个区域获得了最多的数据点。

Answer 1

RTrees 正是您正在寻找的。 Shapely 现在支持 RTree 的构建和使用。但问题是匀称的 rtree 只支持矩形。因此，使用它的配方可以如下：

1. 构建包含多边形的最小矩形的 rtree。
2. 对于给定的点，使用 RTree 了解哪些矩形包含该点。该点需要稍微膨胀一点，因为 rtree 需要一个框来检查交叉点。这会给你一些误报，但没关系。
3. 过滤误报以获得最终结果。

除非您的多边形非常奇怪，否则这将为您提供平均值的对数时间（每个点的多边形数量的对数）。 当然，您需要为构建 RTree 结构付出代价，但如果多边形以某种方式均匀分布，那么这也是在对数时间内完成的，并且可以序列化 rtree 以供以后使用。 （也就是说，向树中添加一个新框与树中已经存在的框的数量成对数，因此总复杂度小于 n*log(n)。

实现如下。

from rtree import index
from shapely.geometry import Polygon, Point

def get_containing_box(p):
    xcoords = [x for x, _ in p.exterior.coords]
    ycoords = [y for _, y in p.exterior.coords]
    box = (min(xcoords), min(ycoords), max(xcoords), max(ycoords))   
    return box

def build_rtree(polys):
    def generate_items():
        pindx = 0
        for pol in polys:
            box = get_containing_box(pol)
            yield (pindx, box,  pol)
            pindx += 1
    return index.Index(generate_items())


def get_intersection_func(rtree_index):
    MIN_SIZE = 0.0001
    def intersection_func(point):
        # Inflate the point, since the RTree expects boxes:
        pbox = (point[0]-MIN_SIZE, point[1]-MIN_SIZE, 
                 point[0]+MIN_SIZE, point[1]+MIN_SIZE)
        hits = rtree_index.intersection(pbox, objects='raw')
        #Filter false positives:
        result = [pol for pol in hits if pol.intersects(Point(point)) ]
        return result
    return intersection_func

我用来测试的例子：

triangles = [
 [(8.774239095627754, 32.342041683826224),
  (8.750148703126552, 32.899346596303054),
  (4.919576457288677, 35.41040289384488)],
 [(8.485955148136222, 32.115258769399446),
  (9.263360720698277, 30.065319757618354),
  (4.562638192761559, 38.541192819415855)],
 [(2.613649959824923, 38.14802347408093),
  (7.879211442172622, 35.97223726358858),
  (0.9726266770834235, 32.12523430143625)]
]

polys = [Polygon(t) for t in triangles]
my_rtree = build_rtree(polys)
my_func = get_intersection_func(my_rtree)


my_intersections =  my_func((5, 35))
for pol in my_intersections:
    print pol.exterior.coords[:]