📚 课程目录第2章 GeoPandas环境搭建与地质数据准备

2. GeoPandas环境搭建与地质数据准备

批处理地质矢量数据的第一步不是写循环，而是把运行环境和数据底座搭稳。GeoPandas 本身只是“表格 + 几何”的封装层，它依赖 Shapely 处理几何对象、pyproj 处理坐标参考系、pyogrio 或 Fiona 读写矢量文件、matplotlib 出图。只要其中任一依赖版本错配，后续批量读取就可能出现几何类型识别错误、CRS 丢失或驱动缺失。因此本节按“环境 → 数据 → 目录 → 批量读取”的顺序展开，先保证单文件能稳定读入，再谈批量化。

一、用 conda 建立可复现的环境

GeoPandas 涉及大量带 C 扩展的二进制包，直接用 pip 在系统 Python 上安装容易与已有的 GDAL 冲突。推荐用 conda（或 mamba）创建独立环境，让 GDAL、PROJ、GEOS 由同一渠道统一解析版本：

conda create -n geobatch python=3.11 -y
conda activate geobatch
conda install -c conda-forge geopandas pyogrio shapely pyproj fiona matplotlib pandas -y

安装完成后必须做一次自检，确认版本与驱动可用，而不是只看 import 是否报错：

import geopandas as gpd, pyproj, shapely
print(gpd.__version__, shapely.__version__, pyproj.__version__)
print(gpd.list_layers("data/raw/geo_map.shp"))   # 查看图层名与几何类型
print(gpd.io.file.fiona.drvsupport.supported_drivers.keys())

若 list_layers 能返回图层信息，说明矢量驱动链路正常。GeoPandas 0.14 之后默认使用 pyogrio 引擎，读取速度明显优于 Fiona；如遇个别格式兼容问题，可在读取时显式指定 engine="fiona" 回退。建议把环境导出为 environment.yml 随项目一起保存，保证换机器后能重建同一套依赖。

二、准备典型地质矢量数据

地质批处理的对象通常包含四类图层，它们的几何类型和属性结构差异很大，正好用来检验批处理脚本的健壮性：

地质图面图层：Polygon 或 MultiPolygon，属性含地层代号（如 stratum_code）、岩性名称、时代代号，是面积统计与专题制图的主体。
断层线图层：LineString，属性含断层性质（正断层/逆断层/平移断层）、倾向、倾角、长度，用于构造分析和缓冲区计算。
地层界线图层：LineString，表示不同地层单元的分界，常与面图层共享 stratum_code 以便拓扑核对。
钻孔点图层：Point，属性含孔号、孔口高程、孔深、终孔层位，用于与面图层做空间连接取地层属性。

教学示例中可假设数据来自某 1:5 万图幅，统一采用 CGCS2000 地理坐标（EPSG:4490）存放原始数据。需要强调的是，原始数据应保持“只读”状态，任何裁剪、投影转换、字段重命名都写入派生目录，避免批处理中途污染源数据而无法回溯。

三、数据目录组织

批处理脚本最容易失控的地方是路径。建议按“原始 / 中间 / 成果 / 配置”四层组织，并让脚本只通过配置读取路径，而不是在代码里硬编码：

geology_project/
├─ data/
│  ├─ raw/          # 原始矢量，只读
│  │  ├─ geo_map.shp
│  │  ├─ faults.shp
│  │  ├─ boundaries.shp
│  │  └─ boreholes.shp
│  ├─ interim/      # 投影转换、裁剪后的中间结果
│  └─ output/       # 统计表、成图、合并后的 GeoPackage
├─ config/
│  └─ layers.yaml   # 图层名、目标 CRS、关键字段
└─ scripts/
   └─ batch_read.py

用一份 YAML 描述每个图层的元信息，脚本据此决定读取后做什么处理，这比在代码里写 if-else 判断文件名更易维护：

target_crs: EPSG:4547        # CGCS2000 / 3-degree Gauss-Kruger CM 117E
layers:
  - name: geo_map
    path: data/raw/geo_map.shp
    geom_type: polygon
    key_field: stratum_code
  - name: faults
    path: data/raw/faults.shp
    geom_type: line
    key_field: fault_type
  - name: boreholes
    path: data/raw/boreholes.shp
    geom_type: point
    key_field: hole_id

四、批量读取策略

批量读取的核心是把“逐个文件读入”变成“统一入口 + 统一校验”。下面这段示例按配置遍历图层，读入后立即检查几何类型、CRS 和空几何，并统一转换到目标投影：

import yaml, geopandas as gpd
from pathlib import Path

cfg = yaml.safe_load(open("config/layers.yaml", encoding="utf-8"))
target = cfg["target_crs"]
store = {}

for item in cfg["layers"]:
    gdf = gpd.read_file(item["path"])
    if gdf.crs is None:
        raise ValueError(f'{item["name"]} 缺少 CRS，需先确认坐标系')
    if gdf.crs.to_string() != target:
        gdf = gdf.to_crs(target)
    gdf = gdf[~gdf.geometry.is_empty & gdf.geometry.notna()]
    store[item["name"]] = gdf
    print(item["name"], len(gdf), gdf.geom_type.unique())

运行后应观察三类信号：记录数是否与源文件一致（差异通常来自空几何被剔除）、geom_type 是否与配置声明相符（面图层混入 LineString 往往意味着源数据有误）、CRS 是否已统一为 EPSG:4547。只有这三项都通过，才进入后续的叠加、缓冲与统计。对于成果输出，推荐统一写入 GeoPackage 而非 Shapefile，因为后者字段名限 10 字符、不支持中文编码且单文件容量受限，而 GeoPackage 可在一个文件中存放多图层并保留完整字段。

完成上述准备后，判断环境与数据是否就绪的标准很具体：四个图层都能被 read_file 读入且记录数非零；几何类型与配置声明一致；CRS 全部统一到 EPSG:4547；空几何已被剔除；成果能成功写入 GeoPackage 并被重新读回。满足这些条件，后续的批量裁剪、空间连接与面积统计才有可靠输入。