用R树加速空间检索

几年前，我正在做一个业余项目。 我想创建一个网络应用程序来推荐当地的瑰宝，例如咖啡馆、书店或隐藏的酒吧。 这个想法是在地图上显示用户触手可及的所有此类兴趣点。 由于数据集中有数十万个点，我必须巧妙地过滤用户给定范围内的数据点。 简单的方法是计算用户和每个兴趣点之间的距离，并丢弃指定范围之外的所有点。 特别是对于像我这样的大数据集，这种方法通常会导致处理时间很长。

当然，必须有更好的方法，因为响应时间在交互式应用程序中很重要。 这是我接触到数据结构 R 树的时候。 这些树用于快速空间访问和搜索。 使用 R 树，我能够快速隔离靠近用户位置的兴趣点并将其显示在地图上。 这使我的网络应用程序的响应时间得到了巨大的提升——只需要额外的四行代码！

在本文中，我将解释什么是 R 树以及它们如何工作。 前两节中以纽约市行道树为例对此进行了说明。 第三部分演示了如何在 Python 中使用此数据结构来加速地理空间数据处理例程。

1、通过分析纽约市的树木来学习 R 树

假设我们被要求分析纽约市的社区与其树木的健康状况是否存在相关性。 纽约市开放数据门户提供街道树木普查数据集，其中包括每棵树的物种、直径、健康感知和地理位置。

首先，我们要计算上东区行道树的数量。 下面的伪代码片段迭代数据集树并检查树是否落在 upper_east_side 边界内：

total_count, tree_count = 0, 0
for tree in trees:
    total_count += 1
    if upper_east_side.contains(tree):
        tree_count += 1
print_results(num_tests=total_count, num_trees=tree_count)
>>> Total number of trees tested: 683,788
>>> Number of trees in Upper East Side: 8,807

我们发现上东区大约有 9,000 棵树。 然而，我们必须测试总共 684k 棵树才能达到这个目标。 下面的动画形象地展示了我们测试的树木距离我们的目标邻居数英里，因此很容易被忽视。 但是我们如何才能将远处的树排除在昂贵的计算之外以实现显着的性能提升呢？

我们实际上免费获得的一项信息是多边形的边界框（它可以由其节点的最小值和最大值确定）。 此外，测试一个点是否落在矩形内很简单，只需要四次比较操作（该点必须大于或等于左下角，并且小于或等于右上角）。 现在，假设 bounding_box是一个数据集，其中包含上东区周围紧密矩形中的所有树木（在下一节中我们将学习如何轻松获得这样的矩形）。 考虑到这一点，得出：

total_count, tree_count = 0, 0
for tree in bounding_box:
    total_count += 1
    if upper_east_side.contains(tree):
        tree_count += 1
print_results(num_tests=total_count, num_trees=tree_count)
>>> Total number of trees tested: 10,768
>>> Number of trees in Upper East Side: 8,807

动画的右侧表明我们现在只测试潜在的候选者。 这些是紧邻多边形的树，即落在其边界框内的点。 通过忽略远处的树木，我们能够将测试数量从 684k 减少到 11k — 60 倍！ 在下一节中，我们将看到 R 树正是利用了这个想法。

2、用于空间搜索的数据结构：R 树

R 树是基于树的数据结构，用于以有效的方式创建空间索引。 R 树通常用于快速空间查询或加速最近邻搜索 [1]。 一个常见的用例可能是存储兴趣点的空间信息（例如餐馆、加油站、街道等）。在 R 树的帮助下，人们可以快速检索到某个位置一定距离内的所有兴趣点。 作为回报，这些结果可以显示在地图或导航系统上。

R 树的基本思想很简单：树的叶节点保存空间数据，而分支节点对应于包含其所有子节点的最小边界框。 通过这种结构，R 树将空间划分为矩形，随着树的生长，矩形变得更加细化。 下面的示例对此进行了说明。

在 R 树中查询矩形，即我们想要检索此搜索窗口中包含的所有数据。 请记住，每个非叶节点对应于一个包含其所有子节点的边界框。 为了完成搜索查询，我们只需沿着树的分支行进，并沿着与给定矩形相交的路径，直到到达叶节点。 这些叶节点以及我们的数据点都包含在搜索矩形中并满足查询。 下面的动画演示了我们可以通过忽略不符合搜索条件的整个分支来大大减少搜索操作的数量。

3、Python 中的 R 树

Python 包 Rtree 提供了 R 树数据结构的实现，并附带许多方便的功能，例如最近邻搜索、交集搜索或多维索引。

我们可以使用Python的包管理器pip方便地安装包： pip install Rtree。

3.1 基本用法

在处理点或多边形等几何图形之前，我们先介绍 Rtree 包的基本用法。

from rtree import index
idx = index.Index()
bbox_0 = (0.0, 0.0, 1.0, 1.0) # (left, bottom, right, top)
bbox_1 = (3.0, 3.0, 6.0, 6.0)
idx.insert(0, bbox_0)
idx.insert(1, bbox_1)

index模块帮助我们构建空间索引。 该索引是通过插入对象的边界框自动建立的。 边界框是通过指定其左、下、右和上坐标来定义的。 请注意，我们将边界框与标识符一起插入（在上面的示例中为 0 和 1）。 ID 将帮助我们在执行查询时识别边界框：

search_window = (-1.0, -1.0, 2.0, 2.0) # ex. 1
list(idx.intersection(search_window))
>>> [0]
search_window = (4.0, 4.0, 5.0, 5.0)   # ex. 2
list(idx.intersection(search_window))
>>> [1]
search_window = (0.0, 0.0, 6.0, 6.0)   # ex. 3
list(idx.intersection(search_window))
>>> [0, 1]
search_window = (1.01, 1.01, 2.0, 2.0) # ex. 4
list(idx.intersection(search_window))
>>> []

查询给定矩形的索引，同样由其左、下、右和上坐标指定。 交集方法的结果是搜索窗口中包含的对象的 ID（示例 1-3）。 如果搜索窗口超出索引中的数据范围，则结果为空（示例 4）。 类似地，我们使用最近的方法来查找距离给定搜索窗口最近的 k 个对象：

k = 1
search_window = (1.01, 1.01, 2.0, 2.0)
list(idx.nearest(search_window, k))
>>> [0]

3.2 处理点、线和多边形

在上一节中，我们了解了如何通过插入对象的边界框来构造索引。 我们现在想要继续对这些对象使用点、线和多边形。 Shapely 包提供了一种在 Python 中处理此类几何图形的简单方法：

from shapely.geometry import Point, LineString, Polygon
point = Point(0, 0)
line = LineString([(1, 0), (2, 1)])
polygon = Polygon([(3, 0), (3, 1), (4, 1)])
idx = index.Index()
idx.insert(0, point.bounds)   # bounds: (0.0, 0.0, 0.0, 0.0)
idx.insert(1, line.bounds)    # bounds: (1.0, 0.0, 2.0, 1.0)
idx.insert(2, polygon.bounds) # bounds: (3.0, 0.0, 4.0, 1.0)

上面，我们首先创建了一个点、一条线和一个多边形。 接下来，这些对象的边界框被插入到使用 ID 0、1 和 2 的索引中。我们现在查询不同搜索窗口的索引：

search_window = (0.0, 0.0, 0.5, 0.5)  # ex. 1
list(idx.intersection(search_window))
>>> [0]
search_window = (0.0, 0.0, 2.0, 1.0)  # ex. 2
list(idx.intersection(search_window))
>>> [0, 1]
search_window = (0.0, 0.0, 4.0, 2.0)  # ex. 3
list(idx.intersection(search_window))
>>> [0, 1, 2]
k = 1
search_window = (1.0, 0.0, 2.5, 1.5)  # ex. 4
list(idx.nearest(search_window, k))
>>> [1]

下图显示了几何图形和搜索窗口：

3.3 搜索上东区的所有树木

我们终于拥有了提取上东区所有树木所需的一切！ 我们将浏览下面的代码片段，但可以在此处找到完整版本。

首先，我们使用 GeoPandas 包加载所有必需的几何图形：

import geopandas as gpd
from rtree import index
df_trees = gpd.read_file('2015 Street Tree Census.geojson')
df_upper_east_side = gpd.read_file('upper_east_side.geojson')
geom_upper_east_side = df_upper_east_side.loc[0, 'geometry']

接下来，我们创建一个包含纽约市所有树木的 R 树索引：

idx = index.Index()
for id, geom in zip(df_trees.index, df_trees.geometry):
    idx.insert(id, geom.bounds)

现在，我们生成潜在候选者的列表，即上东区边界框中的所有树木：

search_window = geom_upper_east_side.bounds
potential = list(idx.intersection(search_window))

最后，我们遍历所有潜在的候选者，以提取完全位于上东区的候选者：

trees = []
for tree in df_trees.loc[potential, 'geometry']:
    if geom_upper_east_side.contains(tree):
        trees.append(tree)

print_results(num_data=len(df_trees), num_tests=len(potential), num_trees=len(trees))
>>> Total number of trees in dataset: 683,788
>>> Total number of trees tested: 10,768
>>> Number of trees in Upper East Side: 8,807

4、结束语

在本文中，我们了解了 R 树如何通过将底层空间划分为矩形来组织地理信息。 这种结构使得 R 树的空间查找速度非常快。 在我们的纽约市街道树示例中，利用 R 树将操作数量减少了 60 倍。我们还了解了如何在 Python 中使用 R 树。 我们的示例中的加速仅用四行代码即可实现：初始化索引（1 行）、构建索引（2 行）以及使用交集函数查找附近的候选者（1 行）。

那么为什么 R 树没有到处使用呢？ 虽然我们通过减少搜索操作次数获得了时间，但我们通过构建索引而损失了时间。 对于后者，我们实际上必须迭代整个数据集。 这使得 R 树不适合仅需要少量搜索的应用程序或索引经常更改的应用程序（由于树重新平衡）。

自 Antonin Guttman 于 1984 年发明以来，R 树已经取得了长足的进步。如今，它们出现在各种应用中，例如计算机图形学 [2]、视频游戏 [3]、交通控制系统 [4]、 最突出的是空间数据管理数据库[5]。 也许在你的下一次地理空间数据分析中也是如此！

原文链接：Speed up your Geospatial Data Analysis with R-Trees

BimAnt翻译整理，转载请标明出处