ARTICLE DETAIL

资讯详情

深耕郑州网站建设与运营推广的一线实战洞察。

内蒙古乡道路网矢量数据提取与ArcSWAT分析实战

内蒙古乡道路网矢量数据提取与ArcSWAT分析实战 简介内蒙古自治区道路矢量数据包是以最新道路分级成果为核心的RAR压缩包数据精确到乡道面向GIS分析师、交通规划工程师以及地图制图与空间数据研究者。压缩包内收录16种道路矢量数据既有按城市等级划分的一级至四级道路线也有按行政等级划分的高速、国道、省道、县道、乡道线还包含OpenStreetMap来源的铁路、轻轨、窄轨、地铁、有轨电车及主干道、次干道、支道、人行道、住宅街道、自行车道等图层覆盖面广部分数据相互重叠又能交叉验证方便考证。RAR文件整体大小39.26MB便携易用。目前已有216人浏览学习该数据集。使用这份路网数据可进行区域路网结构分析、交通可达性评估、专题制图与空间可视化分类细致、字段完整适合科研、规划、教学及GIS二次开发时作为基础底图。1. 内蒙古自治区道路数据拿到手这包矢量数据到底能做什么做内蒙古这边的路网分析、ArcSWAT 小流域分析或者交通可达性的时候最怕的不是没有数据而是拿到一份“似乎能用”的矢量数据解压后根本不知道里面是什么坐标系、有没有属性表、乡道字段到底叫什么。这份内蒙古自治区道路数据最新分级精确到乡道矢量数据.rar本质上是一份全自治区范围的线状路网矢量数据覆盖高速、国道、省道、县道和乡道精度到乡镇一级。对规划、水利、环保、农牧业方向的从业者来说它的核心价值在于不用再拿 OSM 那种噪点很多的道路网凑合跑水文模型也不用去不同旗县拼凑路网边界。常见做法是直接把它当成 ArcGIS / QGIS 的基础底图做“道路密度 vs 小流域产流”“路网通达性 vs 草地退化”这类空间统计或者把它转成栅格成本数据喂给 ArcSWAT 做汇流路径分析。适合谁适合手上只有 DEM、土地利用、气象站数据但缺一套统一标准路网的人。这篇笔记先讲怎么把 .rar 里的矢量数据顺利解压、识别、分层再讲字段和坐标系的坑最后给出 ArcSWAT 场景下的落地用法。2. 从 .rar 到可用图层解压、目录识别与数据体检2.1 解压前先做三件事看文件列表、看压缩注释、核对压缩格式拿到这份 .rar不要急着双击解压到桌面。城市规划和 GIS 数据工程师的习惯是先把压缩包当做一个“黑匣子”做一次体检。RAR 压缩包里有可能是单层文件夹也有可能是嵌套四五层的目录树甚至压缩包本身是分卷包的一部分直接解压会得到一堆半截文件。我的建议是分三步走。第一步用 WinRAR 或 7-Zip 打开压缩包先看文件列表确认里面是 .shp、.dbf、.shx、.prj还是一个 File Geodatabase 文件夹或者干脆是 .gdb 的压缩形式。第二步看压缩包的注释信息和修改时间这能帮你判断数据是哪个版本、谁发布的。内蒙古旗县一级的路网数据更新频繁名称里带“最新”不代表时间戳真的靠前解压后确认内部文件的修改时间是更可靠的依据。第三步确认压缩格式是不是标准 RAR有些下载站点会用自解压格式 .exe 冒充 .rar或者压缩时用了高版本 RAR5 算法老版本 WinRAR 解出来会提示“未知方式”这时候需要换用 7-Zip 21.0 以上版本。# 在 Windows 命令行下用 7-Zip 列出压缩包内容不实际解压 7z l 内蒙古自治区道路数据最新分级精确到乡道矢量数据.rar这条命令的关键是l参数它只列出清单不解压能看到包内文件的完整路径和大小。如果输出直接是一堆 .shp/.dbf/.shx说明是传统 Shapefile如果列出的是一个 .gdb 文件夹则说明是全要素地理数据库。观察路径长度也很重要Shapefile 对路径敏感超过 256 字符在部分工具里会读取失败。解压时我一般会在 D 盘建一个短路径目录例如D:\road_data\避免“解压后发现 ArcGIS 无法添加图层”这种最常见翻车。2.2 打开后怎么判断这套数据的“底细”Shp、GDB、要素类与属性表解压完成后先别急着拖进 ArcMap。先看根目录下都有什么。如果是 Shapefile 格式至少应该有 .shp几何、.shx索引、.dbf属性、.prj坐标系四个文件。缺了 .prj 的矢量数据基本等于坐标参考黑匣子必须马上在 ArcToolbox 里用 Define Projection 去查证坐标系。如果看到的是 .gdb 文件夹那就简单很多gdb 是一个完整的数据库容器用 ArcGIS Pro 目录面板直接拖进地图即可不需要逐文件操作。需要留意的还有数据是线要素还是面要素。标题写的是“道路数据”道路通常是以线要素表达的但有的数据源会把乡道、村道做成面状路网尤其在农村公路数据库中比较常见。如果拖进地图后看到的是细长多边形别吃惊先检查要素类类型再做转换。用 Python 写一段体检脚本是最快的# 体检脚本读取shp或gdb要素类的基本信息 import geopandas as gpd # 如果是shapefile直接读文件路径 # 如果是gdb路径格式示例D:/road_data/data.gdb/road_2024 gdf gpd.read_file(rD:/road_data/road.shp, encodingutf-8) print(要素数量, len(gdf)) print(几何类型, gdf.geom_type.unique()) print(坐标范围, gdf.total_bounds) print(字段列表, list(gdf.columns)) # 检查坐标系 print(坐标系, gdf.crs)这段脚本能一次告诉你五件事要素数量是不是准确反映了全区路网、几何类型是线还是面、四至范围是不是覆盖内蒙古全境、属性字段是不是符合常识、坐标系是不是 CGCS2000。注意encoding参数内蒙古的 shp 属性表很多用的是 GBK 编码如果读出来中文字段名乱码大概率要把utf-8改成gbk再读。这一步不是选择题而是必考题因为后续筛选乡道全靠字段名和字段值。2.3 用统计命令确认分级字段完整度不打开 GIS 也能看属性很多人解压后就打开 ArcMap 看图层样式这太慢了。更快的方式是直接用 Python 对属性表做分组统计看看“等级”字段里到底有哪几类道路数量各是多少。这样不用等地图渲染完成就能马上判断这份数据是否真的“精确到乡道”。继续用 geopandas# 统计道路等级字段的分布判断数据是否覆盖乡道 import geopandas as gpd gdf gpd.read_file(rD:/road_data/road.shp, encodinggbk) # 常见字段名可能是 fld_layer / type / 道路等级不同数据源命名不同 # 先打印所有字段再挑一个可能是分级字段的列 print(gdf.columns.tolist()) level_col type # 按实际字段名修改这里以 type 为例 level_counts gdf[level_col].value_counts() print(level_counts) # 输出每个等级的长度单位是米方便核对数据量 gdf[length_m] gdf.geometry.length stats gdf.groupby(level_col)[length_m].agg([count, sum]) stats[sum_km] stats[sum] / 1000 print(stats)这段代码背后有三个逻辑。第一value_counts只看类别数量能快速暴露“只有高速和国道没有乡道”这种数据不完整的情况。第二geometry.length得到的是投影坐标系下的长度单位是米如果坐标系是地理坐标经纬度这个长度单位是度数值无意义所以在读统计前先确认坐标系。第三按等级汇总里程可以用来和官网公开的全区公路里程做交叉验证。内蒙古的乡道和村道总里程在公开资料里能查到大致的量级如果你算出来乡道总里程只有几百公里那这份数据大概率只抽取了部分旗县并不是全区覆盖。数据体检做到这一步花不到两分钟省下的是后面浪费时间跑模型却发现结果不可信的成本。3. 字段与坐标系把乡道从属性表里正确捞出来3.1 看懂分级字段高速、国道、省道、县道、乡道的典型编码内蒙古道路数据的分级字段不同批次数据命名差异很大。常见的有fld_layer、type、kind、道路等级、class这几种。字段值也不统一有的用中文“乡道”有的用拼音缩写“XD”有的用数字编码“5”。很多新手会直接用“乡道”字符串去筛选结果一条都选不出来然后怀疑数据有问题。其实大部分情况下数据没毛病是字段值用了拼音或编码。我的建议是打开属性表后先对分级字段做一次唯一值去重把每个类别的样本看一遍。不要猜。最快的方法是直接用 pandas 的unique()或者用 QGIS 的字段统计。看到唯一值之后再建立一套等级映射规则。import pandas as pd import geopandas as gpd gdf gpd.read_file(rD:/road_data/road.shp, encodinggbk) # 假设字段名是 type先打印唯一值 print(gdf[type].unique()) # 建立映射字典把各种写法归一化成统一等级 mapping { G: 国道, 国道: 国道, national: 国道, S: 省道, 省道: 省道, X: 县道, 县道: 县道, 乡道: 乡道, XD: 乡道, C: 村道, 村道: 村道, } # 新增一列 level gdf[level] gdf[type].map(mapping) # 检查映射后还有多少空值空值说明有未覆盖到的情况 print(gdf[level].isna().sum()) print(gdf[gdf[level].isna()][type].unique())这段代码的精髓在打印空值这一步。映射之后出现空值说明字段里存在你没预料到的写法必须回头再看一眼原始值。常见情况是高速被写成“GX”或者“G45”这种具体编号国道路线编号绝对不能和“国道”这个等级混在一列里路线编号是属性等级是类别。建议把“等级”和“路线编号”分两个字段处理因为 ArcSWAT 环境里做成本栅格时只关心道路等级对应的道路宽度和行车速度不需要路线编号参与计算。3.2 坐标系的坑CGCS2000、西安1980、投影坐标系与经纬度内蒙古的道路矢量数据常见有三种坐标系。第一种是 CGCS2000 地理坐标系单位是度代号 4490第二种是 CGCS2000 / 3-degree Gauss-Kruger按带号区分例如 CM 117E代号带 ESRI 的字符串第三种是最头痛的——西安 1980 坐标系。如果 .prj 文件丢失数据在 ArcGIS 里会显示为 Unknown Units这时候叠加行政区划边界会发现错位几十公里到几百米不等。拿到数据的第一步就是看crs。如果是 4490直接to_crs转投影。如果是西安 1980内蒙古范围跨度大需要先确定它用的是 6 度分带还是 3 度分带。内蒙东西经度跨度大不能简单用同一个中央经线处理建议用gdf.total_bounds计算中心经度再选择合适的投影带而不是盲目选择“全国统一”的投影。下面这段代码展示了如何做坐标系转换和验证import geopandas as gpd gdf gpd.read_file(rD:/road_data/road.shp, encodinggbk) # 查看原始坐标系如果返回 None 或空说明 .prj 缺失 print(原始坐标系, gdf.crs) # 如果缺失坐标系先用边界四至判断大概是哪种坐标系 # 内蒙古东西向跨度大min_x 接近 97max_x 接近 126是经纬度范畴 # 如果 min_x 是几百到几千数量级那已经是投影坐标系 bounds gdf.total_bounds print(四至minx,miny,maxx,maxy, bounds) # 经纬度坐标下统一转成 CGCS2000 3度分带投影中央经线取120度 # EPSG:4547 是 CGCS2000 / 3-degree Gauss-Kruger CM 120E if bounds[0] 90 and bounds[2] 140: gdf gdf.to_crs(epsg4547) else: # 如果已经是投影坐标先确保和理解坐标系一致不要反复转换 print(已经是投影坐标系, gdf.crs) # 重新计算长度这时候单位是米 gdf[length_m] gdf.geometry.length print(gdf[[type, length_m]].groupby(type).sum() / 1000)参数说明EPSG 4547 是 CGCS2000 三度带 120° 中央经线的投影坐标适合内蒙古中东部如果数据主要集中在阿拉善盟、乌海一带中央经线 105° 的投影带更合适对应 EPSG 4545。别嫌麻烦投影选错后面所有面积和长度统计全部失真。ArcSWAT 小流域分析中道路密度计算如果用了不同投影的数据结果偏差可能超过 20%。3.3 按等级提取乡道的两种做法QGIS 表达式与 Python 筛选体质筛查完成后开始提取乡道。两种主流工具二选一如果你已经打开 QGIS用表达式更直接如果你还在 Python 环境里用属性筛选再导出。先讲 QGIS 的做法。在图层上右键打开属性表点击“按表达式选择要素”输入表达式type IN (乡道, XD)这里有一个坑属性表字段名如果是中文字段在 QGIS 表达式里必须用双引号括起来字符串值用单引号。如果你在 QGIS 里写表达式时中文字段输入不了检查图层编码设置矢量图层的编码在“图层属性 - 源 - 编码”里改成 GBK否则中文值永远匹配不上。筛选结果出来后右键图层 - 导出 - 保存所选要素为选择 GeoPackage 格式这样避免 Shapefile 对字段名长度和编码的历史包袱。Python 的做法则适合批处理场景import geopandas as gpd gdf gpd.read_file(rD:/road_data/road.shp, encodinggbk) # 归一化之后直接筛选乡道 gdf[level] gdf[type].map(mapping) township_roads gdf[gdf[level] 乡道].copy() print(f乡道要素数量{len(township_roads)}) # 导出为GeoPackage保留坐标系和属性字段 township_roads.to_file(rD:/road_data/township_roads.gpkg, driverGPKG, encodingutf-8)这里解释一下为什么导出用 GeoPackage 而不是 Shapefile。GeoPackage 是 SQLite 容器字段名长度限制宽松不会出现中文字段变成fld_1、fld_2的问题。另外导出后随手再读一次check gpd.read_file(rD:/road_data/township_roads.gpkg) print(check.crs) print(check[level].value_counts())这一步是验证闭环因为你不能保证刚才的筛选表达式一定正确重新读回来检查属性分布才是最稳的收尾。4. 避坑内蒙古道路矢量数据的五个经典翻车现场4.1 压缩包内嵌套多层目录导致路径太深无法加载现象解压后把 .shp 拖进 ArcGIS提示“无法添加数据”或者图层能添加但符号系统全是灰色。原因压缩包内目录嵌套过深shapefile 全路径加文件名超过 Windows 最长路径限制更常见的是 .shp、.dbf、.shx 分散在多层目录里而用户只拷贝了其中一个文件到新目录。解决先解压到短路径如D:\road_data\如果解压后是三层以上目录直接把最底层的 shp 文件移动到根目录层级并确认 .shp、.dbf、.shx、.prj 四个兄弟文件在同一级目录且文件名完全一致。4.2 数据叠加后错位几十公里坐标系不明是黑匣子现象道路图层和旗县行政区划边界叠加后整体偏移有的路段跑到山脊上。原因道路数据是西安 1980 坐标系或北京 1954而底图是 CGCS2000且原始 .prj 缺失软件默认按无投影显示。解决先用total_bounds判断是经纬度还是投影再通过控制点方式做坐标纠正。内蒙古常用的纠偏手段是找几个已知的旗县行政中心点做位移计算计算横向和纵向偏移量常量然后用translate平移。注意这种平移只能解决同椭球下的投影带偏移不能解决不同椭球系统的高精度对齐真正要做数据转换时优先用 ArcGIS 的投影工具而不是手动平移。4.3 乡道字段值不统一筛选时漏掉大半现象用“乡道”筛选结果只筛出几百条但属性表里明明有很多乡级道路。原因字段值存在混写同一条路有的记录是“乡道”有的记录是“XD”有的是“5”。解决不要直接筛先unique()全量枚举值再建映射字典。一个字都不能省。“县道”可能写成“X”、“县道”、“XianDao”“乡道”可能写成“C”或“VillageRoad”。这个坑在新手区出现频率极高因为很多人直接拿“乡道”两个字去选结果做出错误结论还以为数据不行。4.4 属性表中文乱码字段名变成未知字符现象用 QGIS 打开后属性表里的中文显示成乱码或者字段名变成???。原因Shapefile 的 .dbf 编码是 GBK用 UTF-8 打开当然乱。ArcGIS 环境有时候能自动识别QGIS 则经常猜错。解决在 QGIS 图层属性里手动指定编码为 GBKPython 读取时设置encodinggbk重新导出后统一转成 UTF-8。这里还容易栽一个跟头加了encoding之后to_file导出时还需要再设置一次编码两个环节都设置对才能保证中文入库后不乱。4.5 道路几何存在断裂做网络分析前必须拓扑检查现象乡道提取完进行网络分析时提示“没有连接到任何源”或者路径规划完全找不到路。原因线要素之间存在微小空隙或悬挂点放大后才会发现道路网并不是一个连通图。解决使用 ArcGIS 的拓扑修复或者 QGIS 的v.clean工具用 GRASS 的break和snap参数处理。最常见参数组合是snap0.0001和threshold0.0001这个阈值在投影坐标系下是 0.1 毫米足够修复航摄误差造成的微小断裂又不至于把平行双车道粘连成一条线。处理后再用网络分析工具做一次连通性测试而不是用眼睛目测。5. 把这份数据用出价值从乡道提取到 ArcSWAT 流域分析的落地技巧5.1 做 ArcSWAT 小流域分析时乡道矢量数据到底以什么方式进入模型ArcSWAT 做小流域分析需要输入 DEM、土壤、土地利用和气象数据道路数据本身通常不直接作为 SWAT 的输入参数。但它是制作“汇流路径”和“道路密度”的关键辅助数据。常见做法是把道路按等级转换成栅格成本面用来模拟地表径流被道路截断、加速汇流的过程。乡道和村道在内蒙古草原地区往往是地表径流的“人工排水通道”路基两侧的边坡和边沟会改变汇流方向。我一般先把乡道提取出来按等级对道路宽度和曼宁系数赋值然后转栅格作为 HRU 生成时的土地利用附属因子参与叠加。操作路径是ArcToolbox - Conversion Tools - To Raster - Polyline to Raster字段选“等级”像元大小设为 30 米和 DEM 重采样分辨率一致。然后重分类成成本栅格高速成本低、乡道成本高再用成本路径工具修正 SWAT 默认的河道输出位置。注意 ArcSWAT 里的“流域”是基于 DEM 水文分析生成的道路数据只起辅助修正作用不能替代 DEM 的流向计算。5.2 把乡道矢量数据转成栅格成本面的参数设置与验证结果用Polyline to Raster时有三个必调参数。第一个是像元大小不能比 DEM 小太多否则计算量爆炸30 米是一个平衡折中。第二个是字段必须用分级字段不能把路线编号用作栅格值否则每个路号都成独立类别。第三个是像元类型选 8 bit unsigned保证重分类时类别数量够用且文件体积小。转完栅格后做一次成本面合理性验证。打开栅格属性表统计每一类像元数量再和矢量属性表里各等级道路的里程做交叉对比。用 Python 打印一段验证结果import rasterio import numpy as np with rasterio.open(rD:/road_data/road_cost.tif) as src: data src.read(1) # 统计不同等级像元数 unique, counts np.unique(data, return_countsTrue) res dict(zip(unique, counts)) print(不同等级栅格值统计, res) print(有效像元数, np.count_nonzero(data)) print(栅格坐标系, src.crs) print(像元大小, src.transform[0])这个输出能帮你判断栅格是否和矢量范围契合。如果有效像元数太少大概率是几何转换时范围没对齐需要用 DEM 的栅格范围去裁剪道路成本面再参与 ArcSWAT 叠加。还有一种坑投影坐标系不一致导致栅格转出来是倾斜的打印src.crs就能一眼看出问题。5.3 用里程统计反向验证数据可信度三个数字对上了才算数数据可信度验证不能只看“有没有数据”要看数值量级是否合理。内蒙古自治区总面积约 118 万平方公里乡道和村道的总里程在公开统计中体量很大如果你的乡道提取结果只有两三百公里数据覆盖率一定有问题。我常用三个指标互相验证第一全区道路总里程和你手上分层结果的加和第二乡道里程占县道以下道路总里程的比例第三单位面积乡道里程密度是否和当地旗县的公路网密度一致。写一段短脚本算这三个指标import geopandas as gpd gdf gpd.read_file(rD:/road_data/road.shp, encodinggbk) gdf[length_km] gdf.geometry.length / 1000 # 指标1各省道、县道、乡道总里程 print(gdf.groupby(level)[length_km].sum()) # 指标2乡道在县道以下路网中的比例 total_lower gdf.loc[gdf[level].isin([县道, 乡道]), length_km].sum() township_km gdf.loc[gdf[level] 乡道, length_km].sum() print(乡镇道路占比, township_km / total_lower) # 指标3按旗县代码统计乡道密度观察各地差异 area_km2 1180000 # 内蒙古总面积任意区域细化分析时改用具体旗县面积 total_km gdf[length_km].sum() print(全区路网总密度, total_km / area_km2)这里的目的是为 ArcSWAT 提供一张可解释的道路成本面而不是为了追求报告里一个“精确”的数字。如果你发现乡道占县道以下路网比例超过 60%反而要警惕可能县道被错误归入乡道类别。别把验证当作走过场ArcSWAT 的小流域分析结果对道路密度非常敏感成本面错了流域划分可能整体偏掉。5.4 进阶用法把乡道网络转成道路密度面参与 HRU 的划分决策完成成本面之后还有一个进阶操作用核密度工具生成道路密度栅格叠加到 ArcSWAT 的 HRU 分类中。具体做法是把乡道转成点要素每 100 米生成一个点再用 ArcGIS 的 Kernel Density 工具搜索半径设为 3000 米输出栅格像元大小仍然 30 米。这一步解决的实际问题是SWAT 的 HRU 生成过程默认不看重道路因素只依据土地利用、土壤和坡度三重叠加。你把道路密度作为一个辅助分类条件就能在草原牧区把邻近道路的高强度放牧区和远离道路的天然草地区区分开。这在内蒙古很多旗县的水文模拟中非常实用因为道路周边往往是牲畜集中活动区域土壤压实程度和地表糙率差异极大。道路密度面叠加进 HRU 分类的具体操作是在 ArcSWAT 的 Land Use / Soil / Slope 定义界面把重分类后的道路密度栅格加入“HRU 附加分类层”并赋一个阈值例如密度大于 0.001 的像元划入“牧道影响区”。阈值可以按你的研究尺度调节草原流域建议先试 2000 米搜索半径因为乡道间距大半径太小密度面会碎成斑点。这一步做完SWAT 的模拟结果会增加一个解释维度也能让你的分析在答辩和评审时更有说服力——因为你不再只是拿一套通用路网凑数而是真正把乡道的空间分布写进了流域划分逻辑里。最后想提醒你一句道路矢量数据不管多新都不能直接当精度验证标准拿它做分析之前先跑到野外抽几条乡道拍拍照、对一下 GPS 轨迹这份数据在你手里才算真正活起来。希望帮到你。本文还有配套的精品资源点击获取
返回列表