简介这份数据面向GIS分析、城市规划与交通研究从业者提供河北省最新分级道路矢量成果精度可达乡道级别。压缩包共97个文件以Shapefile格式为主包含.shp几何、.dbf属性、.prj投影、.shx索引等配套文件整包约73.95MB解压后可直接加载至ArcGIS或QGIS使用。资源涵盖城市一至四级道路、高速、国道、省道、县道、乡道以及OSM来源的铁路、地铁、轻轨、主干道、住宅街道等多种路网类型共16套道路数据不同来源数据相互补充既可进行专题制图也能用于缓冲区分析、路径规划等场景。目前已有312人学习下载适合需要按行政等级或功能等级分类调用道路数据的项目制作者参考。1. 一份 rar 解开后为什么全是同名文件拿到“河北省道路数据最新分级精确到乡道矢量数据.rar”的时候第一反应多半是里面怎么有几十个“同名”文件比如高速_河北省.shp、高速_河北省.dbf、高速_河北省.shx……这不是重复文件而是一整套 ESRI Shapefile。一个完整要素类本来就该由 .shp、.shx、.dbf、.prj 等一堆文件组成缺一个都可能打不开或者丢投影。这套数据把河北省道路按行政等级分成高速、国道、省道、县道、乡道又按城市内部等级分成一至四级道路还叠加了一份 OSM 来源的 roads 和 railways 线层总共 16 组线要素。对做路网分析、地图制图、城市规划和交通物流的人来说它解决的是“市级道路上哪找、县道乡道有没有、数据能不能直接进数据库”这一类问题。后文我会先从 Shapefile 文件家族说起再给出一套可以用命令行和 Python 直接跑的拆分、比对、出图流程并说明每个环节的参数怎么调、坑在哪。2. Shapefile 文件家族从 .shp 到 .sbn 一个个说清2.1 同一目录下每个扩展名都是干什么的Shapefile 不是单文件而是一个文件家族。这套河北道路数据里出现的扩展名基本把常见的都覆盖全了。先看一张我整理的文件与用途对照表扩展名作用缺失时的影响.shp要素几何本体记录每个线要素的形状文件打不开核心文件.shx形状索引用于快速定位 .shp 中的要素某些软件无法读取显示损坏.dbf属性表存道路名称、等级、类型等字段属性全丢只剩几何.prj投影与坐标系描述文件文本格式软件按 WGS84 猜测坐标可能错位.cpg属性表编码声明例如 UTF-8、GBK中文乱码的风险大幅上升.sbn / .sbx空间索引由 ArcGIS 生成可重建删除后下次打开会自动生成.shp.xml元数据XML 格式不影响使用只丢描述信息注意一个细节很多从网盘或公众号下回来的矢量数据.sbn和.sbx是成对出现的。它俩是 ArcGIS 在创建要素类时附带的空间索引你用 QGIS 或 GDAL 读取时并不会依赖它真正决定能否读取的是 .shp、.shx、.dbf 三个文件。所以如果解压后报缺失先把这三个文件确认齐了再谈其他。2.2 用 ogrinfo 快速验证数据完整性与坐标系我拿到数据的第一步从来不是直接拉进 QGIS而是在终端用ogrinfo看头部信息成本最低且能避开图形界面缓存问题。以其中一份乡道_河北省.shp为例ogrinfo -ro -so -al 乡道_河北省.shp参数说明-ro表示只读打开避免误写-so是 summary only只打印图层概要而不展开每条要素-al表示列出所有图层。输出里重点看这几行Layer name: 乡道_河北省 Geometry: Line String Feature Count: 27194 Extent: (113.608000, 36.051000) - (119.857000, 42.612000) Layer SRS: GEOGCS[WGS 84]如果Feature Count是 0多半是.shx和.dbf不配套如果Layer SRS显示未知说明缺失 .prj 或 .prj 为空。此时我用另一条命令重新指定坐标系ogrinfo -ro -al -so -oo ENCODINGUTF-8 乡道_河北省.shp-oo ENCODING是 open option告诉 GDAL 按 UTF-8 解码 dbf能解决一部分中文乱码。实际项目中我遇到过 .cpg 写的是 UTF-8但内部实际是 GBK 的混乱情况单看扩展名没用要靠输出里属性字段的值是否可读来判断。2.3 用 Python 列出字段和统计要素数ogrinfo 适合快速看但要做批量检查我更喜欢用 Python 的fiona或pyogrio一次遍历整个目录。下面脚本会扫描所有*_河北省.shp打印每个图层的要素数量、几何类型以及 dbf 里有哪些字段import glob import fiona for shp in glob.glob(*_河北省.shp): with fiona.open(shp, encodingUTF-8) as src: schema src.schema props list(schema[properties].keys()) print(f{shp}: {len(src)} 个要素, 类型{schema[geometry]}) print(f 字段: {props[:8]}{... if len(props) 8 else })逻辑说明fiona.open的encoding参数对应 dbf 的读入编码和 GDAL 的ENCODING等价schema[properties]是一个有序字典字段顺序与 dbf 物理顺序一致。这套数据里的“各级道路OSM_河北省.shp”字段较多通常包含highway、name、ref、oneway等 OSM 标准字段而行政等级道路层字段相对精简常见的是name、type、length。先看清字段名后面做拆分筛选才有依据。3. 按行政等级拆路网高速、国道、省道、县道、乡道一次拆干净3.1 确认分级字段的值域忻州、保定这种多等级混合的区域道路层如果已经按文件拆分那还相对省事。但实际解压后你会发现既有高速_河北省.shp又有各级道路OSM_河北省.shp后者是把所有等级混在一个图层里的。所以第一步永远是确认分层文件里的等级字段到底是什么值。用 2.3 节的脚本先看一遍然后对某个混编文件执行ogrinfo -ro -al -so -sql SELECT DISTINCT highway FROM 各级道路OSM_河北省 各级道路OSM_河北省.shp这里如果highway字段不存在会直接报no such column。报错时先跑ogrinfo -ro -al -so 各级道路OSM_河北省.shp把真实字段名打出来常见的是type、fclass、kind。OSM 来源的数据里highway的典型值包括motorway、trunk、primary、secondary、tertiary、residential、track等。行政等级和 OSM 等级是两个体系不能简单画等号。3.2 用 ogr2ogr 按属性拆出单独图层当你确认了分级字段之后拆分就很简单。下面命令把 OSM 全量路网按highwaymotorway拆出来另存为独立 Shpogr2ogr -f ESRI Shapefile 河北高速_from_osm.shp 各级道路OSM_河北省.shp \ -where highway IN (motorway,motorway_link) \ -nln 河北高速_from_osm \ -t_srs EPSG:3857 -lco ENCODINGUTF-8参数说明-where是 SQL 筛选条件注意字符串值用单引号整个条件用双引号包起来-nln指定输出图层名如果省略会默认用文件名-t_srs把坐标从原始 WGS84 转成 Web 墨卡托便于后续与在线底图叠加-lco ENCODINGUTF-8是创建选项让输出的 dbf 用 UTF-8 保存中文属性否则在 Linux 下很容易变成乱码。同样的命令可以把trunk、primary、secondary分别拆出或者反过来用-where highway NOT IN (...)排除掉人行道、自行车道、私家路。这里要留意 OSM 的motorway_link这类匝道如果不剔除后面统计里程时会和主线反复重叠导致长度虚高。3.3 用 Python 批量拆分并统计里程行政等级的五个文件本来就是按高速、国道、省道、县道、乡道拆开的你直接读取即可。但如果想统一汇总里程或者把某个混编层的不同级别批量输出成独立文件下面这段 geopandas 脚本更顺手import geopandas as gpd gdf gpd.read_file(各级道路OSM_河北省.shp, encodingUTF-8) if highway not in gdf.columns: raise KeyError(字段不存在请先确认属性字段名) # 先做必要清洗剔除缺失几何和空等级 gdf gdf[gdf.geometry.notnull()].copy() gdf gdf[gdf[highway].notna()] # 按等级批量写出 for level, subset in gdf.groupby(highway): if subset.empty: continue subset.to_file(fosm_{level}_河北省.shp, encodingUTF-8) print(f{level}: {len(subset)} 条, 总长 {subset.geometry.length.sum() / 1000:.1f} km)逻辑说明geopandas.read_file的encoding参数直接传给底层 fionagroupby(highway)会按等级分组每组写一个 Shp。geometry.length的单位取决于当前坐标系如果数据仍是 WGS84 经纬度这个值是“度”而不是米算出来的里程数字是没法直接用的。正确做法是先投影到适合河北的 CGCS2000 / UTM 50NEPSG:4547 或 32650再用投影坐标算长度。河北区域我常用 CGCS2000 3-degree Gauss-Kruger zone 39EPSG:4547精度好且是标准分带。gdf_utm gdf.to_crs(EPSG:4547) gdf_utm[length_km] gdf_utm.geometry.length / 1000注意to_crs之后长度就算正确了但这条结果只用于统计不要覆盖原始数据。保留一份 WGS84 经纬度版本方便后续和 OSM 原数据做空间连接。4. 公路等级与 OSM 分类映射空间对比和数据质量核查4.1 两套分级体系的关系行政等级道路高速、国道、省道、县道、乡道是交通行业按管理权属划分的OSM 的highway标签是按道路功能划分的比如motorway对应封闭快速路trunk对应主干快速路primary对应主干道。两者不矛盾但不能直接替换。项目里我见过有人拿 OSM 的motorway直接当高速用结果漏掉了一部分限速 100 的省级高速因为它们在 OSM 里被标成了trunk。反过来也有把trunk全当国道匹配导致国道里程虚高的情况。比较可靠的映射方式是把两套数据叠在一起做长度匹配率统计。下面是常见映射参考OSM highway 类别通常对应行政等级举例注意事项motorway / motorway_link高速主线及匝道有些高速路在 OSM 是 trunktrunk / trunk_link国道或城市快速路需要结合 name 和 ref 判断primary省道或城市主干道东部地区 primary 很多穿城secondary县道或城市次干道部分县道也会标 secondarytertiary乡道或城市支路tertiary 覆盖范围最杂residential / unclassified村道或内部道路和乡道有重叠但不等价4.2 用空间连接算重叠长度先来看一个具体疑问国道_河北省.shp和 OSM 里的trunk到底一对一吗我通常先算国道层与 OSM 各等级的叠加长度比。下面脚本把国道层按 OSM 等级打断后统计长度import geopandas as gpd national gpd.read_file(国道_河北省.shp, encodingUTF-8) osm gpd.read_file(各级道路OSM_河北省.shp, encodingUTF-8) # 统一坐标系避免几何类型歧义 national national.to_crs(EPSG:4547) osm osm.to_crs(EPSG:4547) # 空间连接重叠的线会被组合 joined gpd.overlay(national, osm[[highway, geometry]], howintersection) joined[len_km] joined.geometry.length / 1000 # 按 OSM 等级汇总被覆盖的长度 summary joined.groupby(highway)[len_km].sum().sort_values(ascendingFalse) print(summary)逻辑说明gpd.overlay是求两个面或线图层的几何交集对线图层来说它会把相交部分切分成新的短线。这里两次to_crs(EPSG:4547)的目的是让长度统计基于平面投影而不是经纬度。输出结果如果显示国道层和trunk重叠长度占据国道总长 75% 以上那说明该省国道和 OSMtrunk的对应关系较稳如果占比很低就要检查是不是坐标系偏了或者 OSM 在这片区域只画了primary。这类空间连接非常吃内存河北全省几万条线叠加后可能产生几十万个分段。数据量大的时候我建议先把两个输入spatial index打开geopandas 会自动用rtree不用手动干预。真正卡住时改用osmnx里的图论方法先把国道层转成网络图再按节点归属到 OSM 边复杂度从 O(N×M) 降到 O(NM)。4.3 检查拓扑错误悬挂点和自相交从网盘下载的 Shp 经常带着拓扑问题最常见的两种悬挂点一段乡道没有和任何其他道路连接。自相交同一条线自己穿过自己造成长度统计重复。检查悬挂点可以用 QGIS 的拓扑规则也可以在终端用 GDAL 的ogr2ogr配合 SQL 做一步傻瓜排查。我个人的偏好是先把数据导入 PostGIS然后跑一条 SQL 统计孤立线段SELECT a.gid, a.name, ST_IsSimple(a.geom) AS is_simple, ST_Length(a.geom) AS len FROM xiangdao a WHERE NOT ST_IsSimple(a.geom) LIMIT 20;说明ST_IsSimple返回 false 即表示存在自相交ST_Length返回投影系下的长度。如果是悬挂点则需要和邻近线段做ST_DWithin判断通常会设置 1 米容差SELECT count(*) FROM xiangdao a WHERE NOT EXISTS ( SELECT 1 FROM xiangdao b WHERE a.gid b.gid AND ST_DWithin(a.geom, b.geom, 1) );这里 1 的单位取决于当前坐标参考。如果表里存的是经纬度1 就是 1 度距离非常大必须改 0.00001 或先ST_Transform成投影坐标。和 3.3 节同样的原则任何拓扑检查前先统一坐标系。PostGIS 里我用 4547 或 3857前者按米算后者也是米但变形略大。4.4 OSM 与行政等级交叉修正方法发现两套数据存在差异后要以精度更高的为准。我的修正顺序通常是保留高速_河北省.shp里已经明确的类型不轻信 OSM 的motorway分类。对国道_河北省.shp与 OSMtrunk的差异段打开卫星影像目视确认。把确认无误的 OSMprimary中因新建而未收入行政等级文件的路段作为乡道以上补集纳入。这一步没有纯自动解法但可以用上述脚本先把差异段导成 kml叠加 Google Earth 看。输出用to_file写 kml 不如直接用ogr2ogr转ogr2ogr -f KML 差异路段.kml 差异路段.shp生成后放进 Google Earth 或 QGIS 的在线影像服务里核对比直接在代码里调影像 API 省事很多。5. 出图与分享前先把样式、编码和打包这三件事做对5.1 用 QGIS 样式文件统一展示规则一份路网数据如果直接叠加到地图上所有层都是同一种颜色根本没法看。我在 QGIS 里会给每类道路配一套 .qml 样式高速用红色实线国道用橙色省道用黄色县道用浅橙乡道用灰色细线城市一至四级道路用不同深浅的城镇色。具体操作上在图层右键进入Properties-Symbology选Categorized然后按字段type或highway分类。完成后Layer - Save Style To File保存为.qml。下次导入这批数据时直接加载 .qml就能省掉重复配色时间。QML 是一个 XML 文件内部不保存数据只存符号化规则。团队协作时把这个 .qml 和 Shp 一起放进同一个 rar 包拿到的人就不会再问“颜色为什么全是默认的”。所以要发原始数据给别人我建议目录结构这样安排河北道路数据/ ├── shp/ │ ├── 高速_河北省.shp │ ├── 国道_河北省.shp │ ├── ... (其余 shp 及配套文件) ├── style/ │ ├── 河北道路分级.qml │ └── 河北OSM路网.qml └── 数据说明.md5.2 重投影和编码统一避免交付后再修在 Windows 下用 ArcGIS 处理的 .shpdbf 编码常是 GBK。你如果直接拿去做 Web 服务会遇到中文乱码。解决办法是在导出时强制 UTF-8就像 2.2 节提到的那样。批量重投影并统一编码时我习惯写一个 Shell 循环for f in *_河北省.shp; do base${f%.shp} ogr2ogr -f ESRI Shapefile utf8_${base}.shp $f \ -t_srs EPSG:4547 \ -lco ENCODINGUTF-8 \ -lco PRECISIONNO done说明PRECISIONNO让 dbf 数字字段改用更紧凑的 NUMBER 类型避免小数点位数过多导致文件膨胀-t_srs EPSG:4547会把所有图层转到同一坐标系这样后续空间连接不用再统一。如果你要交付给 Web 端投影会改成 EPSG:3857如果要打印制图保留 4547 或 4326 更合适。5.3 边线简化做 Web 地图前的一个小技巧河北乡道 27000 多条全量塞给 Leaflet 会卡到崩溃。用mapshaper做简化是成本最低的路径。mapshaper是 Node 工具只需要一句安装命令npm i -g mapshaper然后用它简化并输出 GeoJSONmapshaper 乡道_河北省.shp \ -simplify dp 10% \ -o 乡道_简化.json formatgeojson-simplify dp指定使用 Douglas-Peucker 算法10%表示保留原图约 10% 的节点可视范围在乡镇级基本够用。如果你觉得 10% 还是太粗可以改成15%或20%但文件体积会线性增长。简化之前记得先投影地图服务器如果要求 WGS84最后再做一次mapshaper 乡道_河北省.shp \ -proj wgs84 \ -simplify dp 15% \ -o 乡道_简化.shp这样交付出去的数据在 QGIS、Web 端都能稳定显示。整套流程跑下来你会发现这套河北道路数据真正的价值在于它同时提供了行政等级和 OSM 两套参照系而能不能把它用好取决于你是否能对文件族、编码、投影和拓扑规则有足够的控制力。本文还有配套的精品资源点击获取
河北道路Shapefile数据处理：从文件家族到拆分出图完整指南
简介这份数据面向GIS分析、城市规划与交通研究从业者提供河北省最新分级道路矢量成果精度可达乡道级别。压缩包共97个文件以Shapefile格式为主包含.shp几何、.dbf属性、.prj投影、.shx索引等配套文件整包约73.95MB解压后可直接加载至ArcGIS或QGIS使用。资源涵盖城市一至四级道路、高速、国道、省道、县道、乡道以及OSM来源的铁路、地铁、轻轨、主干道、住宅街道等多种路网类型共16套道路数据不同来源数据相互补充既可进行专题制图也能用于缓冲区分析、路径规划等场景。目前已有312人学习下载适合需要按行政等级或功能等级分类调用道路数据的项目制作者参考。1. 一份 rar 解开后为什么全是同名文件拿到“河北省道路数据最新分级精确到乡道矢量数据.rar”的时候第一反应多半是里面怎么有几十个“同名”文件比如高速_河北省.shp、高速_河北省.dbf、高速_河北省.shx……这不是重复文件而是一整套 ESRI Shapefile。一个完整要素类本来就该由 .shp、.shx、.dbf、.prj 等一堆文件组成缺一个都可能打不开或者丢投影。这套数据把河北省道路按行政等级分成高速、国道、省道、县道、乡道又按城市内部等级分成一至四级道路还叠加了一份 OSM 来源的 roads 和 railways 线层总共 16 组线要素。对做路网分析、地图制图、城市规划和交通物流的人来说它解决的是“市级道路上哪找、县道乡道有没有、数据能不能直接进数据库”这一类问题。后文我会先从 Shapefile 文件家族说起再给出一套可以用命令行和 Python 直接跑的拆分、比对、出图流程并说明每个环节的参数怎么调、坑在哪。2. Shapefile 文件家族从 .shp 到 .sbn 一个个说清2.1 同一目录下每个扩展名都是干什么的Shapefile 不是单文件而是一个文件家族。这套河北道路数据里出现的扩展名基本把常见的都覆盖全了。先看一张我整理的文件与用途对照表扩展名作用缺失时的影响.shp要素几何本体记录每个线要素的形状文件打不开核心文件.shx形状索引用于快速定位 .shp 中的要素某些软件无法读取显示损坏.dbf属性表存道路名称、等级、类型等字段属性全丢只剩几何.prj投影与坐标系描述文件文本格式软件按 WGS84 猜测坐标可能错位.cpg属性表编码声明例如 UTF-8、GBK中文乱码的风险大幅上升.sbn / .sbx空间索引由 ArcGIS 生成可重建删除后下次打开会自动生成.shp.xml元数据XML 格式不影响使用只丢描述信息注意一个细节很多从网盘或公众号下回来的矢量数据.sbn和.sbx是成对出现的。它俩是 ArcGIS 在创建要素类时附带的空间索引你用 QGIS 或 GDAL 读取时并不会依赖它真正决定能否读取的是 .shp、.shx、.dbf 三个文件。所以如果解压后报缺失先把这三个文件确认齐了再谈其他。2.2 用 ogrinfo 快速验证数据完整性与坐标系我拿到数据的第一步从来不是直接拉进 QGIS而是在终端用ogrinfo看头部信息成本最低且能避开图形界面缓存问题。以其中一份乡道_河北省.shp为例ogrinfo -ro -so -al 乡道_河北省.shp参数说明-ro表示只读打开避免误写-so是 summary only只打印图层概要而不展开每条要素-al表示列出所有图层。输出里重点看这几行Layer name: 乡道_河北省 Geometry: Line String Feature Count: 27194 Extent: (113.608000, 36.051000) - (119.857000, 42.612000) Layer SRS: GEOGCS[WGS 84]如果Feature Count是 0多半是.shx和.dbf不配套如果Layer SRS显示未知说明缺失 .prj 或 .prj 为空。此时我用另一条命令重新指定坐标系ogrinfo -ro -al -so -oo ENCODINGUTF-8 乡道_河北省.shp-oo ENCODING是 open option告诉 GDAL 按 UTF-8 解码 dbf能解决一部分中文乱码。实际项目中我遇到过 .cpg 写的是 UTF-8但内部实际是 GBK 的混乱情况单看扩展名没用要靠输出里属性字段的值是否可读来判断。2.3 用 Python 列出字段和统计要素数ogrinfo 适合快速看但要做批量检查我更喜欢用 Python 的fiona或pyogrio一次遍历整个目录。下面脚本会扫描所有*_河北省.shp打印每个图层的要素数量、几何类型以及 dbf 里有哪些字段import glob import fiona for shp in glob.glob(*_河北省.shp): with fiona.open(shp, encodingUTF-8) as src: schema src.schema props list(schema[properties].keys()) print(f{shp}: {len(src)} 个要素, 类型{schema[geometry]}) print(f 字段: {props[:8]}{... if len(props) 8 else })逻辑说明fiona.open的encoding参数对应 dbf 的读入编码和 GDAL 的ENCODING等价schema[properties]是一个有序字典字段顺序与 dbf 物理顺序一致。这套数据里的“各级道路OSM_河北省.shp”字段较多通常包含highway、name、ref、oneway等 OSM 标准字段而行政等级道路层字段相对精简常见的是name、type、length。先看清字段名后面做拆分筛选才有依据。3. 按行政等级拆路网高速、国道、省道、县道、乡道一次拆干净3.1 确认分级字段的值域忻州、保定这种多等级混合的区域道路层如果已经按文件拆分那还相对省事。但实际解压后你会发现既有高速_河北省.shp又有各级道路OSM_河北省.shp后者是把所有等级混在一个图层里的。所以第一步永远是确认分层文件里的等级字段到底是什么值。用 2.3 节的脚本先看一遍然后对某个混编文件执行ogrinfo -ro -al -so -sql SELECT DISTINCT highway FROM 各级道路OSM_河北省 各级道路OSM_河北省.shp这里如果highway字段不存在会直接报no such column。报错时先跑ogrinfo -ro -al -so 各级道路OSM_河北省.shp把真实字段名打出来常见的是type、fclass、kind。OSM 来源的数据里highway的典型值包括motorway、trunk、primary、secondary、tertiary、residential、track等。行政等级和 OSM 等级是两个体系不能简单画等号。3.2 用 ogr2ogr 按属性拆出单独图层当你确认了分级字段之后拆分就很简单。下面命令把 OSM 全量路网按highwaymotorway拆出来另存为独立 Shpogr2ogr -f ESRI Shapefile 河北高速_from_osm.shp 各级道路OSM_河北省.shp \ -where highway IN (motorway,motorway_link) \ -nln 河北高速_from_osm \ -t_srs EPSG:3857 -lco ENCODINGUTF-8参数说明-where是 SQL 筛选条件注意字符串值用单引号整个条件用双引号包起来-nln指定输出图层名如果省略会默认用文件名-t_srs把坐标从原始 WGS84 转成 Web 墨卡托便于后续与在线底图叠加-lco ENCODINGUTF-8是创建选项让输出的 dbf 用 UTF-8 保存中文属性否则在 Linux 下很容易变成乱码。同样的命令可以把trunk、primary、secondary分别拆出或者反过来用-where highway NOT IN (...)排除掉人行道、自行车道、私家路。这里要留意 OSM 的motorway_link这类匝道如果不剔除后面统计里程时会和主线反复重叠导致长度虚高。3.3 用 Python 批量拆分并统计里程行政等级的五个文件本来就是按高速、国道、省道、县道、乡道拆开的你直接读取即可。但如果想统一汇总里程或者把某个混编层的不同级别批量输出成独立文件下面这段 geopandas 脚本更顺手import geopandas as gpd gdf gpd.read_file(各级道路OSM_河北省.shp, encodingUTF-8) if highway not in gdf.columns: raise KeyError(字段不存在请先确认属性字段名) # 先做必要清洗剔除缺失几何和空等级 gdf gdf[gdf.geometry.notnull()].copy() gdf gdf[gdf[highway].notna()] # 按等级批量写出 for level, subset in gdf.groupby(highway): if subset.empty: continue subset.to_file(fosm_{level}_河北省.shp, encodingUTF-8) print(f{level}: {len(subset)} 条, 总长 {subset.geometry.length.sum() / 1000:.1f} km)逻辑说明geopandas.read_file的encoding参数直接传给底层 fionagroupby(highway)会按等级分组每组写一个 Shp。geometry.length的单位取决于当前坐标系如果数据仍是 WGS84 经纬度这个值是“度”而不是米算出来的里程数字是没法直接用的。正确做法是先投影到适合河北的 CGCS2000 / UTM 50NEPSG:4547 或 32650再用投影坐标算长度。河北区域我常用 CGCS2000 3-degree Gauss-Kruger zone 39EPSG:4547精度好且是标准分带。gdf_utm gdf.to_crs(EPSG:4547) gdf_utm[length_km] gdf_utm.geometry.length / 1000注意to_crs之后长度就算正确了但这条结果只用于统计不要覆盖原始数据。保留一份 WGS84 经纬度版本方便后续和 OSM 原数据做空间连接。4. 公路等级与 OSM 分类映射空间对比和数据质量核查4.1 两套分级体系的关系行政等级道路高速、国道、省道、县道、乡道是交通行业按管理权属划分的OSM 的highway标签是按道路功能划分的比如motorway对应封闭快速路trunk对应主干快速路primary对应主干道。两者不矛盾但不能直接替换。项目里我见过有人拿 OSM 的motorway直接当高速用结果漏掉了一部分限速 100 的省级高速因为它们在 OSM 里被标成了trunk。反过来也有把trunk全当国道匹配导致国道里程虚高的情况。比较可靠的映射方式是把两套数据叠在一起做长度匹配率统计。下面是常见映射参考OSM highway 类别通常对应行政等级举例注意事项motorway / motorway_link高速主线及匝道有些高速路在 OSM 是 trunktrunk / trunk_link国道或城市快速路需要结合 name 和 ref 判断primary省道或城市主干道东部地区 primary 很多穿城secondary县道或城市次干道部分县道也会标 secondarytertiary乡道或城市支路tertiary 覆盖范围最杂residential / unclassified村道或内部道路和乡道有重叠但不等价4.2 用空间连接算重叠长度先来看一个具体疑问国道_河北省.shp和 OSM 里的trunk到底一对一吗我通常先算国道层与 OSM 各等级的叠加长度比。下面脚本把国道层按 OSM 等级打断后统计长度import geopandas as gpd national gpd.read_file(国道_河北省.shp, encodingUTF-8) osm gpd.read_file(各级道路OSM_河北省.shp, encodingUTF-8) # 统一坐标系避免几何类型歧义 national national.to_crs(EPSG:4547) osm osm.to_crs(EPSG:4547) # 空间连接重叠的线会被组合 joined gpd.overlay(national, osm[[highway, geometry]], howintersection) joined[len_km] joined.geometry.length / 1000 # 按 OSM 等级汇总被覆盖的长度 summary joined.groupby(highway)[len_km].sum().sort_values(ascendingFalse) print(summary)逻辑说明gpd.overlay是求两个面或线图层的几何交集对线图层来说它会把相交部分切分成新的短线。这里两次to_crs(EPSG:4547)的目的是让长度统计基于平面投影而不是经纬度。输出结果如果显示国道层和trunk重叠长度占据国道总长 75% 以上那说明该省国道和 OSMtrunk的对应关系较稳如果占比很低就要检查是不是坐标系偏了或者 OSM 在这片区域只画了primary。这类空间连接非常吃内存河北全省几万条线叠加后可能产生几十万个分段。数据量大的时候我建议先把两个输入spatial index打开geopandas 会自动用rtree不用手动干预。真正卡住时改用osmnx里的图论方法先把国道层转成网络图再按节点归属到 OSM 边复杂度从 O(N×M) 降到 O(NM)。4.3 检查拓扑错误悬挂点和自相交从网盘下载的 Shp 经常带着拓扑问题最常见的两种悬挂点一段乡道没有和任何其他道路连接。自相交同一条线自己穿过自己造成长度统计重复。检查悬挂点可以用 QGIS 的拓扑规则也可以在终端用 GDAL 的ogr2ogr配合 SQL 做一步傻瓜排查。我个人的偏好是先把数据导入 PostGIS然后跑一条 SQL 统计孤立线段SELECT a.gid, a.name, ST_IsSimple(a.geom) AS is_simple, ST_Length(a.geom) AS len FROM xiangdao a WHERE NOT ST_IsSimple(a.geom) LIMIT 20;说明ST_IsSimple返回 false 即表示存在自相交ST_Length返回投影系下的长度。如果是悬挂点则需要和邻近线段做ST_DWithin判断通常会设置 1 米容差SELECT count(*) FROM xiangdao a WHERE NOT EXISTS ( SELECT 1 FROM xiangdao b WHERE a.gid b.gid AND ST_DWithin(a.geom, b.geom, 1) );这里 1 的单位取决于当前坐标参考。如果表里存的是经纬度1 就是 1 度距离非常大必须改 0.00001 或先ST_Transform成投影坐标。和 3.3 节同样的原则任何拓扑检查前先统一坐标系。PostGIS 里我用 4547 或 3857前者按米算后者也是米但变形略大。4.4 OSM 与行政等级交叉修正方法发现两套数据存在差异后要以精度更高的为准。我的修正顺序通常是保留高速_河北省.shp里已经明确的类型不轻信 OSM 的motorway分类。对国道_河北省.shp与 OSMtrunk的差异段打开卫星影像目视确认。把确认无误的 OSMprimary中因新建而未收入行政等级文件的路段作为乡道以上补集纳入。这一步没有纯自动解法但可以用上述脚本先把差异段导成 kml叠加 Google Earth 看。输出用to_file写 kml 不如直接用ogr2ogr转ogr2ogr -f KML 差异路段.kml 差异路段.shp生成后放进 Google Earth 或 QGIS 的在线影像服务里核对比直接在代码里调影像 API 省事很多。5. 出图与分享前先把样式、编码和打包这三件事做对5.1 用 QGIS 样式文件统一展示规则一份路网数据如果直接叠加到地图上所有层都是同一种颜色根本没法看。我在 QGIS 里会给每类道路配一套 .qml 样式高速用红色实线国道用橙色省道用黄色县道用浅橙乡道用灰色细线城市一至四级道路用不同深浅的城镇色。具体操作上在图层右键进入Properties-Symbology选Categorized然后按字段type或highway分类。完成后Layer - Save Style To File保存为.qml。下次导入这批数据时直接加载 .qml就能省掉重复配色时间。QML 是一个 XML 文件内部不保存数据只存符号化规则。团队协作时把这个 .qml 和 Shp 一起放进同一个 rar 包拿到的人就不会再问“颜色为什么全是默认的”。所以要发原始数据给别人我建议目录结构这样安排河北道路数据/ ├── shp/ │ ├── 高速_河北省.shp │ ├── 国道_河北省.shp │ ├── ... (其余 shp 及配套文件) ├── style/ │ ├── 河北道路分级.qml │ └── 河北OSM路网.qml └── 数据说明.md5.2 重投影和编码统一避免交付后再修在 Windows 下用 ArcGIS 处理的 .shpdbf 编码常是 GBK。你如果直接拿去做 Web 服务会遇到中文乱码。解决办法是在导出时强制 UTF-8就像 2.2 节提到的那样。批量重投影并统一编码时我习惯写一个 Shell 循环for f in *_河北省.shp; do base${f%.shp} ogr2ogr -f ESRI Shapefile utf8_${base}.shp $f \ -t_srs EPSG:4547 \ -lco ENCODINGUTF-8 \ -lco PRECISIONNO done说明PRECISIONNO让 dbf 数字字段改用更紧凑的 NUMBER 类型避免小数点位数过多导致文件膨胀-t_srs EPSG:4547会把所有图层转到同一坐标系这样后续空间连接不用再统一。如果你要交付给 Web 端投影会改成 EPSG:3857如果要打印制图保留 4547 或 4326 更合适。5.3 边线简化做 Web 地图前的一个小技巧河北乡道 27000 多条全量塞给 Leaflet 会卡到崩溃。用mapshaper做简化是成本最低的路径。mapshaper是 Node 工具只需要一句安装命令npm i -g mapshaper然后用它简化并输出 GeoJSONmapshaper 乡道_河北省.shp \ -simplify dp 10% \ -o 乡道_简化.json formatgeojson-simplify dp指定使用 Douglas-Peucker 算法10%表示保留原图约 10% 的节点可视范围在乡镇级基本够用。如果你觉得 10% 还是太粗可以改成15%或20%但文件体积会线性增长。简化之前记得先投影地图服务器如果要求 WGS84最后再做一次mapshaper 乡道_河北省.shp \ -proj wgs84 \ -simplify dp 15% \ -o 乡道_简化.shp这样交付出去的数据在 QGIS、Web 端都能稳定显示。整套流程跑下来你会发现这套河北道路数据真正的价值在于它同时提供了行政等级和 OSM 两套参照系而能不能把它用好取决于你是否能对文件族、编码、投影和拓扑规则有足够的控制力。本文还有配套的精品资源点击获取
