在处理大规模地理空间数据时，空间查询的性能是一个关键问题。空间索引通过对几何对象的空间位置建立索引结构，可以显著加速空间查询操作。本章将全面介绍 GeoPandas 中空间索引的原理、使用方法和优化策略。
在没有空间索引的情况下，执行空间查询（如"找出与某区域相交的所有要素"）需要逐一检查每个要素，时间复杂度为 O(n)。当数据量达到数十万甚至数百万时，这种暴力搜索将变得非常缓慢。
空间索引通过预先组织数据的空间分布信息，将查询的平均时间复杂度降低到 O(log n)。
R-tree（Rectangle Tree）是最常用的空间索引数据结构，由 Antonin Guttman 于 1984 年提出。
R-tree 的基本思想：
[根节点 MBR]
/ \
[分支 MBR_1] [分支 MBR_2]
/ | \ / |
[叶1] [叶2] [叶3] [叶4] [叶5]
这个"两阶段"过程大大减少了需要进行精确几何计算的要素数量。
GeoPandas（通过 Shapely）使用的是 STRtree（Sort-Tile-Recursive tree），它是一种批量加载的 R-tree 变体：
GeoPandas 的空间索引基于 Shapely 的 STRtree 实现，而 Shapely 又基于 GEOS 库的 C++ 实现。
import geopandas as gpd
from shapely.geometry import Point
import numpy as np
# 创建大量随机点
np.random.seed(42)
n = 10000
points = [Point(x, y) for x, y in zip(
np.random.uniform(0, 100, n),
np.random.uniform(0, 100, n)
)]
gdf = gpd.GeoDataFrame(geometry=points)
# 查看空间索引类型
sindex = gdf.sindex
print(f"空间索引类型: {type(sindex)}")
print(f"索引大小: {sindex.size}")
许多 GeoPandas 操作会自动使用空间索引：
sjoin() - 空间连接clip() - 空间裁剪overlay() - 空间叠加分析nearest() - 最近邻查找import geopandas as gpd
from shapely.geometry import Point, box
# 创建 GeoDataFrame
gdf = gpd.GeoDataFrame(geometry=[
Point(0, 0), Point(1, 1), Point(2, 2),
Point(3, 3), Point(4, 4), Point(5, 5),
])
# 访问空间索引
sindex = gdf.sindex
print(f"索引中的对象数: {sindex.size}")
print(f"索引是否有效: {sindex.is_empty == False}")
# GeoSeries 也有 sindex
gs = gdf.geometry
sindex_gs = gs.sindex
print(f"GeoSeries 索引大小: {sindex_gs.size}")
空间索引不会在创建 GeoDataFrame 时自动构建，而是在首次访问 .sindex 属性时才创建（延迟初始化）。
import geopandas as gpd
from shapely.geometry import Point
gdf = gpd.GeoDataFrame(geometry=[Point(i, i) for i in range(1000)])
# 此时索引尚未创建
print(f"索引是否存在: {gdf.has_sindex}") # False
# 首次访问时创建
sindex = gdf.sindex
print(f"索引是否存在: {gdf.has_sindex}") # True
has_sindex 属性可以检查空间索引是否已创建，而不会触发索引构建。
# 不触发创建的检查
if not gdf.has_sindex:
print("索引尚未创建，首次空间查询时将自动创建")
# 显式创建索引
_ = gdf.sindex # 触发创建
print(f"现在索引存在: {gdf.has_sindex}")
sindex.query(geometry, predicate=None, sort=False)
| 参数 | 说明 |
|---|---|
geometry |
查询几何对象 |
predicate |
空间谓词（如 'intersects', 'contains' 等） |
sort |
是否对结果排序 |
import geopandas as gpd
from shapely.geometry import Point, box
import numpy as np
# 创建数据
np.random.seed(42)
gdf = gpd.GeoDataFrame(
geometry=[Point(np.random.uniform(0, 100), np.random.uniform(0, 100))
for _ in range(10000)]
)
# 查询区域
query_box = box(20, 20, 40, 40)
# 仅 MBR 过滤（粗过滤）
candidates = gdf.sindex.query(query_box)
print(f"MBR 候选数: {len(candidates)}")
# 返回 MBR 与查询框相交的所有要素索引
# 使用谓词进行精确查询
from shapely.geometry import Polygon
query_circle = Point(50, 50).buffer(15)
# intersects 谓词
intersecting = gdf.sindex.query(query_circle, predicate='intersects')
print(f"相交的要素数: {len(intersecting)}")
# within 谓词
within_idx = gdf.sindex.query(query_circle, predicate='contains')
print(f"被圆包含的要素数: {len(within_idx)}")
| 谓词 | 说明 |
|---|---|
'intersects' |
相交 |
'within' |
在查询几何内部 |
'contains' |
包含查询几何 |
'overlaps' |
重叠 |
'crosses' |
交叉 |
'touches' |
接触 |
'covers' |
覆盖 |
'covered_by' |
被覆盖 |
'contains_properly' |
严格包含 |
'dwithin' |
距离范围内 |
None |
仅 MBR 过滤 |
# 查询结果是索引数组，可以直接用于筛选
idx = gdf.sindex.query(query_circle, predicate='intersects')
result = gdf.iloc[idx]
print(f"查询到 {len(result)} 个要素")
query_bulk() 对多个查询几何同时执行空间查询，返回配对的索引数组。这比循环调用 query() 更高效。
注意： 在较新版本的 GeoPandas 中，推荐使用 query() 直接传入几何数组。
import geopandas as gpd
from shapely.geometry import Point, box
import numpy as np
# 创建大量点
np.random.seed(42)
n = 50000
gdf = gpd.GeoDataFrame(
geometry=[Point(np.random.uniform(0, 100), np.random.uniform(0, 100))
for _ in range(n)]
)
# 多个查询区域
query_geoms = gpd.GeoSeries([
box(10, 10, 20, 20),
box(30, 30, 50, 50),
box(70, 70, 90, 90),
])
# 批量查询
input_idx, tree_idx = gdf.sindex.query(query_geoms, predicate='intersects')
print(f"总匹配数: {len(tree_idx)}")
# input_idx 表示查询几何的索引
# tree_idx 表示匹配的数据要素索引
for i in range(len(query_geoms)):
mask = input_idx == i
count = mask.sum()
print(f"查询区域 {i}: 匹配 {count} 个要素")
sjoin() 内部使用类似的批量查询机制：
# sjoin 本质上是批量查询 + 属性合并
result = gpd.sjoin(
gdf,
gpd.GeoDataFrame(geometry=query_geoms),
predicate='intersects'
)
nearest() 使用空间索引高效地查找最近的几何对象。
import geopandas as gpd
from shapely.geometry import Point
import numpy as np
# 创建数据点
np.random.seed(42)
gdf = gpd.GeoDataFrame(
geometry=[Point(np.random.uniform(0, 100), np.random.uniform(0, 100))
for _ in range(10000)]
)
# 查询点
query = Point(50, 50)
# 查找最近的要素
input_idx, tree_idx = gdf.sindex.nearest(query, return_all=False)
print(f"最近的要素索引: {tree_idx[0]}")
print(f"最近的要素: {gdf.iloc[tree_idx[0]].geometry}")
print(f"距离: {query.distance(gdf.iloc[tree_idx[0]].geometry):.4f}")
# 使用 sjoin_nearest 查找最近邻并保留距离
target = gpd.GeoDataFrame(geometry=[Point(50, 50)])
result = gpd.sjoin_nearest(target, gdf, max_distance=10, distance_col='dist')
result_sorted = result.sort_values('dist')
print(f"10 单位内最近的 5 个要素:")
print(result_sorted.head())
# 为每个居民区找到最近的医院
residential = gpd.read_file("residential.shp")
hospitals = gpd.read_file("hospitals.shp")
# 确保使用投影坐标系
residential_proj = residential.to_crs(epsg=32650)
hospitals_proj = hospitals.to_crs(epsg=32650)
# 最近邻连接
result = gpd.sjoin_nearest(
residential_proj, hospitals_proj,
distance_col='hospital_distance',
how='left'
)
print(f"平均到最近医院的距离: {result['hospital_distance'].mean():.0f} 米")
print(f"最远到最近医院的距离: {result['hospital_distance'].max():.0f} 米")
import geopandas as gpd
from shapely.geometry import Point
gdf = gpd.GeoDataFrame(geometry=[Point(0, 0)])
sindex = gdf.sindex
# 查看所有可用的查询谓词
print("可用谓词:")
for pred in sorted(sindex.valid_query_predicates):
print(f" - {pred}")
| 场景 | 推荐谓词 |
|---|---|
| 查找区域内的点 | 'intersects' 或 'within' |
| 查找重叠区域 | 'overlaps' |
| 查找相邻要素 | 'touches' |
| 查找包含关系 | 'contains' |
| 距离范围查询 | 'dwithin' |
sjoin() 内部使用空间索引加速空间连接：
import geopandas as gpd
from shapely.geometry import Point, Polygon
import numpy as np
# 创建点数据
np.random.seed(42)
points = gpd.GeoDataFrame({
'value': np.random.randint(1, 100, 10000),
'geometry': [Point(np.random.uniform(0, 10), np.random.uniform(0, 10))
for _ in range(10000)]
})
# 创建网格多边形
from shapely.geometry import box
grid = []
for x in range(10):
for y in range(10):
grid.append({'cell': f'{x}_{y}', 'geometry': box(x, y, x+1, y+1)})
grid_gdf = gpd.GeoDataFrame(grid)
# 空间连接（自动使用空间索引）
joined = gpd.sjoin(points, grid_gdf, predicate='within')
print(f"连接结果行数: {len(joined)}")
# 统计每个格子内的点数
cell_counts = joined.groupby('cell').size()
print(f"每格平均点数: {cell_counts.mean():.1f}")
# 最近邻空间连接
from shapely.geometry import Point
# 设施点
facilities = gpd.GeoDataFrame({
'type': ['医院', '学校', '超市'] * 10,
'geometry': [Point(np.random.uniform(0, 10), np.random.uniform(0, 10))
for _ in range(30)]
})
# 查询点
queries = gpd.GeoDataFrame({
'name': [f'位置{i}' for i in range(5)],
'geometry': [Point(np.random.uniform(0, 10), np.random.uniform(0, 10))
for _ in range(5)]
})
# 为每个查询点找到最近的设施
result = gpd.sjoin_nearest(queries, facilities, distance_col='distance')
print(result[['name', 'type', 'distance']])
当 GeoDataFrame 的几何列发生变化时，空间索引会自动失效，在下次访问时重建。
import geopandas as gpd
from shapely.geometry import Point
gdf = gpd.GeoDataFrame(geometry=[Point(i, i) for i in range(100)])
# 触发索引创建
_ = gdf.sindex
print(f"索引存在: {gdf.has_sindex}") # True
# 修改几何列
gdf.geometry = gdf.geometry.buffer(1)
# 索引自动失效
print(f"修改后索引存在: {gdf.has_sindex}") # False
# 下次访问自动重建
_ = gdf.sindex
print(f"重建后索引存在: {gdf.has_sindex}") # True
# 如果确定不再需要索引，可以避免不必要的创建
# 例如，在一系列几何修改操作之间，不要访问 sindex
# 批量修改时的最佳实践
gdf_temp = gdf.copy() # 在副本上操作
# 进行多次几何修改...
gdf_temp['geometry'] = gdf_temp.buffer(1)
gdf_temp['geometry'] = gdf_temp.simplify(0.1)
# 修改完成后，索引只需构建一次
result = gdf_temp.sindex.query(some_geometry, predicate='intersects')
import geopandas as gpd
from shapely.geometry import Point, box
import numpy as np
import time
# 创建不同规模的数据集
def benchmark(n_points):
np.random.seed(42)
points = gpd.GeoDataFrame(
geometry=[Point(np.random.uniform(0, 1000), np.random.uniform(0, 1000))
for _ in range(n_points)]
)
query = box(400, 400, 600, 600)
# 方法一：不使用空间索引（逐个检查）
start = time.time()
result1 = points[points.intersects(query)]
time_no_index = time.time() - start
# 方法二：使用空间索引
start = time.time()
idx = points.sindex.query(query, predicate='intersects')
result2 = points.iloc[idx]
time_with_index = time.time() - start
print(f"数据量: {n_points:>8d} | "
f"无索引: {time_no_index:.4f}s | "
f"有索引: {time_with_index:.4f}s | "
f"加速比: {time_no_index/max(time_with_index, 0.0001):.1f}x | "
f"结果数: {len(result1)}")
# 测试不同规模
for n in [1000, 5000, 10000, 50000, 100000]:
benchmark(n)
| 数据量 | 无索引 | 有索引 | 加速比 |
|---|---|---|---|
| 1,000 | 0.005s | 0.001s | 5x |
| 10,000 | 0.05s | 0.001s | 50x |
| 100,000 | 0.5s | 0.002s | 250x |
| 1,000,000 | 5s | 0.005s | 1000x |
数据量越大，空间索引的优势越明显。
import time
# 索引构建时间也需要考虑
def measure_index_build(n):
np.random.seed(42)
gdf = gpd.GeoDataFrame(
geometry=[Point(np.random.uniform(0, 1000), np.random.uniform(0, 1000))
for _ in range(n)]
)
start = time.time()
_ = gdf.sindex
build_time = time.time() - start
print(f"数据量 {n:>8d}: 索引构建时间 {build_time:.4f}s")
for n in [1000, 10000, 100000, 500000]:
measure_index_build(n)
import geopandas as gpd
import numpy as np
def process_in_chunks(large_gdf, query_geom, chunk_size=50000):
"""分块处理大数据集"""
results = []
n_chunks = (len(large_gdf) + chunk_size - 1) // chunk_size
for i in range(n_chunks):
start_idx = i * chunk_size
end_idx = min((i + 1) * chunk_size, len(large_gdf))
chunk = large_gdf.iloc[start_idx:end_idx]
# 每个块独立使用空间索引
idx = chunk.sindex.query(query_geom, predicate='intersects')
if len(idx) > 0:
results.append(chunk.iloc[idx])
if results:
return gpd.GeoDataFrame(pd.concat(results))
return gpd.GeoDataFrame()
import geopandas as gpd
from shapely.geometry import box
import numpy as np
def spatial_partition(gdf, n_partitions=10):
"""按空间范围分区"""
bounds = gdf.total_bounds # [minx, miny, maxx, maxy]
x_step = (bounds[2] - bounds[0]) / n_partitions
partitions = {}
for i in range(n_partitions):
x_min = bounds[0] + i * x_step
x_max = bounds[0] + (i + 1) * x_step
partition_box = box(x_min, bounds[1], x_max, bounds[3])
# 使用空间索引快速筛选
idx = gdf.sindex.query(partition_box, predicate='intersects')
partitions[i] = gdf.iloc[idx]
return partitions
# 简化几何以加速索引查询
gdf_simplified = gdf.copy()
gdf_simplified['geometry'] = gdf.simplify(tolerance=10)
# 用简化几何做粗过滤
coarse_idx = gdf_simplified.sindex.query(query_geom, predicate='intersects')
# 用原始几何做精确过滤
candidates = gdf.iloc[coarse_idx]
final_result = candidates[candidates.intersects(query_geom)]
# 不好的做法：每次循环都修改几何，导致重复构建索引
for buffer_dist in [100, 200, 500, 1000]:
buffered = gdf.buffer(buffer_dist)
# 每次 buffer 后创建新的 GeoSeries，索引需要重建
# 好的做法：预计算所有缓冲区
buffers = {d: gdf.buffer(d) for d in [100, 200, 500, 1000]}
| 策略 | 适用场景 | 预期加速 |
|---|---|---|
| 使用空间索引 | 所有空间查询 | 10-1000x |
| 分块处理 | 内存受限的大数据集 | 避免内存溢出 |
| 空间分区 | 并行处理 | 与核数成正比 |
| 几何简化 | 复杂几何的查询 | 2-10x |
| 预计算缓冲区 | 多次使用的缓冲区 | 避免重复计算 |
| 投影优化 | 跨投影的分析 | 减少投影转换 |
本章全面介绍了 GeoPandas 中的空间索引与查询优化。主要内容回顾：
| 概念 | 说明 |
|---|---|
| R-tree | 基于最小外接矩形的层次索引结构 |
| STRtree | 批量加载的 R-tree 变体，GeoPandas 使用的实现 |
.sindex |
访问空间索引的属性 |
has_sindex |
检查索引是否已创建（不触发创建） |
| 方法 | 说明 |
|---|---|
query(geom, predicate) |
按几何和谓词查询 |
query(geom_array, predicate) |
批量查询 |
nearest(geom) |
最近邻查询 |
valid_query_predicates |
查看可用谓词 |
| 策略 | 说明 |
|---|---|
| 使用空间索引 | 大数据集必须使用 |
| 使用 sjoin | 自动利用空间索引 |
| 分块处理 | 避免内存溢出 |
| 空间分区 | 支持并行处理 |
| 避免重复构建 | 注意索引失效机制 |
sjoin() 自动处理空间索引，比手动查询更简洁至此，我们已经完成了 GeoPandas 空间分析核心功能的学习。通过空间索引的优化，您可以高效地处理大规模地理空间数据，充分发挥 GeoPandas 的分析能力。
此内容由惯性聚合(RSS阅读器)自动聚合整理，仅供阅读参考。 原文来自 — 版权归原作者所有。