【算法分享】R树索引

在车联网、POI 分析、地理围栏、行政区划判断等业务中，我们经常会遇到一个非常基础、但极其高频的问题：

👉给定一个经纬度点，它属于哪个（或哪些）多边形？

比如：

• 一个 GPS 点属于哪个行政区？
• 一个车辆轨迹点是否落在某个物流园围栏内？
• 一个 POI 是否在规划的服务区域中？

这个问题在 GIS 领域被称为：Point in Polygon（PIP）问题。

本文将结合一份完整可运行的 Python 代码，从：

• 算法原理
• 工程实现
• 性能瓶颈
• 加速手段（包围盒、R 树）

几个层面，系统性地讲清楚：点在哪个多边形里，到底该怎么算、怎么选、怎么快。

一、最原始的方法：逐个多边形判断

1️⃣ 思路

最直观的方式是：

• 对每一个多边形
• 判断点是否在这个多边形内部

数学上常见的实现有：

• 射线法（Ray Casting）

  它计算从点P开始的射线穿过多边形边界的次数。当交叉数是偶数时，点在外面；当交叉数是奇数时，点在里面。

• 环绕数法（Winding Number）

  他计算多边形围着点P旋转地次数，只有当圈数wn=0时，点在外面；否则，点在里面。

在工程中，我们通常不会自己手写几何算法，而是直接使用成熟库。

2️⃣ 代码实现

defpoint_in_polygon(point_lon,point_lat,polygon_coords):""" 判断一个点是否在多边形内 参数: point_lon: 点的经度（float） point_lat: 点的纬度（float） polygon_coords: 多边形顶点列表，格式为 [(lon1, lat1), (lon2, lat2), ...] 返回: True: 点在多边形内 False: 点不在多边形内 """# 方式一：基于射线法importmatplotlib.pathasmplPathimportnumpyasnp polygon=mplPath.Path(np.array(polygon_coords))returnpolygon.contains_point(point)# 方式二：基于射线法importcv2 cv2.pointPolygonTest(np.array(polygon_coords),point,measureDist=False)# 方式三：基于环绕数法fromshapely.geometryimportPoint,Polygon point=Point(point_lon,point_lat)polygon=Polygon(polygon_coords)returnpolygon.contains(point)

如果有多个多边形：

deffind_point_in_polygons(point_lon,point_lat,polygons_dict):""" 查找点所在的多边形（不使用索引，适合小数据集） 参数: point_lon: 点的经度 point_lat: 点的纬度 polygons_dict: 多边形字典 返回: 包含该点的多边形名称列表 """point=Point(point_lon,point_lat)result=[]forname,coordsinpolygons_dict.items():polygon=Polygon(coords)ifpolygon.contains(point):result.append(name)returnresult

3️⃣ 性能分析

• 时间复杂度：O(N)
• N = 多边形数量

如果：

• 多边形数量 < 100
• 查询点数量不多

👉完全没问题，简单、清晰、好维护。

但一旦进入真实业务场景：

• 上千 / 上万围栏
• 每秒上百 / 上千点

性能会迅速成为瓶颈。

二、第一步加速：包围盒（Bounding Box）过滤

1️⃣ 核心思想

绝大多数点，
连多边形的“外接矩形”都不在里面。

那我们就可以：

1. 先判断点是否在多边形的包围盒内（极快）
2. 再对“可能命中”的少量多边形做精确判断

2️⃣ 什么是包围盒？

对一个多边形：

(minx, miny, maxx, maxy)

只要：

minx ≤ lon ≤ maxx miny ≤ lat ≤ maxy

才有可能在多边形内。

3️⃣ 代码结构

classPolygonIndexWithBounds:""" 使用包围盒快速过滤的空间索引 """def__init__(self,polygons_dict):"""初始化包围盒索引"""self.polygons={}self.bounds_map={}forname,coordsinpolygons_dict.items():polygon=Polygon(coords)self.polygons[name]=polygon self.bounds_map[name]=polygon.boundsdeffind_point_in_polygons(self,point_lon,point_lat):"""查找点所在的多边形"""point=Point(point_lon,point_lat)result=[]# 通过包围盒过滤candidates=[]forname,boundsinself.bounds_map.items():minx,miny,maxx,maxy=boundsifminx<=point_lon<=maxxandminy<=point_lat<=maxy:candidates.append(name)# 精确检查fornameincandidates:ifself.polygons[name].contains(point):result.append(name)returnresultdeffind_points_in_polygons_batch(self,points_list):"""批量查询"""results=[]forlon,latinpoints_list:polygons=self.find_point_in_polygons(lon,lat)results.append({"point":(lon,lat),"polygons":polygons})returnresults

4️⃣ 性能效果

• 包围盒判断：纯数值比较，极快
• 精确几何判断次数大幅减少

在测试中：

🚀通常可提速 10 ~ 100 倍

5️⃣ 适用场景

• 多边形数量：100 ~ 1000
• 实现成本低
• 不依赖额外空间索引库

👉性价比极高的工程方案

三、终极方案：R 树（Spatial Index）

1️⃣ 问题再升级

当你遇到：

• 上万 / 数十万多边形
• 实时或准实时查询

包围盒 O(N) 扫描也开始吃力。

这时候，就需要真正的空间索引。

2️⃣ 什么是 R 树？

R 树是一种：

为多维空间查询设计的树结构

步骤：

1. 对每个空间对象（点、线、面）构建最小外接矩形MBR。（如图D1|R8、D2|R9、D3|R10）
2. 建树过程：自底向上分组
   1. 把若干个对象的 MBR→ 组合成一个“节点”
   2. 对一个节点，再算一个更大的 MBR→ 覆盖所有子节点
   3. 节点数量超过阈值→ 进行 节点分裂（split）
3. 查询过程：从根开始逐层过滤，向下递归，到达叶子节点

特点：

• 多层包围盒，类似 B+ 树的分层结构
• 专为“空间相交 / 包含”优化，查询复杂度接近O(log N)

你可以把它理解为：

空间世界里的索引（像数据库索引一样）

3️⃣ Shapely 中的 R 树

Shapely 2.x 提供了STRtree：

fromshapely.strtreeimportSTRtree tree=STRtree(geometry_list)

查询时：

candidates=tree.query(point,predicate='intersects')

再进行一次精确判断即可。

classPolygonIndexWithRtree:""" 使用R树空间索引的多边形查询 性能：最快（100-1000x） Rtree使用动态边界体积树，专门为空间查询优化 """def__init__(self,polygons_dict):"""初始化R树索引"""try:fromshapely.strtreeimportSTRtreeexceptImportError:print("⚠️ 需要 shapely >= 2.0，使用 pip install --upgrade shapely")raiseself.polygons={}self.geometry_list=[]self.name_list=[]# 构建几何体列表用于R树forname,coordsinpolygons_dict.items():polygon=Polygon(coords)self.polygons[name]=polygon self.geometry_list.append(polygon)self.name_list.append(name)# 创建R树索引self.tree=STRtree(self.geometry_list)deffind_point_in_polygons(self,point_lon,point_lat):""" 使用R树查询点所在的多边形 """point=Point(point_lon,point_lat)result=[]# R树查询：获取可能包含该点的多边形索引candidates_idx=self.tree.query(point,predicate='intersects')# 精确检查（因为intersects只是粗略检查）foridxincandidates_idx:ifself.geometry_list[idx].contains(point):result.append(self.name_list[idx])returnresultdeffind_points_in_polygons_batch(self,points_list):"""批量查询"""results=[]forlon,latinpoints_list:polygons=self.find_point_in_polygons(lon,lat)results.append({"point":(lon,lat),"polygons":polygons})returnresults

4️⃣ 性能表现

在示例代码中，使用：

• 1000 个多边形
• 500 个查询点

代码链接：点击查看

测试结果：

多边形越多，R 树优势越明显。

四、三种方案的对比总结

五、工程选型建议

• 多边形 < 100 个： 直接 contains 判断
• 多边形 100 ~ 1000： 包围盒过滤
• 多边形 > 1000： R 树索引（强烈推荐）

📬 关注我 · 获取更多内容

📌 南墨的技术小栈

这里是我的个人知识分享空间。我会定期整理和分享工作与学习中积累的经验与资源，内容涵盖：

• 算法分享 —— 深入讲解算法原理、实现思路与代码示例。
• 工具分享 —— 推荐实用工具与脚本，包括我个人开发的小工具和精选开源工具。
• 开源项目 —— 精选 GitHub 上高星项目，拆解原理、使用方法和最佳实践。
• 数据分享 —— 工作学习中收集整理的数据资源。

无论你是技术爱好者、算法研究者，还是对数据与开源感兴趣的朋友，这里都希望能成为你学习、探索和实践的参考空间。

若在阅读或使用过程中有任何疑问，欢迎在公众号私信我，我会尽快与您交流。