空间数据不会Python处理？GIS二次开发与地理处理脚本实战手册（含：代码模板）

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

如果你正在找一份能直接上手的《空间数据不会Python处理？GIS二次开发与地理处理脚本实战手册（含：代码模板）》，这篇文章会从最常见的矢量数据读取、坐标系检查、字段处理、空间分析、批量导出开始，带你搭出一套可复用的 Python GIS 脚本框架。

很多 GIS 同学会用 QGIS 或 ArcGIS Pro 点工具，但一到批量处理几十个 Shapefile、自动裁剪行政区、批量计算面积、把结果写入 GeoPackage，就不知道怎么用 Python 处理空间数据。其实 GIS 二次开发并不一定要从复杂插件开始，先把常用地理处理流程写成脚本，才是最稳的入门路线。

引言：为什么 GIS 二次开发要先学地理处理脚本

GIS 二次开发通常有三类任务：软件扩展、自动化处理、WebGIS 服务开发。对大多数 GIS 学生和初级工程师来说，最先遇到的不是写一个完整系统，而是把重复的地理处理工作自动化。

例如：

• 每天收到一批 Shapefile，需要统一坐标系并合并。
• 项目要求按行政区批量裁剪道路、水系、建筑物数据。
• 需要把点位数据批量缓冲，并统计落入每个区域的数量。
• 要将 Excel 经纬度表转换为空间图层，再导出 GeoPackage 或 GeoJSON。
• ArcGIS Pro 模型构建器能跑，但希望改成 ArcPy 脚本方便复用。

这些场景都属于典型的空间数据 Python 处理问题。只要掌握读取、坐标系、属性字段、空间关系、批处理、结果验证这几件事，就可以写出相当实用的 GIS 二次开发脚本。

背景：空间数据 Python 处理常见工具怎么选

Python GIS 生态里工具很多，新手最容易卡在“我到底该用哪个库”。下面先把常用工具放在同一个表里，方便你根据工作场景选择。

工具	适合场景	典型用途	注意点
GeoPandas	矢量数据批处理、空间分析入门	读取 Shapefile、GeoPackage、GeoJSON，做叠加、缓冲、空间连接	处理超大数据时内存压力较大
Shapely	几何对象计算	判断相交、包含、距离、缓冲	通常配合 GeoPandas 使用
PyProj	坐标系转换	经纬度转投影坐标、EPSG 检查	面积和距离计算前必须确认坐标系
Rasterio	栅格数据处理	读取 GeoTIFF、裁剪栅格、重投影	适合栅格，不适合矢量属性表处理
GDAL/OGR	数据格式转换、底层地理处理	格式转换、投影转换、批量导出	安装和参数较复杂，但能力很强
ArcPy	ArcGIS Pro 自动化	调用 ArcGIS 工具箱、批量制图、企业环境脚本	依赖 ArcGIS Pro 授权和 Python 环境
QGIS Processing	QGIS 自动化处理	在 PyQGIS 中调用 QGIS 算法	需要正确加载 QGIS Python 环境

如果你刚开始做空间数据 Python 处理，建议先从 GeoPandas 开始；如果单位主要使用 ArcGIS Pro，就学习 ArcPy；如果你要做格式转换和大批量数据清洗，再逐步补 GDAL/OGR。

原理：写地理处理脚本前必须理解的 5 件事

1. 空间数据不是普通表格

空间数据除了普通字段，还有几何字段。几何字段保存点、线、面等空间对象。Python 脚本处理空间数据时，不能只看属性表，还要检查几何类型、坐标系、空间范围和拓扑质量。

2. 坐标系决定面积、距离和叠加结果是否可靠

很多新手脚本能运行，但结果是错的，根源通常是坐标系。经纬度坐标系适合表达位置，但不适合直接计算面积和距离。做缓冲区、面积统计、距离分析前，最好先投影到合适的平面坐标系。

3. 空间关系需要先保证数据在同一坐标系

空间连接、裁剪、相交、包含等分析，本质上是在比较几何对象的位置关系。如果两个图层坐标系不一致，结果可能为空、错位，或者看起来能跑但空间关系错误。

4. 批处理脚本要考虑输入、输出和日志

真正可用的 GIS 二次开发脚本，不应该只处理一个文件。它应该能批量遍历文件夹，自动检查格式，输出处理结果，并记录失败原因。

5. 结果必须验证

地理处理脚本不是运行成功就结束。至少要验证输出图层数量、坐标系、字段、几何有效性、空间范围是否符合预期。

步骤：用 GeoPandas 写一个可复用的空间数据处理脚本

下面以常见任务为例：读取一个面图层，统一坐标系，计算面积，筛选面积大于指定阈值的要素，并导出为 GeoPackage。

步骤 1：准备 Python 环境

推荐使用 Conda 创建独立环境，避免 GeoPandas、GDAL、Fiona、PyProj 之间的依赖冲突。

conda create -n gis-python python=3.11
conda activate gis-python
conda install -c conda-forge geopandas pyogrio shapely pyproj pandas

如果你使用的是 ArcGIS Pro，也可以使用 ArcGIS Pro 自带的 Python 环境，但不要随意破坏默认环境。更稳妥的方式是克隆环境后再安装依赖。

步骤 2：读取空间数据并查看基本信息

import geopandas as gpd

input_path = r"D:gis_projectdatadistrict.shp"

gdf = gpd.read_file(input_path)

print("要素数量：", len(gdf))
print("字段列表：", list(gdf.columns))
print("几何类型：", gdf.geom_type.unique())
print("坐标系：", gdf.crs)
print("空间范围：", gdf.total_bounds)

这一步不要省。很多脚本错误不是算法问题，而是输入数据本身有问题，比如没有坐标系、几何为空、字段名不一致、图层不是面数据。

步骤 3：检查并统一坐标系

如果要计算面积，建议转换为适合本地的投影坐标系。下面示例使用 EPSG:3857 只是演示，正式项目中应根据研究区选择合适的投影坐标系，例如 CGCS2000 高斯克吕格分带或当地常用投影。

target_crs = "EPSG:3857"

if gdf.crs is None:
    raise ValueError("输入数据缺少坐标系，请先在 GIS 软件中定义正确坐标系。")

gdf_proj = gdf.to_crs(target_crs)

print("转换后坐标系：", gdf_proj.crs)

这里要区分两个概念：定义坐标系和投影转换。定义坐标系只是告诉软件数据原来是什么坐标系；投影转换是把坐标真正转换到另一个坐标系。原始坐标系不清楚时，不能随便用 set_crs 硬指定。

步骤 4：计算面积并筛选要素

gdf_proj["area_m2"] = gdf_proj.geometry.area
gdf_proj["area_km2"] = gdf_proj["area_m2"] / 1_000_000

result = gdf_proj[gdf_proj["area_km2"] >= 10].copy()

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

如果你的数据仍然是经纬度坐标系，geometry.area 得到的不是平方米，而是“度的平方”，没有直接业务意义。这是空间数据 Python 处理中最常见的错误之一。

步骤 5：导出为 GeoPackage

output_path = r"D:gis_projectoutputdistrict_area_result.gpkg"

result.to_file(
    output_path,
    layer="district_area_result",
    driver="GPKG"
)

print("导出完成：", output_path)

相比 Shapefile，GeoPackage 更适合现代 GIS 工作流。它支持较长字段名、多个图层、中文路径兼容性更好，也更适合在 QGIS 和 Python 之间来回使用。

完整 GeoPandas 代码模板

import geopandas as gpd
from pathlib import Path

def process_polygon_area(input_path, output_path, target_crs, min_area_km2):
    input_path = Path(input_path)
    output_path = Path(output_path)

    if not input_path.exists():
        raise FileNotFoundError(f"输入文件不存在：{input_path}")

    gdf = gpd.read_file(input_path)

    if gdf.empty:
        raise ValueError("输入图层为空。")

    if gdf.crs is None:
        raise ValueError("输入图层缺少坐标系，请先定义正确坐标系。")

    if not all(gdf.geom_type.isin(["Polygon", "MultiPolygon"])):
        raise ValueError("当前脚本只适用于面图层。")

    gdf = gdf[gdf.geometry.notna()].copy()
    gdf = gdf[gdf.geometry.is_valid].copy()

    gdf_proj = gdf.to_crs(target_crs)

    gdf_proj["area_m2"] = gdf_proj.geometry.area
    gdf_proj["area_km2"] = gdf_proj["area_m2"] / 1_000_000

    result = gdf_proj[gdf_proj["area_km2"] >= min_area_km2].copy()

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

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

    return {
        "input_count": len(gdf),
        "output_count": len(result),
        "output_path": str(output_path)
    }

if __name__ == "__main__":
    info = process_polygon_area(
        input_path=r"D:gis_projectdatadistrict.shp",
        output_path=r"D:gis_projectoutputdistrict_area_result.gpkg",
        target_crs="EPSG:3857",
        min_area_km2=10
    )
    print(info)

步骤：批量处理一个文件夹里的 Shapefile

GIS 项目里很少只处理一个图层。下面这个模板用于批量读取文件夹中的 Shapefile，并统一输出为 GeoPackage。

import geopandas as gpd
from pathlib import Path

input_dir = Path(r"D:gis_projectshp")
output_dir = Path(r"D:gis_projectgpkg")
output_dir.mkdir(parents=True, exist_ok=True)

target_crs = "EPSG:3857"

for shp_path in input_dir.glob("*.shp"):
    try:
        print(f"正在处理：{shp_path.name}")

        gdf = gpd.read_file(shp_path)

        if gdf.empty:
            print(f"跳过空图层：{shp_path.name}")
            continue

        if gdf.crs is None:
            print(f"跳过无坐标系图层：{shp_path.name}")
            continue

        gdf = gdf[gdf.geometry.notna()].copy()
        gdf = gdf.to_crs(target_crs)

        output_path = output_dir / f"{shp_path.stem}.gpkg"

        gdf.to_file(
            output_path,
            layer=shp_path.stem,
            driver="GPKG"
        )

        print(f"完成：{output_path}")

    except Exception as e:
        print(f"处理失败：{shp_path.name}，原因：{e}")

这个模板的重点不是代码有多复杂，而是形成批处理脚本的基本习惯：遍历、检查、处理、导出、异常捕获。

步骤：用 Python 做空间连接统计

空间连接是空间数据 Python 处理中的高频任务。典型需求是：统计每个行政区内有多少个点位。

import geopandas as gpd

points_path = r"D:gis_projectdatapoi.shp"
polygons_path = r"D:gis_projectdatadistrict.shp"
output_path = r"D:gis_projectoutputdistrict_poi_count.gpkg"

points = gpd.read_file(points_path)
polygons = gpd.read_file(polygons_path)

if points.crs is None or polygons.crs is None:
    raise ValueError("点图层或面图层缺少坐标系。")

points = points.to_crs(polygons.crs)

joined = gpd.sjoin(
    points,
    polygons[["district_id", "district_name", "geometry"]],
    how="inner",
    predicate="within"
)

count_table = joined.groupby("district_id").size().reset_index(name="poi_count")

result = polygons.merge(count_table, on="district_id", how="left")
result["poi_count"] = result["poi_count"].fillna(0).astype(int)

result.to_file(output_path, layer="district_poi_count", driver="GPKG")

print("空间连接统计完成。")

这里的 predicate="within" 表示点在面内。常见空间关系还包括 intersects、contains、touches。选择哪个关系，取决于你的业务定义。

步骤：ArcGIS Pro 用户的 ArcPy 脚本模板

如果你的工作环境以 ArcGIS Pro 为主，ArcPy 是更直接的选择。ArcPy 可以调用 ArcGIS Pro 的地理处理工具，适合企业项目、模型自动化和批量制图。

import arcpy
from pathlib import Path

arcpy.env.overwriteOutput = True

input_fc = r"D:gis_projectdata.gdbdistrict"
output_fc = r"D:gis_projectoutput.gdbdistrict_buffer"
buffer_distance = "1000 Meters"

desc = arcpy.Describe(input_fc)

if desc.spatialReference.name == "Unknown":
    raise ValueError("输入数据坐标系未知，请先定义投影。")

arcpy.analysis.Buffer(
    in_features=input_fc,
    out_feature_class=output_fc,
    buffer_distance_or_field=buffer_distance,
    dissolve_option="NONE"
)

print("缓冲区分析完成：", output_fc)

ArcPy 的优势是能稳定调用 ArcGIS Pro 工具箱，缺点是运行环境依赖 ArcGIS Pro 授权。对于已有 ArcGIS 工作流的单位，ArcPy 通常比纯开源 Python 更容易接入生产流程。

常见坑：空间数据 Python 处理最容易出错的地方

1. 没有检查坐标系就计算面积和距离

这是最常见的问题。经纬度数据直接计算面积，结果通常没有实际意义。解决方式是先确认原始坐标系，再转换到合适的投影坐标系。

2. 把定义坐标系当成投影转换

set_crs 是定义坐标系，to_crs 是坐标转换。如果原数据坐标已经是经纬度，但你用 set_crs 强行指定为投影坐标系，数据会被错误解释。

3. Shapefile 字段名被截断

Shapefile 字段名通常存在长度限制，长字段名可能被截断。建议项目结果优先使用 GeoPackage、FileGDB 或 PostGIS。

4. 中文路径和编码导致读取失败

部分旧数据和旧环境对中文路径、中文字段、编码支持不好。建议项目目录使用英文路径，字段名尽量使用英文，中文名称放在字段值中。

5. 几何无效导致叠加分析失败

自相交面、空几何、重复节点都可能导致 overlay、clip、intersection 出错。处理前应检查 geometry.is_valid，必要时修复几何。

gdf["is_valid"] = gdf.geometry.is_valid
invalid = gdf[~gdf["is_valid"]]
print("无效几何数量：", len(invalid))

gdf["geometry"] = gdf.geometry.buffer(0)

buffer(0) 有时可以修复简单面几何问题，但不是万能方法。正式项目中仍建议在 QGIS、ArcGIS Pro 或 PostGIS 中进一步验证。

6. 空间连接结果重复

点落入多个重叠面时，空间连接会产生多条记录。这不是代码错误，而是面图层本身存在重叠或业务规则没有定义清楚。

方法比较：GeoPandas、ArcPy、QGIS Processing 怎么选

方法	优点	缺点	推荐使用场景
GeoPandas	开源、语法接近 Pandas、适合数据分析	超大数据性能有限，复杂拓扑处理需谨慎	学习 Python GIS、批量矢量处理、空间统计
ArcPy	直接调用 ArcGIS Pro 工具，企业兼容性好	依赖授权，跨平台能力弱	ArcGIS Pro 项目自动化、批量制图、企业生产环境
QGIS Processing	可调用 QGIS 大量算法，开源免费	独立脚本环境配置较麻烦	QGIS 工作流自动化、插件开发、开源桌面 GIS 项目
GDAL/OGR	格式支持强、命令行和脚本都好用	参数学习成本高	数据转换、投影转换、批量导入导出
PostGIS	适合大数据和多人协作，空间索引强	需要数据库基础	海量空间数据查询、WebGIS 后端、企业空间库

如果你的目标是入门 GIS 二次开发，推荐路线是：先学 GeoPandas 写脚本，再学 ArcPy 或 QGIS Processing 接入桌面软件，最后根据项目需要学习 PostGIS 和 WebGIS 服务。

检查清单：写完地理处理脚本后这样验收

• 输入数据是否存在：路径、文件名、图层名是否正确。
• 坐标系是否明确：原始 CRS 是否存在，目标 CRS 是否适合业务分析。
• 几何类型是否正确：点、线、面是否与脚本逻辑一致。
• 几何是否有效：是否存在空几何、自相交面、异常要素。
• 字段是否完整：连接字段、统计字段、导出字段是否存在。
• 空间范围是否正常：输出结果是否落在预期区域内。
• 数量是否合理：输入要素数、输出要素数、统计总数是否能对上。
• 单位是否正确：面积是平方米还是平方公里，距离是米还是度。
• 输出格式是否合适：临时结果可用 GeoPackage，生产环境可用 FileGDB 或 PostGIS。
• 异常是否可追踪：批处理时是否记录失败文件和失败原因。

FAQ：空间数据 Python 处理常见问题

Q1：不会 Python，可以直接学 GIS 二次开发吗？

可以，但建议先掌握 Python 基础语法，包括变量、列表、字典、函数、循环、异常处理和文件路径。GIS 二次开发不是只背 API，更重要的是能把业务流程拆成可执行步骤。

Q2：空间数据 Python 处理应该先学 GeoPandas 还是 ArcPy？

如果你没有 ArcGIS Pro 授权，或者想走开源路线，先学 GeoPandas。如果你所在单位主要使用 ArcGIS Pro，并且数据都在 FileGDB 或企业地理数据库中，先学 ArcPy 更合适。

Q3：为什么我的面积计算结果特别小？

大概率是因为数据仍然是经纬度坐标系。经纬度的单位是度，不是米。计算面积前应使用 to_crs 转换到合适的投影坐标系。

Q4：Python 读取 Shapefile 中文乱码怎么办？

优先检查数据编码和伴随的 .cpg 文件。实际项目中建议把 Shapefile 转为 GeoPackage，减少字段名、编码和多文件管理带来的问题。

Q5：GeoPandas 能替代 ArcGIS Pro 吗？

不能简单替代。GeoPandas 很适合脚本化矢量分析和数据清洗，但 ArcGIS Pro 在制图、企业工具箱、栅格分析、三维和专业扩展模块方面仍有优势。更现实的做法是两者配合使用。

Q6：脚本运行很慢怎么办？

先检查数据量、空间索引、坐标系转换次数和叠加分析复杂度。对于大数据，可以考虑分块处理、使用 PostGIS、减少不必要字段、提前裁剪研究区，或者用更适合大规模处理的数据库方案。

Q7：Python GIS 脚本如何和 QGIS 配合？

常见方式有两种：一种是在 Python 中用 GeoPandas 处理数据，再把结果加载到 QGIS 检查和制图；另一种是在 PyQGIS 中调用 QGIS Processing 算法，实现完整的 QGIS 自动化流程。

Q8：地理处理脚本需要写日志吗？

建议写。尤其是批量处理空间数据时，日志能告诉你哪个文件成功、哪个文件失败、失败原因是什么。项目越大，日志越重要。

结论：从一个脚本开始建立 GIS 二次开发能力

空间数据不会 Python 处理并不可怕，真正需要避免的是一开始就追求复杂系统。更有效的学习方式，是先把一个具体地理处理任务写成脚本：读取数据、检查坐标系、处理几何、执行分析、导出结果、验证成果。

当你能稳定完成批量转换、面积计算、空间连接、缓冲区分析、裁剪叠加这些任务后，GIS 二次开发就不再只是“写代码”，而会变成解决实际空间数据问题的工具箱。

建议你把本文中的代码模板保存成自己的项目脚手架。以后遇到新的空间数据处理任务，只需要替换输入路径、目标坐标系、分析逻辑和输出格式，就能快速搭建一套可复用的 Python GIS 工作流。