还在用ArcGIS？GeoPandas官方文档实操详解（附：完整代码）

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

还在用ArcGIS？GeoPandas官方文档实操详解（附：完整代码）这篇文章面向已经会一点 GIS、但想把重复制图和空间处理流程迁移到 Python 的读者。我们不讨论“ArcGIS 是否应该被替代”这种泛泛话题，而是用 GeoPandas 官方文档中最常用的一组能力，完成一次可复现的矢量数据读取、坐标系检查、属性筛选、空间连接、面积计算和结果导出流程。

引言：为什么 GIS 用户要学 GeoPandas

很多 GIS 学生和初级 GIS 工程师习惯在 ArcGIS 或 QGIS 里点工具箱：裁剪、叠加、空间连接、字段计算、导出结果。这个方式直观，但当任务需要重复处理几十个城市、几百个图层，或者每天自动更新一次结果时，纯 GUI 操作就会变得低效且难以复查。

GeoPandas 是 Python 生态中最常用的矢量 GIS 数据处理库之一。它把 Pandas 的表格处理能力和 Shapely 的几何计算能力结合起来，让我们可以像处理 Excel 表一样处理 Shapefile、GeoPackage、GeoJSON 等空间数据。

本文的核心目标很明确：用 GeoPandas 完成一套与 ArcGIS 常见矢量处理工具类似的工作流，并给出完整代码，方便你直接改成自己的项目脚本。

背景：这次实操要解决什么问题

假设我们有两类常见 GIS 数据：

• 一个行政区面图层，例如 city_boundary.gpkg，表示城市或区县边界。
• 一个兴趣点点图层，例如 poi.gpkg，表示学校、医院、商场、公交站等点位。

我们希望完成下面几个任务：

1. 读取矢量数据并查看字段、几何类型和坐标系。
2. 把数据统一到适合面积和距离计算的投影坐标系。
3. 筛选目标类型的 POI，例如只保留医院。
4. 统计每个行政区内有多少个目标 POI。
5. 计算行政区面积，并得到单位面积 POI 密度。
6. 把结果导出为 GeoPackage 或 Shapefile，供 ArcGIS、QGIS 或 WebGIS 后续使用。

这个流程对应 ArcGIS 里的多个工具：投影、按属性选择、空间连接、字段计算、导出要素。使用 GeoPandas 后，它们可以被组织成一段清晰的 Python 脚本。

原理：GeoPandas 和 ArcGIS 工具的对应关系

GeoPandas 的核心数据结构是 GeoDataFrame。你可以把它理解为“带 geometry 字段的 Pandas 表格”。每一行是一条空间要素，每一列是属性字段，其中 geometry 列保存点、线、面几何对象。

对于从 ArcGIS 转过来的用户，可以先建立下面这张对应表：

ArcGIS 常见操作	GeoPandas 对应方法	说明
添加数据	gpd.read_file()	读取 Shapefile、GeoPackage、GeoJSON 等矢量数据
查看图层属性表	gdf.head()	查看前几行属性和 geometry
定义/查看坐标系	gdf.crs	查看数据当前坐标参考系统
投影	gdf.to_crs()	把数据转换到另一个坐标系
按属性选择	gdf[gdf["字段"] == "值"]	用 Pandas 条件筛选属性记录
空间连接	gpd.sjoin()	按 contains、within、intersects 等空间关系连接属性
字段计算器	gdf["新字段"] = ...	新增字段并计算面积、密度、分类值
导出要素	gdf.to_file()	导出为 GeoPackage、Shapefile、GeoJSON 等格式

理解这个对应关系后，GeoPandas 官方文档中的很多示例就不再抽象。它本质上是在用代码表达 GIS 软件工具箱中的处理逻辑。

步骤：GeoPandas官方文档实操完整流程

步骤一：准备 Python 环境

建议使用 Conda 环境安装 GeoPandas。GeoPandas 依赖 GEOS、PROJ、GDAL、pyogrio、Shapely 等空间库，直接用 pip 在某些系统上可能遇到编译或动态库问题。

conda create -n geopandas-demo python=3.11
conda activate geopandas-demo
conda install -c conda-forge geopandas matplotlib pyogrio

如果你习惯用 Jupyter Notebook，可以继续安装：

conda install -c conda-forge jupyterlab

安装完成后，在 Python 中验证：

import geopandas as gpd

print(gpd.__version__)

步骤二：读取矢量数据

GeoPandas 读取矢量数据使用 read_file()。这里假设你的数据目录结构如下：

project/
  data/
    city_boundary.gpkg
    poi.gpkg
  output/
  analysis.py

读取代码如下：

import geopandas as gpd

boundary_path = "data/city_boundary.gpkg"
poi_path = "data/poi.gpkg"

boundary = gpd.read_file(boundary_path)
poi = gpd.read_file(poi_path)

print(boundary.head())
print(poi.head())

print(boundary.geometry.geom_type.value_counts())
print(poi.geometry.geom_type.value_counts())

print("boundary CRS:", boundary.crs)
print("poi CRS:", poi.crs)

这里最重要的检查不是字段名，而是 geometry 类型 和 CRS。行政区边界应当是 Polygon 或 MultiPolygon，POI 应当是 Point。如果几何类型不符合预期，后面的空间连接结果就可能完全错误。

步骤三：统一坐标系

在 GeoPandas 中，crs 表示坐标参考系统。常见的 WGS 84 经纬度坐标系是 EPSG:4326，单位是度，不适合直接计算面积和距离。

如果要计算面积、长度或密度，应转换到投影坐标系。例如中国区域常见做法是根据研究区选择合适的高斯克吕格、UTM 或地方投影。为了示范，这里使用 EPSG:3857 作为通用投影示例，但正式项目中不建议无脑使用它做精确面积统计。

target_crs = "EPSG:3857"

boundary_proj = boundary.to_crs(target_crs)
poi_proj = poi.to_crs(target_crs)

print(boundary_proj.crs)
print(poi_proj.crs)

如果你的数据已经是适合本地分析的投影坐标系，也可以不转换。但必须确认两个图层的 CRS 一致，否则空间连接、叠加分析和裁剪结果都不可信。

步骤四：按属性筛选目标 POI

假设 POI 图层中有一个字段叫 category，其中医院的值是 hospital。可以这样筛选：

hospital = poi_proj[poi_proj["category"] == "hospital"].copy()

print(hospital.shape)
print(hospital.head())

如果字段值不确定，先查看唯一值：

print(poi_proj["category"].value_counts().head(20))

中文数据常见字段值可能是“医院”“综合医院”“医疗卫生”等。实际项目中，不要直接假设字段值，先打印统计结果再写筛选条件。

步骤五：空间连接，统计每个行政区内的 POI 数量

空间连接是 GeoPandas 实操中最常用的功能之一。这里我们要判断每个医院点落在哪个行政区面内。

joined = gpd.sjoin(
    hospital,
    boundary_proj[["name", "geometry"]],
    how="inner",
    predicate="within"
)

print(joined.head())

这段代码表示：把医院点与行政区面做空间连接，只保留落在行政区内部的医院点，并把行政区的 name 字段连接到医院点上。

然后按行政区名称统计数量：

hospital_count = joined.groupby("name").size().reset_index(name="hospital_count")

print(hospital_count.head())

最后把统计表合并回行政区图层：

result = boundary_proj.merge(hospital_count, on="name", how="left")

result["hospital_count"] = result["hospital_count"].fillna(0).astype(int)

print(result[["name", "hospital_count"]].head())

步骤六：计算面积和单位面积密度

投影坐标系下的 geometry.area 通常返回平方米。我们可以把它转换为平方公里，并计算每平方公里医院数量。

result["area_km2"] = result.geometry.area / 1_000_000

result["hospital_density"] = result["hospital_count"] / result["area_km2"]

print(result[["name", "area_km2", "hospital_count", "hospital_density"]].head())

如果某些行政区面积为 0 或几何异常，需要先处理，否则密度字段可能出现无穷大或空值。

result = result[result["area_km2"] > 0].copy()
result["hospital_density"] = result["hospital_count"] / result["area_km2"]

步骤七：导出结果

推荐优先导出 GeoPackage，因为它支持较长字段名、中文路径相对友好，并且能在一个文件中保存多个图层。

output_path = "output/hospital_density.gpkg"

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

print("已导出：", output_path)

如果必须交付 Shapefile，也可以导出：

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

但要注意，Shapefile 对字段名长度、编码、单文件大小和几何类型都有较多限制。正式项目中，GeoPackage 通常更稳妥。

完整代码

下面是把前面所有步骤合并后的完整脚本。你只需要修改数据路径、字段名、目标坐标系和筛选条件。

import geopandas as gpd

# 1. 输入数据
boundary_path = "data/city_boundary.gpkg"
poi_path = "data/poi.gpkg"

# 2. 读取数据
boundary = gpd.read_file(boundary_path)
poi = gpd.read_file(poi_path)

print("边界图层记录数：", len(boundary))
print("POI 图层记录数：", len(poi))
print("边界 CRS：", boundary.crs)
print("POI CRS：", poi.crs)

print("边界几何类型：")
print(boundary.geometry.geom_type.value_counts())

print("POI 几何类型：")
print(poi.geometry.geom_type.value_counts())

# 3. 检查 CRS
if boundary.crs is None:
    raise ValueError("边界图层缺少 CRS，请先确认并定义坐标系。")

if poi.crs is None:
    raise ValueError("POI 图层缺少 CRS，请先确认并定义坐标系。")

# 4. 统一到目标投影坐标系
# 示例使用 EPSG:3857。正式面积统计请替换为适合研究区的投影坐标系。
target_crs = "EPSG:3857"

boundary_proj = boundary.to_crs(target_crs)
poi_proj = poi.to_crs(target_crs)

# 5. 清理无效几何
boundary_proj = boundary_proj[boundary_proj.geometry.notna()].copy()
poi_proj = poi_proj[poi_proj.geometry.notna()].copy()

boundary_proj = boundary_proj[~boundary_proj.geometry.is_empty].copy()
poi_proj = poi_proj[~poi_proj.geometry.is_empty].copy()

# 6. 按属性筛选医院 POI
# 请根据你的真实字段名和值修改 category 和 hospital
print("POI 类型统计：")
print(poi_proj["category"].value_counts().head(20))

hospital = poi_proj[poi_proj["category"] == "hospital"].copy()

print("医院 POI 数量：", len(hospital))

# 7. 空间连接：医院点落入哪个行政区
# 请确认行政区名称字段为 name，如不是，请替换为你的字段名
joined = gpd.sjoin(
    hospital,
    boundary_proj[["name", "geometry"]],
    how="inner",
    predicate="within"
)

print("成功匹配到行政区的医院数量：", len(joined))

# 8. 按行政区统计医院数量
hospital_count = joined.groupby("name").size().reset_index(name="hospital_count")

# 9. 合并回行政区图层
result = boundary_proj.merge(hospital_count, on="name", how="left")
result["hospital_count"] = result["hospital_count"].fillna(0).astype(int)

# 10. 计算面积和密度
result["area_km2"] = result.geometry.area / 1_000_000
result = result[result["area_km2"] > 0].copy()
result["hospital_density"] = result["hospital_count"] / result["area_km2"]

# 11. 查看结果
print(result[["name", "area_km2", "hospital_count", "hospital_density"]].head())

# 12. 导出结果
output_path = "output/hospital_density.gpkg"
result.to_file(output_path, layer="hospital_density", driver="GPKG")

print("处理完成，结果已导出：", output_path)

常见坑：GeoPandas 实操最容易出错的地方

坑一：直接用 EPSG:4326 计算面积

经纬度坐标系的单位是度，不是米。直接执行 geometry.area 会得到平方度，不能当作平方米或平方公里使用。

正确做法是先用 to_crs() 转换到适合研究区的投影坐标系，再计算面积或距离。

坑二：两个图层 CRS 不一致就做空间连接

空间连接要求参与计算的图层位于同一个坐标参考系统。如果一个是 EPSG:4326，另一个是 EPSG:3857，即使图形在 GIS 软件中看起来能叠上，代码里的空间关系也可能错误。

建议每次分析前都打印：

print(boundary.crs)
print(poi.crs)

坑三：把 define projection 和 project 混为一谈

在 GIS 软件中，“定义投影”是给数据补充坐标系标签，“投影”是把坐标值真正转换到另一个坐标系。GeoPandas 中也是类似逻辑。

• set_crs()：给没有 CRS 的数据指定坐标系，不改变坐标值。
• to_crs()：把数据从当前 CRS 转换到目标 CRS，会改变坐标值。

如果数据本身是 WGS 84，却错误使用 set_crs("EPSG:3857")，后面的所有空间分析都会错。

坑四：空间连接 predicate 选错

gpd.sjoin() 中的 predicate 很关键。常见选项包括：

• within：一个几何完全位于另一个几何内部，常用于点落入面。
• contains：一个几何包含另一个几何，常用于面包含点。
• intersects：两个几何相交，适合线面、面面相交判断。

点统计到行政区时，通常写成“点 within 面”更符合直觉。如果改成“点 contains 面”，结果会为空。

坑五：Shapefile 字段名被截断

Shapefile 字段名长度限制较严格，像 hospital_density 这样的字段可能在导出后被截断。ArcGIS 或 QGIS 打开后字段名变化，会影响后续制图和表达式。

如果不是强制要求 Shapefile，建议用 GeoPackage 作为中间成果格式。

方法比较：ArcGIS、QGIS 和 GeoPandas 该怎么选

工具	适合场景	优势	限制
ArcGIS Pro	标准制图、企业 GIS、复杂地理处理模型	工具完整，制图能力强，商业支持好	授权成本高，批量自动化需要额外脚本能力
QGIS	开源桌面 GIS、日常数据检查和制图	免费开源，插件丰富，支持多格式	复杂自动化流程仍需要模型构建器或 Python
GeoPandas	批量矢量处理、自动化统计、与数据分析流程结合	代码可复现，适合批处理，能与 Pandas、Matplotlib、PostGIS 联动	超大数据性能有限，复杂拓扑修复不如专业 GIS 软件直观
PostGIS	海量空间数据、多人协作、服务端空间查询	空间索引强，适合数据库级分析和 WebGIS 后端	需要数据库部署和 SQL 能力

实际工作中不必把 GeoPandas 和 ArcGIS 对立起来。更合理的方式是：用 ArcGIS 或 QGIS 做数据检查、符号化和制图，用 GeoPandas 做可复现的批量处理和统计分析，用 PostGIS 管理较大的空间数据。

检查清单：运行 GeoPandas 脚本前先确认这些项

• 数据路径是否正确：相对路径从脚本运行目录开始计算，不一定是脚本所在目录。
• 图层是否能正常读取：先用 gpd.read_file() 和 head() 检查。
• 字段名是否一致：示例中的 name、category 要替换为你的真实字段。
• CRS 是否存在：gdf.crs 不能是 None。
• CRS 是否统一：空间连接、裁剪、叠加前，参与图层应在同一 CRS 下。
• 投影是否适合面积计算：不要直接用 EPSG:4326 计算面积。
• 几何是否为空或无效：处理前过滤 geometry.notna() 和 is_empty。
• 空间关系是否选对：点落入面一般用 within，面面相交常用 intersects。
• 输出格式是否合适：中间成果优先 GeoPackage，交付给旧系统时再考虑 Shapefile。
• 结果是否抽样核查：不要只看脚本无报错，要在 QGIS 或 ArcGIS 中打开结果检查空间位置和字段值。

FAQ：GeoPandas官方文档实操常见问题

GeoPandas 能完全替代 ArcGIS 吗？

不能简单说完全替代。GeoPandas 很适合矢量数据批处理、属性统计、空间连接、自动化分析；ArcGIS Pro 在专业制图、企业级工具链、栅格分析、三维和行业扩展方面仍然有优势。更推荐把 GeoPandas 作为 ArcGIS 工作流的自动化补充。

GeoPandas 适合处理多大的数据？

GeoPandas 主要在内存中处理数据，适合中小规模矢量数据。如果数据达到数百万要素、几 GB 以上，可能会变慢或占用大量内存。此时可以考虑 PostGIS、DuckDB 空间扩展、分块处理，或先做空间索引和字段裁剪。

为什么我的 sjoin 结果为空？

最常见原因有三个：第一，两个图层 CRS 不一致；第二，predicate 选错，例如把 within 写成不合适的关系；第三，数据实际空间位置不重叠。建议先在 QGIS 或 ArcGIS 中叠加查看，再打印两个图层的 total_bounds 检查范围。

print(boundary_proj.total_bounds)
print(poi_proj.total_bounds)

为什么 GeoPandas 计算面积和 ArcGIS 不一样？

通常是投影坐标系不同、面积算法不同或几何修复状态不同导致的。请确认两边使用相同坐标系，并且都在投影坐标系下计算面积。如果是跨大范围区域，应选择等面积投影，而不是随意使用 Web Mercator。

GeoPandas 读取中文 Shapefile 乱码怎么办？

可以在读取时尝试指定编码：

gdf = gpd.read_file("data/sample.shp", encoding="utf-8")

如果不行，再尝试 gbk 或 gb18030。更稳妥的办法是把数据转换为 GeoPackage，减少 Shapefile 编码和字段名限制带来的问题。

GeoPandas 输出结果能在 ArcGIS Pro 中打开吗？

可以。GeoPandas 导出的 GeoPackage、Shapefile、GeoJSON 通常都可以在 ArcGIS Pro 或 QGIS 中打开。若用于正式交付，建议导出后在目标软件中检查字段名、坐标系、中文编码和几何类型。

结论：把 GIS 工具箱流程变成可复现代码

GeoPandas 的价值不在于“取代所有桌面 GIS 软件”，而在于把重复的 GIS 数据处理流程变成可复现、可检查、可批量运行的 Python 代码。

通过本文的 GeoPandas 官方文档实操，你已经完成了从读取数据、检查 CRS、属性筛选、空间连接、面积计算到导出成果的一整套流程。对于熟悉 ArcGIS 的读者来说，可以把它理解为用 Python 重写了一遍常见矢量工具箱。

如果你刚开始学习，建议先从一个小数据集练习：一个面图层、一个点图层、一个明确统计目标。等这条主线跑通后，再逐步加入裁剪、缓冲区、叠加分析、可视化和 PostGIS 数据库连接。这样学习 GeoPandas 会比单纯浏览 API 文档更高效，也更接近真实 GIS 项目。