znlgis 博客

GIS开发与技术分享

第15章：空间索引与查询优化

在处理大规模地理空间数据时，空间查询的性能是一个关键问题。空间索引通过对几何对象的空间位置建立索引结构，可以显著加速空间查询操作。本章将全面介绍 GeoPandas 中空间索引的原理、使用方法和优化策略。

15.1 空间索引概述 - R-tree 原理

15.1.1 为什么需要空间索引

在没有空间索引的情况下，执行空间查询（如”找出与某区域相交的所有要素”）需要逐一检查每个要素，时间复杂度为 O(n)。当数据量达到数十万甚至数百万时，这种暴力搜索将变得非常缓慢。

空间索引通过预先组织数据的空间分布信息，将查询的平均时间复杂度降低到 O(log n)。

15.1.2 R-tree 数据结构

R-tree（Rectangle Tree）是最常用的空间索引数据结构，由 Antonin Guttman 于 1984 年提出。

R-tree 的基本思想：

  1. 使用最小外接矩形（MBR）来近似表示每个几何对象
  2. 将空间中临近的 MBR 分组，每组用一个更大的 MBR 包围
  3. 递归地对这些分组再次分组，形成树状层次结构
              [根节点 MBR]
             /            \
     [分支 MBR_1]     [分支 MBR_2]
      /    |    \       /    |
   [叶1] [叶2] [叶3] [叶4] [叶5]

15.1.3 R-tree 查询过程

  1. 从根节点开始，用查询几何的 MBR 与节点的 MBR 比较
  2. 只进入 MBR 相交的子节点（剪枝）
  3. 到达叶节点后，返回可能匹配的候选要素
  4. 对候选要素进行精确的几何判断（过滤）

这个”两阶段”过程大大减少了需要进行精确几何计算的要素数量。

15.1.4 STRtree

GeoPandas（通过 Shapely）使用的是 STRtree（Sort-Tile-Recursive tree），它是一种批量加载的 R-tree 变体：

  • 一次性构建索引，比逐个插入更高效
  • 查询性能通常优于标准 R-tree
  • 是只读索引，构建后不能添加或删除元素

15.2 GeoPandas 空间索引 - 基于 Shapely STRtree

15.2.1 底层实现

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}")

15.2.2 索引的自动使用

许多 GeoPandas 操作会自动使用空间索引：

  • sjoin() - 空间连接
  • clip() - 空间裁剪
  • overlay() - 空间叠加分析
  • nearest() - 最近邻查找

15.3 访问空间索引 - .sindex 属性

15.3.1 基本访问

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}")

15.3.2 GeoSeries 的空间索引

# GeoSeries 也有 sindex
gs = gdf.geometry
sindex_gs = gs.sindex
print(f"GeoSeries 索引大小: {sindex_gs.size}")

15.4 延迟创建机制 - 首次访问时创建，has_sindex 检查

15.4.1 延迟创建

空间索引不会在创建 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

15.4.2 has_sindex 属性

has_sindex 属性可以检查空间索引是否已创建，而不会触发索引构建。

# 不触发创建的检查
if not gdf.has_sindex:
    print("索引尚未创建，首次空间查询时将自动创建")

# 显式创建索引
_ = gdf.sindex  # 触发创建
print(f"现在索引存在: {gdf.has_sindex}")

15.5 query() 方法 - 按几何+谓词查询

15.5.1 基本语法

sindex.query(geometry, predicate=None, sort=False)

15.5.2 参数说明

参数 说明
geometry 查询几何对象
predicate 空间谓词（如 ‘intersects’, ‘contains’ 等）
sort 是否对结果排序

15.5.3 无谓词查询（仅 MBR 过滤）

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 与查询框相交的所有要素索引

15.5.4 带谓词查询（精确过滤）

# 使用谓词进行精确查询
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)}")

15.5.5 支持的谓词列表

谓词 说明
'intersects' 相交
'within' 在查询几何内部
'contains' 包含查询几何
'overlaps' 重叠
'crosses' 交叉
'touches' 接触
'covers' 覆盖
'covered_by' 被覆盖
'contains_properly' 严格包含
'dwithin' 距离范围内
None 仅 MBR 过滤

15.5.6 使用查询结果

# 查询结果是索引数组，可以直接用于筛选
idx = gdf.sindex.query(query_circle, predicate='intersects')
result = gdf.iloc[idx]
print(f"查询到 {len(result)} 个要素")

15.6 query_bulk() - 批量查询

15.6.1 基本概念

query_bulk() 对多个查询几何同时执行空间查询，返回配对的索引数组。这比循环调用 query() 更高效。

注意： 在较新版本的 GeoPandas 中，推荐使用 query() 直接传入几何数组。

15.6.2 代码示例

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} 个要素")

15.6.3 与 sjoin 的关系

sjoin() 内部使用类似的批量查询机制：

# sjoin 本质上是批量查询 + 属性合并
result = gpd.sjoin(
    gdf,
    gpd.GeoDataFrame(geometry=query_geoms),
    predicate='intersects'
)

15.7 nearest() - 最近邻查询

15.7.1 基本概念

nearest() 使用空间索引高效地查找最近的几何对象。

15.7.2 代码示例

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}")

15.7.3 查找 k 个最近邻

# 使用 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())

15.7.4 实际应用

# 为每个居民区找到最近的医院
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} 米")

15.8 valid_query_predicates - 可用谓词属性

15.8.1 查看可用谓词

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}")

15.8.2 谓词选择建议

场景 推荐谓词
查找区域内的点 'intersects' 或 'within'
查找重叠区域 'overlaps'
查找相邻要素 'touches'
查找包含关系 'contains'
距离范围查询 'dwithin'

15.9 空间索引在 sjoin 中的应用

15.9.1 sjoin 的工作原理

sjoin() 内部使用空间索引加速空间连接：

  1. 构建右表的空间索引（如果尚未构建）
  2. 对左表每个几何，使用空间索引查找候选匹配
  3. 对候选匹配进行精确的空间谓词判断
  4. 合并匹配要素的属性
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}")

15.9.2 sjoin_nearest

# 最近邻空间连接
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']])

15.10 索引重建策略 - 几何变化后自动重建

15.10.1 自动失效机制

当 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

15.10.2 手动管理

# 如果确定不再需要索引，可以避免不必要的创建
# 例如，在一系列几何修改操作之间，不要访问 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')

15.11 性能基准测试 - 有无空间索引的性能对比

15.11.1 基准测试代码

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)

15.11.2 典型性能对比

数据量 无索引 有索引 加速比
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

数据量越大，空间索引的优势越明显。

15.11.3 索引构建时间

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)

15.12 大数据集优化策略

15.12.1 分块处理

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()

15.12.2 空间分区

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

15.12.3 减少几何复杂度

# 简化几何以加速索引查询
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)]

15.12.4 避免重复索引构建

# 不好的做法：每次循环都修改几何，导致重复构建索引
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]}

15.12.5 综合优化建议

策略 适用场景 预期加速
使用空间索引 所有空间查询 10-1000x
分块处理 内存受限的大数据集 避免内存溢出
空间分区 并行处理 与核数成正比
几何简化 复杂几何的查询 2-10x
预计算缓冲区 多次使用的缓冲区 避免重复计算
投影优化 跨投影的分析 减少投影转换

15.13 本章小结

本章全面介绍了 GeoPandas 中的空间索引与查询优化。主要内容回顾：

空间索引基础

概念 说明
R-tree 基于最小外接矩形的层次索引结构
STRtree 批量加载的 R-tree 变体，GeoPandas 使用的实现
.sindex 访问空间索引的属性
has_sindex 检查索引是否已创建（不触发创建）

查询方法

方法 说明
query(geom, predicate) 按几何和谓词查询
query(geom_array, predicate) 批量查询
nearest(geom) 最近邻查询
valid_query_predicates 查看可用谓词

性能优化

策略 说明
使用空间索引 大数据集必须使用
使用 sjoin 自动利用空间索引
分块处理 避免内存溢出
空间分区 支持并行处理
避免重复构建 注意索引失效机制

使用建议

  1. 始终使用空间索引：对于超过几百个要素的数据集，空间索引都能带来显著加速
  2. 优先使用 sjoin：sjoin() 自动处理空间索引，比手动查询更简洁
  3. 注意索引生命周期：几何修改后索引会失效，避免频繁修改-查询交替
  4. 选择合适的谓词：不同谓词的性能可能不同，选择最精确的谓词
  5. 监控内存使用：超大数据集的空间索引本身也占用内存

至此，我们已经完成了 GeoPandas 空间分析核心功能的学习。通过空间索引的优化，您可以高效地处理大规模地理空间数据，充分发挥 GeoPandas 的分析能力。

↑