空间数据筛选效率低？GeoPandas实战技巧与完整代码案例（附：shp数据处理脚本）

wowwwai GIS研习社 · 工具流程与项目排障

如果你正在处理“空间数据筛选效率低？GeoPandas实战技巧与完整代码案例（附：shp数据处理脚本）”这个问题，通常不是因为 GeoPandas 不能做空间筛选，而是因为数据读取、坐标系、空间索引、筛选顺序和输出格式没有配合好。本文用一个可复用的 shp 数据处理脚本，演示如何把“慢、乱、容易报错”的空间数据筛选流程整理成可检查、可复现的 GeoPandas 实战流程。

引言：GeoPandas空间数据筛选为什么会慢

在 GIS 项目中，常见需求包括：从全国行政区中筛选某个省、市范围内的要素，从道路数据中提取研究区内的路网，或者从 POI 点数据中筛选落在指定边界内的记录。很多同学会直接使用 gpd.read_file() 读入所有数据，再用 within、intersects 或 clip 做筛选。

这种写法在小数据上没问题，但一旦 shp 文件达到几十万甚至上百万条记录，就很容易出现运行慢、内存占用高、结果为空、坐标偏移等问题。

本文重点解决一个具体问题：如何用 GeoPandas 高效筛选 shp 空间数据，并提供一份可直接改路径使用的完整 Python 脚本。

背景：典型的shp数据处理场景

假设我们有两个 shp 文件：

• 目标数据：例如道路、建筑物、POI、地块等，需要从中筛选要素。
• 研究区边界：例如行政区、规划范围、项目红线，用来限制筛选范围。

最常见的目标是：提取所有与研究区相交的要素，并输出为新的 shp 或 GeoPackage 文件。

很多初学者会写出类似代码：

import geopandas as gpd

data = gpd.read_file("roads.shp")
mask = gpd.read_file("study_area.shp")

result = data[data.intersects(mask.geometry.iloc[0])]
result.to_file("result.shp", encoding="utf-8")

这段代码在逻辑上看起来没问题，但实际项目中可能存在以下问题：

• 目标图层和研究区图层坐标系不一致。
• 研究区有多个面，直接取 iloc[0] 会漏掉其他区域。
• 没有先使用边界框或空间索引缩小候选数据。
• 字段编码、字段名长度、几何无效导致输出 shp 报错。
• 数据量太大时，逐条几何判断速度很慢。

原理：提高GeoPandas空间数据筛选效率的关键

GeoPandas 空间数据筛选的核心不是单纯换一个函数，而是减少真正参与几何计算的要素数量。空间计算通常比普通属性筛选更耗时，因此推荐使用“先粗筛、再精筛”的思路。

1. 坐标系必须一致

坐标系通常用 CRS 表示。两个图层如果 CRS 不一致，即使在地图上看起来位置相近，程序判断时也可能完全不相交。空间筛选前必须检查：

print(data.crs)
print(mask.crs)

如果不一致，应使用 to_crs() 将其中一个图层转换到另一个图层的 CRS。

2. 先用边界框减少候选要素

边界框筛选是一种粗筛。它只判断要素的外接矩形是否与研究区范围相交，速度通常比完整几何关系判断更快。GeoPandas 读取文件时可以使用 bbox 参数，尽量减少读入的数据量。

3. 使用空间索引加速几何关系判断

空间索引可以理解为空间数据的“目录”。它会把大量几何对象按位置组织起来，避免每个要素都和研究区逐一比较。GeoPandas 的 sjoin、clip 等操作会在可用时使用空间索引。

4. 属性预筛选应放在空间筛选之前

如果你只需要某类道路、某个年份的数据、某种用地类型，应先做属性筛选，再做空间筛选。这样可以显著减少空间判断的数据量。

步骤：GeoPandas高效筛选shp完整代码案例

下面给出一个完整的 shp 数据处理脚本。这个脚本适合以下任务：

• 读取一个待筛选 shp 文件。
• 读取一个研究区边界 shp 文件。
• 检查并统一坐标系。
• 修复常见无效几何。
• 进行可选属性预筛选。
• 使用空间连接筛选相交要素。
• 导出结果文件。

1. 安装依赖

建议在独立的 Python 环境中安装 GeoPandas。常用安装方式如下：

pip install geopandas pyogrio shapely rtree

如果使用 Conda，也可以使用：

conda install -c conda-forge geopandas pyogrio shapely rtree

说明：不同环境中 GeoPandas 对空间索引的支持依赖可能不同。若运行 sjoin 时提示空间索引相关错误，优先检查 rtree 或 shapely 是否安装正常。

2. 完整shp数据处理脚本

import geopandas as gpd
from pathlib import Path

# =========================
# 1. 参数配置
# =========================

input_shp = Path(r"data/roads.shp")
mask_shp = Path(r"data/study_area.shp")
output_file = Path(r"output/roads_in_study_area.shp")

# 可选：属性预筛选
# 例如只筛选道路等级为 primary 和 secondary 的要素
attribute_filter_field = "road_type"
attribute_filter_values = ["primary", "secondary"]

# 如果不需要属性筛选，设置为 None
# attribute_filter_field = None
# attribute_filter_values = None

# 空间关系：常用 intersects、within、contains
spatial_predicate = "intersects"

# 输出编码
output_encoding = "utf-8"


# =========================
# 2. 函数定义
# =========================

def read_vector(path):
    if not path.exists():
        raise FileNotFoundError(f"文件不存在：{path}")
    gdf = gpd.read_file(path)
    if gdf.empty:
        raise ValueError(f"图层为空：{path}")
    if gdf.crs is None:
        raise ValueError(f"图层缺少坐标系信息，请先定义 CRS：{path}")
    return gdf


def fix_invalid_geometry(gdf):
    gdf = gdf.copy()
    gdf = gdf[~gdf.geometry.is_empty]
    gdf = gdf[gdf.geometry.notnull()]
    gdf["geometry"] = gdf.geometry.make_valid()
    return gdf


def dissolve_mask(mask_gdf):
    mask_gdf = mask_gdf.copy()
    mask_gdf["__dissolve__"] = 1
    dissolved = mask_gdf.dissolve(by="__dissolve__").reset_index(drop=True)
    return dissolved


# =========================
# 3. 读取研究区
# =========================

print("读取研究区边界...")
mask = read_vector(mask_shp)
mask = fix_invalid_geometry(mask)
mask = dissolve_mask(mask)

print("研究区 CRS：", mask.crs)
print("研究区范围：", mask.total_bounds)


# =========================
# 4. 使用 bbox 粗筛读取目标数据
# =========================

print("根据研究区边界框读取目标数据...")
bbox = tuple(mask.total_bounds)

data = gpd.read_file(input_shp, bbox=bbox)

if data.empty:
    raise ValueError("bbox 粗筛后没有读取到目标要素，请检查坐标系或数据范围。")

if data.crs is None:
    raise ValueError(f"目标图层缺少坐标系信息，请先定义 CRS：{input_shp}")

print("目标数据 CRS：", data.crs)
print("bbox 粗筛后要素数：", len(data))


# =========================
# 5. 统一坐标系
# =========================

if data.crs != mask.crs:
    print("目标数据 CRS 与研究区 CRS 不一致，正在转换...")
    data = data.to_crs(mask.crs)

data = fix_invalid_geometry(data)


# =========================
# 6. 属性预筛选
# =========================

if attribute_filter_field and attribute_filter_values:
    if attribute_filter_field not in data.columns:
        raise ValueError(f"字段不存在：{attribute_filter_field}")

    before_count = len(data)
    data = data[data[attribute_filter_field].isin(attribute_filter_values)]
    print(f"属性预筛选：{before_count} -> {len(data)}")

    if data.empty:
        raise ValueError("属性预筛选后结果为空，请检查字段名和筛选值。")


# =========================
# 7. 空间筛选
# =========================

print("执行空间筛选...")
result = gpd.sjoin(
    data,
    mask[["geometry"]],
    how="inner",
    predicate=spatial_predicate
)

# 删除空间连接产生的辅助字段
if "index_right" in result.columns:
    result = result.drop(columns=["index_right"])

# 去重：一个要素可能与多个研究区面相交
result = result.drop_duplicates()

print("空间筛选后要素数：", len(result))

if result.empty:
    raise ValueError("空间筛选结果为空，请检查坐标系、空间关系和研究区范围。")


# =========================
# 8. 导出结果
# =========================

output_file.parent.mkdir(parents=True, exist_ok=True)

print("导出结果...")
result.to_file(output_file, encoding=output_encoding)

print(f"处理完成：{output_file}")

3. 如何改成筛选落在研究区内部的要素

如果你的需求不是“相交”，而是“完全落在研究区内”，可以把参数改成：

spatial_predicate = "within"

需要注意：对于线和面数据，within 比 intersects 更严格。道路只要有一小段超出研究区，就不会被保留。如果你希望按研究区边界裁剪道路，应使用 clip。

4. 如果要按边界裁剪几何

sjoin 只是筛选要素，不会改变几何形状。也就是说，一条道路只要和研究区相交，整条道路都会保留下来。如果你希望输出结果严格被研究区边界切开，可以使用：

clipped = gpd.clip(data, mask)
clipped.to_file("output/roads_clip.shp", encoding="utf-8")

因此，GeoPandas 空间数据筛选要先区分两个概念：

• 筛选：保留满足空间关系的原始要素，几何不被切割。
• 裁剪：用边界切割几何，只保留边界内部分。

常见坑：GeoPandas筛选shp结果为空或很慢

1. 坐标系看起来一样，实际 EPSG 不一致

有些数据在 GIS 软件中能叠加显示，是因为软件进行了动态投影。但 Python 代码不会自动替你修正所有 CRS 问题。运行空间筛选前，应打印并检查：

print(data.crs)
print(mask.crs)

如果一个是 EPSG:4326，另一个是投影坐标系，例如 CGCS2000 高斯投影，就必须统一后再处理。

2. shp字段名被截断

Shapefile 对字段名长度有限制。导出 shp 时，字段名可能被截断，导致后续脚本找不到字段。对于字段较多或字段名较长的数据，建议优先输出 GeoPackage：

result.to_file("output/result.gpkg", layer="result", driver="GPKG")

3. 忽略无效几何

面要素自相交、空几何、损坏几何都会影响空间筛选。脚本中的 make_valid() 可以修复一部分常见问题，但不是万能的。如果修复后仍然失败，需要在 QGIS 或 ArcGIS Pro 中进一步检查拓扑问题。

4. 直接读取全量大文件

如果目标 shp 很大，不建议一开始就完整读入。可以先使用 bbox 参数按照研究区范围粗筛，再进行空间连接。这样可以减少内存压力。

5. 把筛选和裁剪混为一谈

如果你发现筛选后的道路、河流仍然超出研究区边界，不一定是错误。因为 sjoin 和 intersects 默认只是筛选，不会裁剪几何。需要改变几何形状时，应使用 clip。

方法比较：sjoin、clip、cx和直接几何判断怎么选

对于大多数 shp 数据处理任务，推荐组合是：

1. 用研究区 total_bounds 做 bbox 粗筛。
2. 统一 CRS。
3. 属性预筛选。
4. 使用 sjoin 做空间筛选。
5. 如需切割几何，再使用 clip。

检查清单：运行GeoPandas空间数据筛选前先看这些

• 文件是否存在：路径中不要混用中文特殊符号、空格和错误扩展名。
• CRS 是否为空：如果 gdf.crs 是 None，先定义正确坐标系，不要盲目转换。
• 两个图层 CRS 是否一致：不一致时使用 to_crs() 转换。
• 研究区是否有多个面：多个面建议先 dissolve，避免漏选或重复。
• 是否需要属性预筛选：先筛字段，再做空间关系判断。
• 是否需要裁剪几何：需要裁剪就用 clip，只筛选就用 sjoin。
• 输出格式是否合适：字段复杂时优先用 GeoPackage，而不是 shp。
• 结果是否为空：优先检查 CRS、空间范围、筛选字段和值。

FAQ：GeoPandas实战常见问题

1. GeoPandas空间数据筛选效率低，最先应该优化哪里？

优先检查是否读取了全量大文件。如果目标数据很大，应先用研究区边界框进行 bbox 粗筛，再做属性预筛选和空间筛选。不要一开始就对所有要素做完整几何关系判断。

2. GeoPandas筛选shp结果为空怎么办？

先打印两个图层的 CRS 和范围：

print(data.crs, data.total_bounds)
print(mask.crs, mask.total_bounds)

如果 CRS 不一致或范围明显不重叠，空间筛选结果为空是正常的。还要检查属性筛选条件是否过严，例如字段值大小写不一致、字段名被 shp 截断等。

3. sjoin和clip有什么区别？

sjoin 用于筛选满足空间关系的要素，通常不改变原始几何。clip 会用边界裁剪几何，只保留边界内部分。如果你只是想找出研究区内有哪些要素，用 sjoin；如果你要生成严格落在研究区内的数据，用 clip。

4. 为什么导出的shp中文字段或属性乱码？

shp 对编码支持比较复杂，不同软件读取时可能解释不一致。可以在导出时指定编码：

result.to_file("result.shp", encoding="utf-8")

如果仍然乱码，建议导出为 GeoPackage：

result.to_file("result.gpkg", layer="result", driver="GPKG")

5. 空间索引需要手动创建吗？

多数情况下，使用 GeoPandas 的 sjoin、clip 等方法时会自动利用可用的空间索引。但如果环境缺少相关依赖，可能会报错或性能下降。建议安装 rtree，并保持 GeoPandas、Shapely 等依赖处于兼容版本。

结论：把GeoPandas筛选流程写成可复用脚本

GeoPandas 空间数据筛选效率低，通常不是单个函数的问题，而是整个处理流程的问题。正确做法是：先检查 CRS，再用 bbox 粗筛减少数据量，随后进行属性预筛选和空间索引支持的空间筛选，最后根据需求选择输出 shp 或 GeoPackage。

本文提供的 shp 数据处理脚本可以作为项目模板使用。你只需要替换输入路径、研究区路径、字段筛选条件和输出路径，就能快速完成道路、POI、地块、建筑物等常见空间数据筛选任务。

实战建议：如果结果用于后续空间分析，优先使用 GeoPackage 保存中间结果；如果只是与传统 GIS 软件交换数据，再导出 shp。