GeoPandas实战：Shapefile文件解析与坐标系核验完整指南

简介：四川省地表水水质国控断面坐标数据包含93个断面，覆盖省内主要河流与流域，以GIS矢量文件形式提供，面向环境监测、水资源管理及地理信息分析人员，可用于断面精确定位、水质监测网络可视化与区域对比研究。压缩包共8个文件，以.shp文件保存断面点位置，.dbf记录省、市、流域、河流、断面名称等属性，.prj定义坐标投影，.sbn/.sbx与.shx保障空间索引效率，.cpg及.xml补充编码与元数据，整体仅7KB，可直接导入ArcGIS、QGIS等平台使用。目前已有995人学习下载，适合环保从业者、科研人员及高校学生用于流域水环境分析、断面分布制图或作为环境数据底图。利用属性表中的经纬度与流域字段，可快速生成断面分布图，并叠加污染源、水文等图层开展综合分析，为水质变化趋势研判、重点区域识别和环保决策提供基础数据支持。

1. 93个国控断面坐标，先拆开再看点

四川省地表水水质国控断面坐标数据，压缩包里 9 个文件装的就是 93 个断面位置，每个断面对应一个点要素，属性表带省、市、流域、河流和断面名称。反直觉的地方在隐性信息：同一经纬度按 CGCS2000 与 WGS84 两种基准读，叠加底图差出几十米；中文断面名编码读错，表格里就是乱码。真正的工作量不在打开 .shp，而在认出 .prj 的坐标系、用对 .cpg 的编码、弄清 .sbn 和 .sbx 要不要一起备份。下面从九个文件职责拆起，用 GeoPandas 完成读取、清洗、出图和缓冲区分析，最后落在坐标系核验的具体做法。适合做环境信息化数据接入的 GIS 工程师，也适合拿这批点做水质评价的业务人员。

2. 拆开 .rar：Shapefile 九个文件的职责与读取顺序

2.1 .shp 主文件：93 个点要素的二进制存储结构

.shp 是几何主体文件，结构分两层：100 字节文件头，后面是若干条记录。每条记录的头部 8 字节保存记录编号和内容长度，主体保存 Shape Type 与坐标序列。这批断面数据属于 Point 类型，每个记录的实体部分就是一组双精度 (X, Y) 坐标。文件头第 32 字节起的 4 字节小端整数表示整个文件的几何类型，值为 1 是纯点文件，值为 11 表示带 Z 值的点要素，值为 21 表示带 M 测度值的点要素。如果这里读到 11 而业务上只要平面坐标，下游提取 geometry.z 会多出一维，数据库入库时容易多出一个字段。

理解这个结构对排错很有用。比如某个国产 GIS 平台打开只有 .shp 没有 .shx 的文件时，读取端要顺序扫描全部记录才能确定几何范围，如果文件里有损坏记录，往往就在这个阶段抛异常。反过来，.shx 存在时可以直接定位每条记录的偏移量，中间记录损坏也能跳过继续读。所以备份和传输时不要把 .shx 丢掉。也可以用一段小代码直接读出文件头确认类型：

import struct with open("四川省.shp", "rb") as f: header = f.read(100) shape_type = struct.unpack("<i", header[32:36])[0] print("Shape 类型值:", shape_type) # 1 表示 Point

这段代码把文件头 100 字节读出来，再按小端序解析第 32 到 35 字节的整数。shape_type 是文件级几何类型声明，读出来是 1 和预期一致；读出 3 或 5 说明交付的不是纯点数据，后续所有按点要素设计的分段分析都要重新评估。读出 11 则要特别注意第三维坐标的来源。

2.2 .dbf 属性表与 .cpg 编码：中文断面名能不能显示

.dbf 是 dBASE 格式属性表，字段顺序决定表格里看到的列。国控断面数据的典型字段设计是省、市、流域、河流、断面名称、经度、纬度，部分版本带断面代码。要注意一个限制：dBASE 字段名最长 10 个字符，中文表头导出时经常被压成拼音缩写，接手后要先建立字段映射再开始分析。比如断面名称可能存成 SECNAME 或 DM 这类短名，不做映射直接写查询语句很容易取错列。

.cpg 文件只有一行内容，声明字符编码。它的作用是给读取端一个提示，但 GeoPandas 的 encoding 参数优先级高于 .cpg。也就是说，即使 .cpg 写的是 UTF-8，只要 read_file 显式传了 encoding="gbk"，就会按 GBK 解析。乱码的典型表现是断面名称显示为问号或一组不可读符号，UTF-8 数据按 GBK 读会出双字符乱码，GBK 数据按 UTF-8 读常见问号。我一般准备两种编码各读一次，对比断面名称是否正常，一次就能锁定正确编码。老版本 ArcMap 导出的图层默认多为 GBK，新版 ArcGIS Pro 或 QGIS 导出多为 UTF-8，这是经验判断的起点。

2.3 .prj 坐标系字符串：CGCS2000 与 WGS84 怎么分辨

.prj 是 WKT 格式的坐标系描述。打开后看到 GEOGCS["China Geodetic Coordinate System 2000"]，对应 CGCS2000 地理坐标系，EPSG 编码 4490；看到 GCS_WGS_1984，对应 EPSG 4326。两套地心基准定义差异在厘米级，日常全景展示看不出问题，但叠加卫星影像或精确边界时，山区点位偏移可达几十米，必须显式统一基准。

常见的操作误区是看到经纬度就默认 WGS84。水文国土领域这些年统一要求 CGCS2000，所以这批断面数据非常可能是 4490。如果 .prj 里不是地理坐标而是投影坐标系，四川地区最常见的是 CGCS2000 3 度带高斯克吕格投影，中央经线按 102E、105E、108E 分段。判断方法很简单：地理坐标的经度值在 97 到 109 之间、纬度值在 26 到 35 之间；投影坐标的 X 是 7 位以上、Y 是 3 到 4 位。两者混用会让断面点全部堆到坐标系原点附近，出图时一眼就能看出来。

2.4 .shx、.sbn、.sbx 与 .shp.xml：索引与元数据文件的价值

.shx 是几何索引，保存每条记录在 .shp 中的字节偏移，fiona 驱动在缺少 .shx 时会顺序扫描，93 个点不算慢，但换成几千个多边形就能感到差异。.sbn 和 .sbx 是 ESRI 专有空间索引对，只在 ArcGIS 桌面产品里参与空间查询加速，QGIS、GeoPandas、PostGIS 都不读这两个文件，丢了不影响基本使用。.shp.xml 是元数据文件，记录数据来源、更新时间和投影描述，环境执法与归档场景里这份文件的价值高于几何本身，建议长期保留。

下表汇总九个文件的角色和缺失后果，便于交接时一次性说清：

提示：9 个文件保持同名前缀并在同一目录下，任何 GIS 读取库都按前缀组合找配套文件，改了一个文件名整套数据就可能识别失败。解压时单独建目录，不要把文件散落在下载文件夹根目录。

3. GeoPandas 读取与清洗：分析之前先做四项检查

3.1 环境准备：用 conda 建独立 GIS 环境

处理点状断面数据，我习惯直接用 GeoPandas，read_file 一次性把几何、属性和坐标系读完，比用 pandas 拆字段再拼坐标高效得多。环境搭建用 conda，让 conda-forge 把 GDAL、GEOS、PROJ 这些底层库的版本对齐，避免手动编译：

conda create -n water python=3.11 -y conda activate water conda install geopandas matplotlib pyproj -c conda-forge -y

参数说明：python=3.11 是当前多数 GIS 依赖已完成适配的版本；geopandas 负责矢量读写，matplotlib 出图，pyproj 处理坐标系转换和 EPSG 查询。单位内网无法访问 conda 源时，pip install geopandas 也可行，但要求系统已经存在可用的 GDAL 动态库，否则编译阶段直接失败，这是最常见的一个环境坑。

3.2 读入数据并核对四个基础指标

import geopandas as gpd gdf = gpd.read_file("四川省/四川省.shp", encoding="utf-8") print("记录数:", len(gdf)) print("坐标系:", gdf.crs) print("几何类型:", gdf.geometry.geom_type.unique()) print("字段列表:", list(gdf.columns)) print(gdf.head(3))

四个输出对应四个检查点。记录数应为 93，多了说明存在重复，少了说明导出丢要素；crs 应显示 EPSG 相关字符串，为 None 就是 .prj 没被读到；几何类型应全部是 Point；字段列表要和源资料的字段说明比对，锁定断面名称和经纬度对应的真实列名。head(3) 抽查前三条，确认中文渲染正常。这套数据如果是从老 ArcMap 导出的，read_file 的 encoding 要改成 gbk，否则断面名称大概率乱码。

下面是字段映射参考，不同版本命名有差异，以实际 columns 输出为准：

需要区分 set_crs 与 to_crs 的用法：set_crs 只是给数据打上坐标系标签，不改变任何坐标数值；to_crs 才是真正的坐标变换。如果 .prj 缺失但通过资料确认原始基准是 CGCS2000，用 gdf.set_crs(epsg=4490) 补标签；如果数据在 WGS84 下采集而 .prj 错标成 CGCS2000，就要用 to_crs 转换。两者用错，点位偏移方向完全不同。

3.3 几何坐标与属性坐标的一致性校验

这类数据有个隐蔽问题：.dbf 里的经纬度字段与 .shp 里的几何坐标不一定同源。有的版本是先从几何导出坐标，再人工把表格贴进属性表，两个来源之间可能存在系统性偏差。校验做法是把两边坐标都取出来算距离：

gdf["geom_x"] = gdf.geometry.x gdf["geom_y"] = gdf.geometry.y gdf["xy_diff"] = ( (gdf["geom_x"] - gdf["经度"]) ** 2 + (gdf["geom_y"] - gdf["纬度"]) ** 2 ) ** 0.5 print(gdf.loc[gdf["xy_diff"].idxmax(), ["断面名称", "经度", "纬度", "geom_x", "geom_y", "xy_diff"]])

逻辑说明：先利用 point 几何对象的 x 和 y 属性提取几何坐标，再与属性坐标做欧氏距离。经纬度数据里 0.001 度大约对应 110 米，如果某个断面的 xy_diff 超过 0.002 度，基本可以判定属性列是后期补录的。分析时以几何坐标为准，属性坐标只用于与外部系统比对；偏差过大的记录单独建清单，返回数据提供方确认，不自行修改原始值。

3.4 缺失值与重复断面的处理

print(gdf.isna().sum()) print("重复断面数:", gdf["断面名称"].duplicated().sum())

第一行按列统计缺失值，重点关注断面名称和河流。断面名称是国控点编码的基础，缺失会直接导致后续按断面关联水质监测数据失败；河流字段缺失影响按流域口径的统计。第二行统计重复断面名，理论上国控断面编号不会重复。出现重复时先对比两行几何坐标，坐标一致说明导出冗余，保留任意一条；坐标不同则可能是同名取水点和监测断面混入，不能自行取舍，要交回业务方核实。

注意：清洗过程不修改归档原件。先复制一份原始解压目录，所有写操作放在副本上，改坏还能退回，这是环境数据项目里最便宜的容灾手段。

清洗完成的 GeoDataFrame 按惯例导出为 GeoPackage 归档，替代再生成一份 Shapefile：

gdf_clean = gdf.drop(columns=["geom_x", "geom_y", "xy_diff"]) gdf_clean.to_file("sichuan_sections_clean.gpkg", layer="sichuan_sections", driver="GPKG")

drop 把校验用的中间列去掉，避免污染正式属性表。GPKG 是 SQLite 容器的矢量格式，单文件归档、支持空间索引、字段名不受 10 字符限制，作为项目交付格式比 Shapefile 稳妥，后续第 4 章的分析直接读这个 gpkg 即可。

4. 断面点位出图与空间关联：从 93 个点到分析结论

4.1 基础分布图：一张图找出八成数据问题

环境信息化的第一屏基本是空间分布。用 matplotlib 出静态图，既能快速发现问题，又保留坐标轴与样式的完全控制：

import matplotlib.pyplot as plt fig, ax = plt.subplots(figsize=(10, 8)) gdf.plot(ax=ax, markersize=15, color="#d62728", alpha=0.8) ax.set_title("四川省地表水水质国控断面分布") ax.set_xlabel("经度") ax.set_ylabel("纬度") ax.grid(True, linestyle="--", alpha=0.4) plt.savefig("sections_distribution.png", dpi=300) plt.show()

gdf.plot 直接用几何画点，不用手动拆坐标；markersize 控制点面积，15 在 93 个点的规模下清晰且不互相遮盖；alpha 调透明度，重叠点时能看到密度；网格线方便比对经纬度刻度。出图后重点观察断面是否沿主要水系呈线状分布，以及川西高原与成都平原边界处有没有孤立的异常点。个别点落在明显偏离水系的山区，优先回到 3.3 的校验结果里比对该断面的几何坐标与属性坐标。

需要按流域区分断面时，column 参数配合 legend 出分类图：

fig, ax = plt.subplots(figsize=(10, 8)) gdf.plot(ax=ax, column="流域", cmap="tab20", legend=True, markersize=15) ax.set_title("按流域着色的国控断面分布") plt.tight_layout() plt.show()

column 指定分类字段，cmap 用 tab20 保证二十个类别以内的色彩可区分，legend=True 自动在图侧生成图例。如果流域字段存在空值，GeoPandas 会给空值单独一种颜色，这正好帮我们发现属性缺失的断面。保存报告用图时 dpi 至少 300，矢量格式优先考虑 SVG，后续改字大小不损失清晰度。

4.2 叠合河网：sjoin_nearest 核对断面河流归属

断面本身的合理性要靠河网交叉验证。拿到四川省河网线图层后，常见做法是把最近的河流线段匹配到每个断面：

rivers = gpd.read_file("sichuan_rivers.shp") gdf_web = gdf.to_crs(epsg=3857) rivers_web = rivers.to_crs(epsg=3857) check = gpd.sjoin_nearest( gdf_web, rivers_web[["河流名称", "geometry"]], how="left", max_distance=5000, ) print(check[["断面名称", "河流名称"]])

两个关键参数要理解。to_crs(epsg=3857) 把数据转到 Web Mercator 投影，距离单位从此是米；如果不投影，max_distance=5000 会被当成 5000 度，几乎匹配不到任何线段。max_distance 是搜索半径，国控断面就在河道上，5 公里已经足够宽松。匹配结果为空白的断面，要么坐标偏移，要么河网数据缺该段河道，这两类记录都应列入疑点清单。

4.3 缓冲区与空间连接：周边污染源怎么统计

水质分析里常需要回答每个国控断面周边有多少重点排污单位。做法是先投影再缓冲，再把企业点位和缓冲区做空间连接：

buf = gdf.to_crs(epsg=3857).copy() buf["geometry"] = buf.geometry.buffer(3000) enterprises = gpd.read_file("key_enterprises.shp").to_crs(epsg=3857) joined = gpd.sjoin( enterprises, buf[["断面名称", "geometry"]], how="inner", predicate="within", ) rank = joined.groupby("断面名称").size().reset_index(name="周边企业数") print(rank.sort_values("周边企业数", ascending=False))

两个翻车点。第一，buffer(3000) 必须在投影坐标系下执行，经纬度坐标里 3000 被解释成 3000 度，生成的圆覆盖整个区域。第二，predicate="within" 的方向是判断企业点是否落在缓冲区面内，写反成面在点内结果恒为空。groupby 之后按断面分组计数，得到每个断面 3 公里范围内的企业数量，降序排列后可直接作为现场核查与优先排序的依据。缓冲区距离不是固定值，按业务口径调整：

更精确的做法是用河网线做上游汇水区，而不是简单圆缓冲。圆缓冲包含下游区域，实际污染影响只在断面上游，严谨的评价体系里要用流向数据裁剪。这里给出的 3 公里圆缓冲是初筛口径，正式结论需要水文分析支撑。

5. 坐标系核验与投影转换的落地技巧

5.1 用 pyproj 从 .prj 反查 EPSG

手上没有 GIS 桌面软件时，直接用 pyproj 解析 .prj 是最快的确认方式：

from pyproj import CRS wkt = open("四川省.prj", encoding="utf-8").read() crs = CRS.from_wkt(wkt) print(crs.to_epsg()) # 预期 4490 或 4326

CRS.from_wkt 解析 WKT 定义，to_epsg 返回匹配的 EPSG 编码。4490 对应 CGCS2000 地理坐标系，4326 对应 WGS84。返回 None 时说明 .prj 里是自定义或扩展参数，需要把 crs.to_wkt() 完整输出与源数据说明逐项比对，重点看 datum 和 spheroid 两段。

5.2 统一基准之后再算距离

多源数据叠用场景里，断面来自国土部门常用 CGCS2000，污染源点位来自企业申报常用 WGS84，直接叠加前先统一到同一基准：

gdf_4326 = gdf.to_crs(epsg=4326) gdf_4490 = gdf.to_crs(epsg=4490) delta = ( gdf_4326.to_crs(epsg=3857).geometry .distance(gdf_4490.to_crs(epsg=3857).geometry) ) print("基准差异等效距离:", delta.max())

这段验证代码说明一个事实：CGCS2000 与 WGS84 的框架差异换算成距离在厘米级，远小于 GIS 出图的像素误差。真正要防的是把经纬度坐标直接当平面距离相减，在四川地区 1 度纬向长度约 111 公里，未经投影的距离计算结果没有物理意义。所有缓冲区、最近邻、断面间距计算，统一先 to_crs 再运算。

5.3 用已知断面反查做最终验收

所有处理结束后，抽 3 到 5 个有公开位置的断面做最终核验。把经纬度输入在线地图反查，确认点位是否落在河道主流上。比如岷江、沱江干流上的国控断面如果反查落在山腰或建成区，说明原始坐标存在基准错配或投影错误，必须退回原始 .shp 重新检查几何，而不是修改属性列掩盖问题。这一步工程量很小，却是坐标系处理链路里唯一能证明数据真正可用的验证环节。建议把核验结果截图存档，与断面点位图、属性表一起作为数据验收附件，后续再有人质疑这批断面坐标，直接用核验记录回应即可。

本文还有配套的精品资源，点击获取