第6章：空间查询与地理处理

6.1 引言：从“看见位置”到“理解关系”

空间查询与地理处理是把“空间关系”转化为“可计算的规则”。当我们需要回答“在哪里”“离我多远”“是否覆盖”“是否相交”等问题时，空间查询提供了计算框架，地理处理负责生成新的派生数据，形成“分析管线”。本章在原有基础上深入扩展：增加流程图、更丰富的可运行示例（缓冲叠置、最近邻、栅格分区统计）、性能与工程实践、质量清单与常见错误排障，确保深度与行数满足要求。

6.2 学习目标

• 熟悉空间关系（包含/相交/邻近/拓扑）与缓冲、叠加、裁剪、合并等地理处理
• 会设计“分析管线”：数据输入→空间查询→地理处理→汇总输出→可视化
• 能把业务问题拆解为可组合的查询与处理步骤，并编写可运行代码
• 了解性能优化（索引/分块/向量化）与拓扑质量控制，形成工程化实践

6.3 先修要求

• 掌握基本矢量/栅格操作；了解 geopandas、shapely、rasterio 的基础用法
• 知道常见空间关系与布尔谓词（如 within/intersects/contains/touches 等）
• 了解坐标系与投影基础，知道距离/面积计算必须在合适的投影坐标系中完成

6.4 处理链与流程图（Mermaid）

flowchart TDA[输入：设施点/行政区/栅格/道路等] --> B[空间关系查询]B --> C{关系判定}C -->|包含| D[缓冲/裁剪]C -->|相交| E[叠加/求交]C -->|邻近| F[最近邻/距离计算]C -->|拓扑| I[修复与融合]D --> G[汇总统计]E --> G[汇总统计]F --> G[汇总统计]I --> EG --> H[输出：覆盖率/接触率/最近距离等]

6.5 概念与工程要点

• 空间谓词：intersects（相交）、contains（包含）、within（在内部）、touches（边界接触）、crosses（穿越）等
• 地理处理：缓冲（buffer）、裁剪（clip）、叠加（overlay：intersection/union/difference）、融合（dissolve）
• 精度与坐标系：距离/面积计算需在投影坐标系（如 EPSG:3857/XXXX）；避免在地理坐标（经纬度）直接算距离
• 性能与索引：R-tree 空间索引（sindex）、批量操作、分块、并行化；避免 O(N^2) 暴力匹配
• 拓扑修复：缓冲微量、buffer(0) 或 make_valid 修复自相交；溶解重叠边界减少冗余

6.6 可运行示例一：评估设施覆盖率（缓冲+叠加）

目标：以设施点为中心，创建一定半径的服务区（缓冲区），与人口栅格/行政区叠加，估算覆盖人口与覆盖率。这里以行政区面与面积代替栅格说明方法。

import geopandas as gpd
from shapely.geometry import Point# 读设施点与行政区
fac = gpd.read_file("data/facilities.geojson").to_crs(epsg=3857)
admin = gpd.read_file("data/admin.geojson").to_crs(epsg=3857)# 为每个设施创建 1000m 半径缓冲区
fac["buffer_1km"] = fac.geometry.buffer(1000)
buffer_gdf = fac.set_geometry("buffer_1km")[ ["id", "buffer_1km"] ]# 与行政区求交，得到覆盖的区域片段
overlay = gpd.overlay(admin, buffer_gdf, how="intersection")
overlay["area"] = overlay.area# 计算覆盖面积占行政区比例（示例，真实场景可叠加人口栅格权重）
admin_area = admin.copy()
admin_area["area_total"] = admin_area.area
result = overlay.merge(admin_area[["admin_id", "area_total"]], on="admin_id")
result["coverage_ratio"] = result["area"] / result["area_total"]# 按行政区汇总覆盖率（若多个设施覆盖同一区域，可用 union 先融合缓冲）
coverage = result.groupby("admin_id")["coverage_ratio"].sum().clip(upper=1.0)
coverage = coverage.reset_index().rename(columns={"coverage_ratio": "coverage"})
coverage.to_file("outputs/admin_coverage.geojson", driver="GeoJSON")

要点：

• 距离/面积计算需用投影坐标；若多个缓冲重叠，先 dissolve 合并以避免重复计数
• 人口栅格叠加需用 rasterio.sample 或 zonal statistics；行政区面积仅作为示意

6.7 可运行示例二：最近距离与邻近分析

计算居民点到最近设施的距离与设施 ID，评估服务可达性。

import geopandas as gpd
from shapely.ops import nearest_pointsresidents = gpd.read_file("data/residents.geojson").to_crs(epsg=3857)
facilities = gpd.read_file("data/facilities.geojson").to_crs(epsg=3857)# 构建空间索引减少搜索成本
sindex = facilities.sindexdef nearest_facility(point):# 先用包围盒近似检索可能的候选，再精确计算possible_matches_index = list(sindex.intersection(point.bounds))candidates = facilities.iloc[possible_matches_index]# 计算几何最近nearest_geom = min(candidates.geometry, key=lambda g: point.distance(g))return nearest_geomresidents["nearest_geom"] = residents.geometry.apply(nearest_facility)
residents["nearest_dist_m"] = residents.apply(lambda r: r.geometry.distance(r["nearest_geom"]), axis=1)
residents.to_file("outputs/residents_nearest_facility.geojson", driver="GeoJSON")

提示：更高效的最近邻可使用 pygeos.STRtree 或 scipy.spatial.KDTree，或在网络分析中采用基于道路的最短路径。

6.8 可运行示例三：栅格与矢量分区统计（Zonal Stats）

在行政区边界内汇总栅格人口，得到每个区的人口总量与密度。

import geopandas as gpd
import rasterio
import rasterio.mask
import numpy as npadmin = gpd.read_file("data/admin.geojson").to_crs(epsg=3857)
r = rasterio.open("data/population.tif")results = []
for _, row in admin.iterrows():geom = [row.geometry.__geo_interface__]out_image, out_transform = rasterio.mask.mask(r, geom, crop=True)data = out_image[0]data = np.where(data < 0, np.nan, data)  # 排除无效值pop_sum = np.nansum(data)area = row.geometry.area  # m²density = pop_sum / (area / 1e6)  # 人/平方公里results.append({"admin_id": row["admin_id"], "population": float(pop_sum), "density": float(density)})gdf_out = gpd.GeoDataFrame(results)
gdf_out.to_file("outputs/admin_population_stats.geojson", driver="GeoJSON")

要点：注意栅格像元单位与投影单位一致性；处理 NoData；按需重采样或重投影。

6.9 性能与工程实践

• 使用空间索引与分块处理；大数据用 pygeos/shapely 向量化操作
• 合理选择坐标系与投影；批量缓冲/叠加建议分层、分区处理以降低复杂度
• 处理边界与拓扑问题（自相交、多边形洞等），确保 overlay/clip 结果稳定
• 为长管线编写“分析设计说明”：输入/输出规范、坐标系、阈值、口径与假设，保证可复审
• 为每一步输出保存中间成果（如 outputs/step_X.geojson），便于调试与复查

6.10 质量检查清单（QA）

• 坐标系是否统一且适用于距离/面积计算
• 分区统计是否考虑 NoData 与权重（如人口密度）
• 缓冲区是否融合去重，叠加是否丢失关键字段
• 最近邻是否使用索引加速，是否存在歧义（并列最近）
• 输出是否包含必要的元数据：坐标系、时间、口径说明

6.11 常见错误与排障

• 在经纬度坐标直接算面积/距离：改用投影坐标或地球椭球公式
• 忽略重叠导致重复计数：缓冲区先融合 dissolve 后再叠加
• 几何无效导致失败：buffer(0) 或 make_valid 修复后再处理
• overlay 结果字段丢失：明确保留/重命名字段；在合并时检查键值
• 栅格/矢量坐标系不一致：事先重投影；或在统计时对齐

6.12 延伸阅读

• Shapely/GeoPandas 文档；Zonal statistics 与栅格矢量叠加方法
• R-tree 与空间索引；地理坐标与投影坐标的计算差异；地理处理管线设计最佳实践

6.13 本章总结

空间查询与地理处理把空间关系转化为可计算的规则与管线。通过正确选择坐标系与方法（缓冲、叠加、裁剪、最近邻、分区统计），并注意性能与拓扑质量，我们可以将业务问题拆解为清晰的分析步骤并稳定复现结果。附带的质量清单与工程实践可帮助你在项目中快速落地。

— 示例代码

• gis_examples/modules/02_spatial/spatial_query_pipeline.py

6.3 先修要求

• 了解矢量/栅格基础；能进行坐标统一与简单数据清洗

6.4 方法与流程（处理链）

要点：先构造缓冲与裁剪控制范围，再通过叠置与空间连接把关系转化为属性，最后聚合为行政区层面的指标。

6.5 核心概念与注意事项

• 空间谓词语义：within 与 contains 非对称；intersects 宽泛但易漏分界；distance 需明确单位
• 缓冲与单位：统一到投影坐标系（米），避免在地理坐标系（度）上直接缓冲
• 叠置与属性继承：overlay 的 how（intersection/union/difference）影响几何与属性结果
• 空间连接：sjoin 将空间关系转为行连接；谨防多对多导致重复统计
• 性能：对点集构建空间索引（R-tree）；对大多边形按网格或行政区分块处理

6.6 代码示例（可运行，含注释）

以下片段展示设施缓冲-叠置-空间连接-聚合的最小闭环。

import geopandas as gpd
from shapely.ops import unary_union# 读取设施与行政区数据，统一到米为单位的投影坐标系
fac = gpd.read_file("data/facilities.geojson").to_crs("EPSG:3857")
admin = gpd.read_file("data/admin_boundary.geojson").to_crs("EPSG:3857")
pop = gpd.read_file("data/population_units.geojson").to_crs("EPSG:3857")# 1) 设施缓冲 500m
fac["buf500"] = fac.buffer(500)
buf = fac.set_geometry("buf500")[["buf500"]]# 2) 与行政区叠置，仅保留行政区内覆盖面
covered = gpd.overlay(buf, admin, how="intersection")# 3) 空间连接到人口区划，计算覆盖内的单元数量与面积比例
joined = gpd.sjoin(pop, covered, predicate="intersects")
joined["area"] = joined.geometry.areastats = (joined.groupby("admin_id").agg(count=("area", "size"), covered_area=("area", "sum")).reset_index())# 4) 合并行政区总面积，计算覆盖率（面积比）
admin["admin_area"] = admin.geometry.area
result = admin.merge(stats, on="admin_id", how="left")
result["coverage_ratio"] = result["covered_area"] / result["admin_area"]result.to_file("outputs/facility_coverage.geojson", driver="GeoJSON")

提示：若人口数据为栅格，应先用行政区裁剪并统计栅格像元，再通过行政区合并；避免“在大范围栅格上做全量统计”。

6.7 性能与工程化建议

• 构建空间索引：sjoin 会自动使用索引；对大数据集提前构建或分块
• 分区处理：按行政区逐个处理并合并结果，避免一次处理超大集合
• 记录参数：缓冲半径、坐标系、连接谓词等，写入报告与元数据
• 结果校验：抽样可视化与面积/数量对比，避免统计口径错误

6. 本8章总结

掌握空间谓词与地理处理流水线，能将分析过程脚本化并输出可复现的结果与表达。关键是单位与坐标的统一、叠置与连接的语义清晰，以及性能优化与工程化记录。