简介四川省地表水水质国控断面坐标数据包含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. 拆开 .rarShapefile 九个文件的职责与读取顺序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 显式传了 encodinggbk就会按 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 是元数据文件记录数据来源、更新时间和投影描述环境执法与归档场景里这份文件的价值高于几何本身建议长期保留。下表汇总九个文件的角色和缺失后果便于交接时一次性说清文件作用缺失后果.shp几何主体无法读取.shx记录偏移索引读取变慢可重建.dbf属性表断面名称、河流等信息丢失.prj坐标系定义crs 为 None叠加错位.cpg编码声明中文乱码风险.sbnESRI 空间索引仅影响 ArcGIS 部分查询.sbxESRI 空间索引附属同上.shp.xml元数据数据来源与更新时间不可考.rar压缩容器不影响使用需先解压提示9 个文件保持同名前缀并在同一目录下任何 GIS 读取库都按前缀组合找配套文件改了一个文件名整套数据就可能识别失败。解压时单独建目录不要把文件散落在下载文件夹根目录。3. GeoPandas 读取与清洗分析之前先做四项检查3.1 环境准备用 conda 建独立 GIS 环境处理点状断面数据我习惯直接用 GeoPandasread_file 一次性把几何、属性和坐标系读完比用 pandas 拆字段再拼坐标高效得多。环境搭建用 conda让 conda-forge 把 GDAL、GEOS、PROJ 这些底层库的版本对齐避免手动编译conda create -n water python3.11 -y conda activate water conda install geopandas matplotlib pyproj -c conda-forge -y参数说明python3.11 是当前多数 GIS 依赖已完成适配的版本geopandas 负责矢量读写matplotlib 出图pyproj 处理坐标系转换和 EPSG 查询。单位内网无法访问 conda 源时pip install geopandas 也可行但要求系统已经存在可用的 GDAL 动态库否则编译阶段直接失败这是最常见的一个环境坑。3.2 读入数据并核对四个基础指标import geopandas as gpd gdf gpd.read_file(四川省/四川省.shp, encodingutf-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 输出为准资料页字段常见存储列名用途断面名称SECNAME 或 DM国控点唯一标识河流RIVER 或 HL空间归属核验流域BASIN 或 LY业务分类汇总经度LON 或 JD与几何坐标比对纬度LAT 或 WD与几何坐标比对需要区分 set_crs 与 to_crs 的用法set_crs 只是给数据打上坐标系标签不改变任何坐标数值to_crs 才是真正的坐标变换。如果 .prj 缺失但通过资料确认原始基准是 CGCS2000用 gdf.set_crs(epsg4490) 补标签如果数据在 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 归档替代再生成一份 Shapefilegdf_clean gdf.drop(columns[geom_x, geom_y, xy_diff]) gdf_clean.to_file(sichuan_sections_clean.gpkg, layersichuan_sections, driverGPKG)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(axax, markersize15, color#d62728, alpha0.8) ax.set_title(四川省地表水水质国控断面分布) ax.set_xlabel(经度) ax.set_ylabel(纬度) ax.grid(True, linestyle--, alpha0.4) plt.savefig(sections_distribution.png, dpi300) plt.show()gdf.plot 直接用几何画点不用手动拆坐标markersize 控制点面积15 在 93 个点的规模下清晰且不互相遮盖alpha 调透明度重叠点时能看到密度网格线方便比对经纬度刻度。出图后重点观察断面是否沿主要水系呈线状分布以及川西高原与成都平原边界处有没有孤立的异常点。个别点落在明显偏离水系的山区优先回到 3.3 的校验结果里比对该断面的几何坐标与属性坐标。需要按流域区分断面时column 参数配合 legend 出分类图fig, ax plt.subplots(figsize(10, 8)) gdf.plot(axax, column流域, cmaptab20, legendTrue, markersize15) ax.set_title(按流域着色的国控断面分布) plt.tight_layout() plt.show()column 指定分类字段cmap 用 tab20 保证二十个类别以内的色彩可区分legendTrue 自动在图侧生成图例。如果流域字段存在空值GeoPandas 会给空值单独一种颜色这正好帮我们发现属性缺失的断面。保存报告用图时 dpi 至少 300矢量格式优先考虑 SVG后续改字大小不损失清晰度。4.2 叠合河网sjoin_nearest 核对断面河流归属断面本身的合理性要靠河网交叉验证。拿到四川省河网线图层后常见做法是把最近的河流线段匹配到每个断面rivers gpd.read_file(sichuan_rivers.shp) gdf_web gdf.to_crs(epsg3857) rivers_web rivers.to_crs(epsg3857) check gpd.sjoin_nearest( gdf_web, rivers_web[[河流名称, geometry]], howleft, max_distance5000, ) print(check[[断面名称, 河流名称]])两个关键参数要理解。to_crs(epsg3857) 把数据转到 Web Mercator 投影距离单位从此是米如果不投影max_distance5000 会被当成 5000 度几乎匹配不到任何线段。max_distance 是搜索半径国控断面就在河道上5 公里已经足够宽松。匹配结果为空白的断面要么坐标偏移要么河网数据缺该段河道这两类记录都应列入疑点清单。4.3 缓冲区与空间连接周边污染源怎么统计水质分析里常需要回答每个国控断面周边有多少重点排污单位。做法是先投影再缓冲再把企业点位和缓冲区做空间连接buf gdf.to_crs(epsg3857).copy() buf[geometry] buf.geometry.buffer(3000) enterprises gpd.read_file(key_enterprises.shp).to_crs(epsg3857) joined gpd.sjoin( enterprises, buf[[断面名称, geometry]], howinner, predicatewithin, ) rank joined.groupby(断面名称).size().reset_index(name周边企业数) print(rank.sort_values(周边企业数, ascendingFalse))两个翻车点。第一buffer(3000) 必须在投影坐标系下执行经纬度坐标里 3000 被解释成 3000 度生成的圆覆盖整个区域。第二predicatewithin 的方向是判断企业点是否落在缓冲区面内写反成面在点内结果恒为空。groupby 之后按断面分组计数得到每个断面 3 公里范围内的企业数量降序排列后可直接作为现场核查与优先排序的依据。缓冲区距离不是固定值按业务口径调整分析目标缓冲距离投影要求周边污染源初查3000 米必须投影饮用水源风险排查5000 米必须投影省界跨界断面核查按行政边界裁剪使用面积保留投影更精确的做法是用河网线做上游汇水区而不是简单圆缓冲。圆缓冲包含下游区域实际污染影响只在断面上游严谨的评价体系里要用流向数据裁剪。这里给出的 3 公里圆缓冲是初筛口径正式结论需要水文分析支撑。5. 坐标系核验与投影转换的落地技巧5.1 用 pyproj 从 .prj 反查 EPSG手上没有 GIS 桌面软件时直接用 pyproj 解析 .prj 是最快的确认方式from pyproj import CRS wkt open(四川省.prj, encodingutf-8).read() crs CRS.from_wkt(wkt) print(crs.to_epsg()) # 预期 4490 或 4326CRS.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(epsg4326) gdf_4490 gdf.to_crs(epsg4490) delta ( gdf_4326.to_crs(epsg3857).geometry .distance(gdf_4490.to_crs(epsg3857).geometry) ) print(基准差异等效距离:, delta.max())这段验证代码说明一个事实CGCS2000 与 WGS84 的框架差异换算成距离在厘米级远小于 GIS 出图的像素误差。真正要防的是把经纬度坐标直接当平面距离相减在四川地区 1 度纬向长度约 111 公里未经投影的距离计算结果没有物理意义。所有缓冲区、最近邻、断面间距计算统一先 to_crs 再运算。5.3 用已知断面反查做最终验收所有处理结束后抽 3 到 5 个有公开位置的断面做最终核验。把经纬度输入在线地图反查确认点位是否落在河道主流上。比如岷江、沱江干流上的国控断面如果反查落在山腰或建成区说明原始坐标存在基准错配或投影错误必须退回原始 .shp 重新检查几何而不是修改属性列掩盖问题。这一步工程量很小却是坐标系处理链路里唯一能证明数据真正可用的验证环节。建议把核验结果截图存档与断面点位图、属性表一起作为数据验收附件后续再有人质疑这批断面坐标直接用核验记录回应即可。本文还有配套的精品资源点击获取
简介四川省地表水水质国控断面坐标数据包含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. 拆开 .rarShapefile 九个文件的职责与读取顺序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 显式传了 encodinggbk就会按 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 是元数据文件记录数据来源、更新时间和投影描述环境执法与归档场景里这份文件的价值高于几何本身建议长期保留。下表汇总九个文件的角色和缺失后果便于交接时一次性说清文件作用缺失后果.shp几何主体无法读取.shx记录偏移索引读取变慢可重建.dbf属性表断面名称、河流等信息丢失.prj坐标系定义crs 为 None叠加错位.cpg编码声明中文乱码风险.sbnESRI 空间索引仅影响 ArcGIS 部分查询.sbxESRI 空间索引附属同上.shp.xml元数据数据来源与更新时间不可考.rar压缩容器不影响使用需先解压提示9 个文件保持同名前缀并在同一目录下任何 GIS 读取库都按前缀组合找配套文件改了一个文件名整套数据就可能识别失败。解压时单独建目录不要把文件散落在下载文件夹根目录。3. GeoPandas 读取与清洗分析之前先做四项检查3.1 环境准备用 conda 建独立 GIS 环境处理点状断面数据我习惯直接用 GeoPandasread_file 一次性把几何、属性和坐标系读完比用 pandas 拆字段再拼坐标高效得多。环境搭建用 conda让 conda-forge 把 GDAL、GEOS、PROJ 这些底层库的版本对齐避免手动编译conda create -n water python3.11 -y conda activate water conda install geopandas matplotlib pyproj -c conda-forge -y参数说明python3.11 是当前多数 GIS 依赖已完成适配的版本geopandas 负责矢量读写matplotlib 出图pyproj 处理坐标系转换和 EPSG 查询。单位内网无法访问 conda 源时pip install geopandas 也可行但要求系统已经存在可用的 GDAL 动态库否则编译阶段直接失败这是最常见的一个环境坑。3.2 读入数据并核对四个基础指标import geopandas as gpd gdf gpd.read_file(四川省/四川省.shp, encodingutf-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 输出为准资料页字段常见存储列名用途断面名称SECNAME 或 DM国控点唯一标识河流RIVER 或 HL空间归属核验流域BASIN 或 LY业务分类汇总经度LON 或 JD与几何坐标比对纬度LAT 或 WD与几何坐标比对需要区分 set_crs 与 to_crs 的用法set_crs 只是给数据打上坐标系标签不改变任何坐标数值to_crs 才是真正的坐标变换。如果 .prj 缺失但通过资料确认原始基准是 CGCS2000用 gdf.set_crs(epsg4490) 补标签如果数据在 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 归档替代再生成一份 Shapefilegdf_clean gdf.drop(columns[geom_x, geom_y, xy_diff]) gdf_clean.to_file(sichuan_sections_clean.gpkg, layersichuan_sections, driverGPKG)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(axax, markersize15, color#d62728, alpha0.8) ax.set_title(四川省地表水水质国控断面分布) ax.set_xlabel(经度) ax.set_ylabel(纬度) ax.grid(True, linestyle--, alpha0.4) plt.savefig(sections_distribution.png, dpi300) plt.show()gdf.plot 直接用几何画点不用手动拆坐标markersize 控制点面积15 在 93 个点的规模下清晰且不互相遮盖alpha 调透明度重叠点时能看到密度网格线方便比对经纬度刻度。出图后重点观察断面是否沿主要水系呈线状分布以及川西高原与成都平原边界处有没有孤立的异常点。个别点落在明显偏离水系的山区优先回到 3.3 的校验结果里比对该断面的几何坐标与属性坐标。需要按流域区分断面时column 参数配合 legend 出分类图fig, ax plt.subplots(figsize(10, 8)) gdf.plot(axax, column流域, cmaptab20, legendTrue, markersize15) ax.set_title(按流域着色的国控断面分布) plt.tight_layout() plt.show()column 指定分类字段cmap 用 tab20 保证二十个类别以内的色彩可区分legendTrue 自动在图侧生成图例。如果流域字段存在空值GeoPandas 会给空值单独一种颜色这正好帮我们发现属性缺失的断面。保存报告用图时 dpi 至少 300矢量格式优先考虑 SVG后续改字大小不损失清晰度。4.2 叠合河网sjoin_nearest 核对断面河流归属断面本身的合理性要靠河网交叉验证。拿到四川省河网线图层后常见做法是把最近的河流线段匹配到每个断面rivers gpd.read_file(sichuan_rivers.shp) gdf_web gdf.to_crs(epsg3857) rivers_web rivers.to_crs(epsg3857) check gpd.sjoin_nearest( gdf_web, rivers_web[[河流名称, geometry]], howleft, max_distance5000, ) print(check[[断面名称, 河流名称]])两个关键参数要理解。to_crs(epsg3857) 把数据转到 Web Mercator 投影距离单位从此是米如果不投影max_distance5000 会被当成 5000 度几乎匹配不到任何线段。max_distance 是搜索半径国控断面就在河道上5 公里已经足够宽松。匹配结果为空白的断面要么坐标偏移要么河网数据缺该段河道这两类记录都应列入疑点清单。4.3 缓冲区与空间连接周边污染源怎么统计水质分析里常需要回答每个国控断面周边有多少重点排污单位。做法是先投影再缓冲再把企业点位和缓冲区做空间连接buf gdf.to_crs(epsg3857).copy() buf[geometry] buf.geometry.buffer(3000) enterprises gpd.read_file(key_enterprises.shp).to_crs(epsg3857) joined gpd.sjoin( enterprises, buf[[断面名称, geometry]], howinner, predicatewithin, ) rank joined.groupby(断面名称).size().reset_index(name周边企业数) print(rank.sort_values(周边企业数, ascendingFalse))两个翻车点。第一buffer(3000) 必须在投影坐标系下执行经纬度坐标里 3000 被解释成 3000 度生成的圆覆盖整个区域。第二predicatewithin 的方向是判断企业点是否落在缓冲区面内写反成面在点内结果恒为空。groupby 之后按断面分组计数得到每个断面 3 公里范围内的企业数量降序排列后可直接作为现场核查与优先排序的依据。缓冲区距离不是固定值按业务口径调整分析目标缓冲距离投影要求周边污染源初查3000 米必须投影饮用水源风险排查5000 米必须投影省界跨界断面核查按行政边界裁剪使用面积保留投影更精确的做法是用河网线做上游汇水区而不是简单圆缓冲。圆缓冲包含下游区域实际污染影响只在断面上游严谨的评价体系里要用流向数据裁剪。这里给出的 3 公里圆缓冲是初筛口径正式结论需要水文分析支撑。5. 坐标系核验与投影转换的落地技巧5.1 用 pyproj 从 .prj 反查 EPSG手上没有 GIS 桌面软件时直接用 pyproj 解析 .prj 是最快的确认方式from pyproj import CRS wkt open(四川省.prj, encodingutf-8).read() crs CRS.from_wkt(wkt) print(crs.to_epsg()) # 预期 4490 或 4326CRS.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(epsg4326) gdf_4490 gdf.to_crs(epsg4490) delta ( gdf_4326.to_crs(epsg3857).geometry .distance(gdf_4490.to_crs(epsg3857).geometry) ) print(基准差异等效距离:, delta.max())这段验证代码说明一个事实CGCS2000 与 WGS84 的框架差异换算成距离在厘米级远小于 GIS 出图的像素误差。真正要防的是把经纬度坐标直接当平面距离相减在四川地区 1 度纬向长度约 111 公里未经投影的距离计算结果没有物理意义。所有缓冲区、最近邻、断面间距计算统一先 to_crs 再运算。5.3 用已知断面反查做最终验收所有处理结束后抽 3 到 5 个有公开位置的断面做最终核验。把经纬度输入在线地图反查确认点位是否落在河道主流上。比如岷江、沱江干流上的国控断面如果反查落在山腰或建成区说明原始坐标存在基准错配或投影错误必须退回原始 .shp 重新检查几何而不是修改属性列掩盖问题。这一步工程量很小却是坐标系处理链路里唯一能证明数据真正可用的验证环节。建议把核验结果截图存档与断面点位图、属性表一起作为数据验收附件后续再有人质疑这批断面坐标直接用核验记录回应即可。本文还有配套的精品资源点击获取
