★ FEATURED ARTICLE 简介北京城区道路矢量数据包提供主城区范围内的精细化路网覆盖主干道、次干道、支路及部分街巷适合地理信息、城市规划与地图制图人员直接使用。压缩包共22个文件大小约8.08兆字节除核心的矢量图形、索引、属性表、投影信息等Shapefile完整组件外还附带空间索引、坐标系说明及路网预览图属性字段包含道路名称、道路等级、行车方向、车道数量等基础信息。数据采用WGS84或CGCS2000坐标系具体以投影信息文件为准可直接导入ArcGIS、QGIS、SuperMap等平台进行编辑、叠加与空间分析用于交通规划、路径分析、城市建模等场景。已有70人学习下载文件命名统一、结构清晰开箱即用省去自行爬取与清洗路网数据的繁琐步骤。1. 北京路网能干什么一份成熟矢量数据包怎么省掉三天预处理做GIS开发的人大多有这种经历接到一个需求要在城区路网上做可达性分析甲方说“数据你去找”结果下载了三份来源不同的路网一份没投影参数、一份属性字段全是代码、一份路口断成蜘蛛网。北京路网这类道路矢量数据包表面上看就是Shapefile全套文件真正拉开差距的是你能不能把这套文件当成一个完整的数据基础设施来用。这篇文章从文件结构、坐标系校验、属性字段、Python读取到网络分析把一套城区道路数据包的用途、参数和雷区一次拆完新手能照着跑通熟手也能对照校准自己的预处理流程。2. Shapefile底层逻辑读懂四件套文件与投影坐标系再动手2.1 Shapefile不是一个文件是shp/shx/dbf/prj四件套的协同配合很多第一次接触矢量数据包的开发者看到下载目录里同时躺着五六个后缀名不同的文件第一反应是“是不是压缩包损坏了”。实际上Shapefile规范从一开始就是多文件并存的.shp存几何坐标.shx存坐标索引.dbf存属性表.prj存投影信息。四者共用同一个主文件名任意缺失一个数据就可能出现“图层加载失败”或“属性表打不开”的问题。我在拿到一份城区路网数据包后第一步永远是在终端里确认文件完整性而不是直接拖进GIS软件ls -lh /data/beijing_road/ # 预期输出示例 # -rw-r--r-- 1 user user 48M beijing_road.shp # -rw-r--r-- 1 user user 6M beijing_road.shx # -rw-r--r-- 1 user user 32M beijing_road.dbf # -rw-r--r-- 1 user user 1K beijing_road.prj # -rw-r--r-- 1 user user 80 beijing_road.cpg逻辑说明这一步的重点不是看文件大小而是确认四个核心文件都在。.shx缺失时QGIS会尝试自动重建索引但Python的pyshp库会直接报错.prj缺失时几何数据本身还能打开但坐标会完全错乱.dbf缺失则整个属性表不可读。参数说明ls -lh里的-l列出详细属性-h把字节数显示为可读的MB/KB。实际环境中文件体积会根据数据范围浮动但四件套的配套关系是固定的。如果看到.cpg文件也存在说明原始数据在导出时已经指定了字符编码这个后面调字段中文名时会用到。2.2 坐标系是第一步校验CGCS2000 与 WGS84 的经典误解路网数据最容易被忽略、出问题时最难排查的就是坐标系。北京城区的道路数据常见的投影包含两种一种是WGS84经纬度坐标EPSG:4326单位是度另一种是CGCS2000投影坐标例如EPSG:4547或EPSG:4548单位是米。很多做Web开发的人拿到数据后不做任何校验直接把经纬度数据当成投影米制数据做缓冲区分析出来的结果半径会缩小到十分之一以下。.prj文件是一个纯文本文件可以直接查看关键信息cat /data/beijing_road/beijing_road.prj一个典型的CGCS2000 3-degree Gauss-Kruger投影文件内容包含PROJCS[CGCS2000 / 3-degree Gauss-Kruger zone 39, ...]这样的描述。注意关键字如果出现GEOGCS[WGS 84]和UNIT[degree]说明是经纬度如果出现PROJCS和UNIT[metre]说明是投影坐标。逻辑说明同一份道路矢量数据用投影坐标系显示的图层范围是几万到几十万的量级米用经纬度坐标系显示的图层范围是115到117度。如果预览时发现道路横跨整个画布或者缩放到看不到第一反应应该是去查.prj而不是怀疑数据损坏。参数说明EPSG是欧洲石油调查组织的坐标系编码4326代表WGS84经纬度4547/4548这类编码代表CGCS2000下的分带投影。城区路网一般用高斯-克吕格3度分带投影37带覆盖东经111-11438带覆盖114-117北京多数区域在39带覆盖东经117-120附近具体用几带要按数据范围的中经线判断。我不太建议按行政区划猜带号直接在GIS软件里看图层范围更保险。2.3 属性表里藏着哪些可用字段从道路等级到车道数一份能直接用来做分析的城区路网属性表不是随便填的。常见的字段组合包括道路唯一标识、道路名称、道路等级高速/主干/次干/支路、车道数、是否单行道、限速、长度等。读取属性表最快的方式是用ogrinfo命令ogrinfo -so -al /data/beijing_road/beijing_road.shp输出中会列出所有字段名和类型。典型字段如下字段名类型含义常见取值road_idInteger道路唯一编号10001, 10002...nameString道路名称某环路、某大街classInteger道路等级1高速, 2主干, 3次干, 4支路lanesInteger车道数2, 4, 6, 8onewayInteger单行道标志0双向, 1正向, -1反向speedInteger限速(km/h)40, 60, 80, 120lengthDouble线长度(米)随投影单位变化逻辑说明这些字段是否齐全直接决定了后续网络分析能不能做。没有speed字段路径规划就得靠等级字段去映射默认速度没有oneway字段图模型就只能按无向图处理没有lanes通行能力估算就成了拍脑袋。参数说明class字段的取值约定在不同数据源里并不统一有的用0到3有的用1到5。拿到数据后先打印唯一的取值集合确认语义映射这一步永远是先于分析的。length字段如果是从ArcGIS导出的可能已有预计算长度如果是手工编辑的很可能全是0那就需要自己用几何计算重新生成。3. 数据接入实操从QGIS预检到Python批量读取的完整路径3.1 用QGIS完成图层加载与视觉预检QGIS是验证数据完整性的首选工具加载逻辑很简单菜单栏选择“图层 → 添加图层 → 添加矢量图层”然后选中.shp文件。加载完成后第一步是查看图层面板的属性确认坐标系和范围显示是否合理。凭经验快速预检打开图层属性切到“信息”页看“范围”一栏的数值。如果是经纬度坐标系WGS84范围应该是北纬39-41、东经115-117这样的度数值如果是投影坐标系范围会是六到七位数的大数值米。两项预期不符调整方式是在“图层 → 图层属性 → 源”里手动指定正确的坐标系。视觉预检还有一个习惯把样式里的线宽调细用单一颜色显示然后缩小到图斑全貌快速扫一眼有没有孤立的短线头、有没有明显的格网断裂。这个动作在预处理阶段大概花五分钟但能排除大约三成的后续问题。3.2 用GeoPandas完成路网读取与基础统计命令行验证完成了接下来进入Python脚本阶段。GeoPandas是读取和操作Shapefile最顺手的库一条read_file就能把整个矢量图层读成GeoDataFrameimport geopandas as gpd gdf gpd.read_file(/data/beijing_road/beijing_road.shp) print(gdf.crs) # 查看坐标系 print(gdf.shape) # 查看行数和列数 print(gdf.total_bounds) # 查看四至范围 gdf[length_m] gdf.geometry.length print(f路网总长度: {gdf[length_m].sum() / 1000:.1f} km) print(gdf[class].value_counts())逻辑说明用total_bounds直接拿四至判断范围是否落在北京城区坐标区间经纬度坐标系预期在[115.5, 39.5, 117.0, 41.0]附近投影坐标系则要看具体带号如果四至包含(0, 0)这种明显异常点极可能是图层中混入了Null几何或坐标为0的线要素。参数说明gdf.geometry.length返回的是投影单位长度如果CRS是WGS84经纬度计算出的“长度”单位是度直接求和再转公里是错的。正确做法是先gdf gdf.to_crs(epsg4547)把坐标系转换为米制投影再做几何计算。class的取值统计是为了后面速度映射用的建议此时顺便打印完整取值不要想当然认为只有四类。3.3 导出为GeoJSON供Web工具使用的格式细节路网数据最终经常要供Web前端或轻量分析工具使用GeoJSON是通用格式导出过程本身简单但有几个参数值得注意# 先转换到WGS84经纬度再导出前端地图默认用4326 gdf_wgs84 gdf.to_crs(epsg4326) # 保留需要字段缩小文件体积 cols [road_id, name, class, lanes, oneway, speed] gdf_export gdf_wgs84[cols [geometry]] # 导出GeoJSON gdf_export.to_file( /data/beijing_road_4326.geojson, driverGeoJSON, encodingutf-8 )逻辑说明整个代码块表达的工作流是“先转坐标系再裁剪字段最后导出”。两个常见坑分别是忘记转换坐标系直接导出导致前端地图无法对齐底图把所有字段全部导出文件从几MB膨胀到几十MB加载变慢。参数说明to_crs(epsg4326)是GeoPandas重投影的标准写法。to_file里的driverGeoJSON必须显式声明否则GeoPandas默认会按文件后缀推断不至于出错但养成显式声明习惯能避免边界情况翻车。encodingutf-8是为了保证中文路名在GeoJSON里正常显示缺失时通常出现乱码。4. 道路数据包避坑指南坐标系偏差、拓扑断裂与字段编码4.1 图层打开后一片空白缩放到全图才看到几条歪斜的线现象QGIS或ArcGIS里加载图层后地图主界面没有任何内容按“缩放到图层范围”后屏幕边缘出现几条或十几条乱线不像道路网。原因八成是坐标系信息丢失或错误。数据源在导出时没有写入.prj软件默认按WGS84经纬度读取但实际几何坐标是CGCS2000投影米制数值两者数值量级差约十万倍导致图面严重缩放的线根本看不到。解决右键图层打开属性在“源”标签里手动指定正确的CRS。如果知道数据原本是CGCS2000 3度带投影直接选择对应EPSG如果不确定可以先读取一个线的坐标值和道路的真实经纬度做对比反推投影带。从那以后我每次拿到新路网第一件事永远是打印图层范围而不是直接看地图。4.2 属性表中文路名全部变成乱码现象name字段里的中文道路名显示成“鐜嬪簻搴?”在QGIS里看全是问号在Python里设置encodingutf-8读取依然乱码。原因.dbf是dBASE格式属性表的字符编码在导出时被写成了GBK或ANSI而非UTF-8。QGIS加载时默认读取.cpg文件中的编码声明但很多工具导出时不会自动写入正确的.cpg导致软件按默认UTF-8解码GBK字节流出现乱码。解决在QGIS图层属性里把“编码”改为“GBK”或“System”通常能恢复中文在Python里读取时显式指定encoding参数import geopandas as gpd # 尝试用GBK读取属性表规避UTF-8解码乱码 gdf gpd.read_file( /data/beijing_road/beijing_road.shp, encodinggbk ) print(gdf[name].head())逻辑说明代码块里传的关键参数是encodinggbk它会直接告诉底层fiona库用GBK解码.dbf文件中的字符串字段。处理思路是先小范围打印前几行确认路名正常再继续后续分析不要一次性全量处理乱码数据输出的结果会连带污染导出的GeoJSON。参数说明encoding可选值包括utf-8、gbk、gb2312、latin1。对于国内分享的城区路网数据GBK是最常见来源编码如果gbk还乱码就再试gb2312注意两者互不兼容的情况较少。用print预览前几行是避免大批量导出后才发现问题的后悔药。4.3 路口处道路没有相交网络分析时路径断续现象把路网数据导入networkx做最短路径分析时发现原本一个十字路口被拆成了四条互不相连的线段车流量无法跨越路口。原因数据在原始数字化时不要求路口处节点完全重合或多年编辑后线段的顶点发生了微小偏移。Shapefile本身不保证拓扑一致性线要素的端点不严格等于另一个线要素的端点图论建模时就形成了断点。解决在QGIS里可以手动画点连接批量的做法是使用PostGIS的ST_Snap或Python里的shapely缓冲求交。常见做法是构建一个容差范围把线段端点在容差内的坐标对齐import geopandas as gpd from shapely.ops import snap gdf gpd.read_file(/data/beijing_road/beijing_road.shp) gdf_clean gdf.copy() # 容差设为10米通常足够对齐野外采集误差 tolerance 10 geometry [snap(line, gdf.geometry.unary_union, tolerance) for line in gdf.geometry] gdf_clean.geometry geometry逻辑说明snap的作用是把每条线段吸到容差范围内最近的其他几何对象上实现端点对齐。这个操作可能轻微改变道路几何形状所以容差不能设太大10米在城市管线测量里常见但如果你依赖线端点精确位置做后续计算这一步骤需要配合人工检查。参数说明这里tolerance的单位会随数据坐标系变化。如果数据是经纬度坐标系10代表的是10度显然不合理通常要先转到投影坐标系再执行对齐执行后原地验证节点连通性在networkx里统计孤立子图数量子图数量大幅减少说明对齐生效。4.4 道路等级字段语义与预期不一致速度映射出错现象按官方字段说明把class为1的当高速路但统计出的高速里程不超过10公里明显与常识不符。原因国内路网数据的分类方法很杂有按行政等级国道/省道/县道分的有按通行等级高速/主干/次干/支路分的字段都叫class但取值代表的意义完全不同。解决手工交叉验证。先打印class的取值和道路名称对照路名判断真实等级再用一条已知高速和一条已知支路做锚点反推整张表的等级映射。写一套映射字典import geopandas as gpd gdf gpd.read_file(/data/beijing_road/beijing_road.shp) # 先观察实际取值 print(gdf[class].value_counts()) # 将原始等级映射为通行速度 class_speed { 1: 110, 2: 60, 3: 40, 4: 30 } gdf[speed] gdf[class].map(class_speed).fillna(30) print(gdf[speed].value_counts())逻辑说明代码里的映射字典就是避坑的核心它是把不可信的原始字段翻译为可用于计算的通行速度。不做这一步后续所有基于时间的分析全都会偏差做了之后再抽样几条路核对速度值是否合理抽样的路名可以在属性表里按名称筛选。参数说明映射值可以按实际需求调整比如高速默认按110 km/h但如果你想做城区高峰时段分析可能要整体调低到80或60。fillna(30)用于兜底确保个别未知等级的线不会变成0速度导致路径规划绕过整队道路。5. 从路网到网络分析构建城市道路图与等时圈计算5.1 把Shapefile转成networkx图模型的核心操作把矢量路网接入路径规划需要先把线要素集合转成图数据结构。networkx配合shapely就能完成这个转换把每条线段拆成起点和终点两个节点线段本身作为边道路长度或通行时间作为权重。import networkx as nx import geopandas as gpd from shapely.geometry import LineString from itertools import permutations gdf gpd.read_file(/data/beijing_road/beijing_road.shp) gdf gdf.to_crs(epsg4547) # 统一为米制投影 G nx.Graph() for idx, row in gdf.iterrows(): geom row.geometry if not isinstance(geom, LineString): continue coords list(geom.coords) prev_node coords[0] prev_id f{idx}_0 G.add_node(prev_id, posprev_node) for j, cur in enumerate(coords[1:], start1): cur_id f{idx}_{j} G.add_node(cur_id, poscur) seg_len LineString([prev_node, cur]).length G.add_edge(prev_id, cur_id, weightseg_len) prev_node, prev_id cur, cur_id逻辑说明为了让道路中线上的每个转折点都成为图节点这段代码把一条多段线拆成多个折线段每个折线段作为一条边录入图里。这个做法的好处是转弯点保留完整几何精度缺点是图的节点数会膨胀但城市范围内一般不会卡死替换方案是把整条道路的端点作为节点、整条线的长度作为边权重节点量更小但转弯细节丢失。参数说明to_crs(epsg4547)的作用是把所有几何统一为米制投影这样LineString.length的返回值直接就是米的数值无需再做比例转换。weightseg_len是以长度为边权后续如果要做时间成本分析还需要除以车速换算。5.2 道路等级映射速度把长度权重换算成时间权重上一章已经建好了speed字段本节的意义是把它转为时间成本。通行时间等于长度除以速度得到秒数import networkx as nx # 假设G已经由5.1节代码构建完成 # 每条边增加travel_time成本单位是秒 speed_col gdf.set_index(gdf.index)[speed] # 简化示意 for u, v, data in G.edges(dataTrue): speed 30 # 默认值 for idx, row in gdf.iterrows(): if f{idx}_0 in (u, v): speed row[speed] break length_m data[weight] data[travel_time] length_m / (speed / 3.6)逻辑说明最直接的操作是遍历现有图的每一条边从GeoDataFrame里找出边归属的道路取其限速值再把长度除以速度换算成秒。这里的(speed / 3.6)是把km/h换算成m/s。参数说明遍历gdf逐行匹配是低效操作数据量大会慢。实际工程中我一般先在GeoDataFrame里给每条道路一个唯一ID然后把ID写进图的节点名里这样读边时直接从节点名里解析道路ID完全不需要遍历属性表。示例里保留读者偏好的全遍历写法是为了让你先跑通流程。5.3 最短路径与等时圈的近似实现有了带travel_time的图最短路计算就是一行代码的事import networkx as nx origin 10001_0 # 起点的节点ID取自道路10001的起点 target 25088_3 # 终点的节点ID取自道路25088的某个转折点 route nx.shortest_path(G, origin, target, weighttravel_time) route_time nx.shortest_path_length(G, origin, target, weighttravel_time) print(f预计行驶时间 {route_time / 60:.1f} 分钟) print(f途经节点数 {len(route)})逻辑说明从人工选定的起点和终点节点出发nx.shortest_path_length算累计时间shortest_path返回经过的节点序列。真实应用中起终点往往是坐标点需要先投影到最近的图节点改造方法是shapely求最近点投影到路线上再接入图。参数说明weighttravel_time必须与5.2节里设置的属性名一致否则算法默认按无权图计算最短路径通常表现为“距离最短但不走大路”。等时圈近似可用nx.single_source_dijkstra_path_length(G, origin, cutoff1800)3600代表1小时、1800代表30分钟直接得到一个源点到其他节点的可达节点集合再把这些节点围起来做凸包就是粗略的时间圈层。分析类型核心函数参数要点最短路径nx.shortest_path指定weight为travel_time多源最短距离nx.single_source_dijkstra_path_length指定cutoff控制时间阈值路径还原nx.shortest_path保留pos属性便于前端划线6. 最后的细节一套数据的真实边界与我最常做的三项校验一套路网数据包拿到手前面五章讲的是流程这里补上三个我在交付前必做的校验动作也是这套数据真正的边界所在。第一是坐标质量抽查。随机在数据里选五条环路把线段的形状和路名与实际路网底图叠加看有没有走样。路网数据大概率来自公开采集或众包编辑环路的形状偏移常常在视觉上不易发现但叠加高清底图一看就能看出偏差。城区级数据精度大约在10米左右是常态米级精度几乎不可能是免费分享资源里能见到的。第二是属性完整度评估。打印每个字段的空值比例。name字段缺失率超过百分之三十的路网做地名分析基本是亏的但如果只是用来做路网拓扑、几何修复教学name字段缺失不影响使用。lanes和speed字段如果有连续的空值网络分析的可靠性会下降我通常用等级映射值来补空缺同时会在交付文档里注明“含估算字段”。第三是路网拓扑的一次最终验证。把整张图传入networkx统计孤立子图的数量。一个城区完整路网通常只有一个连通主图如果有几十个孤立的小子图大概率是路口断点没处理干净。我会再跑一次snap流程或者用nx.connected_components找出断路的道路簇输出到调查文件中。这三项校验做完这份数据的能做什么、不能做什么基本心里有数。从那以后我每次拿到新的路网数据不管来源多正规都强制走一遍完整流程读.prj确认坐标系打印字段取值做一次拓扑连通性统计全部通过才写入正式分析环境。这个过程已经变成了我的肌肉记忆也帮我拦下了很多本来要在交付后半途返工的麻烦。希望这些拆解和踩坑记录对你处理手头的城区路网数据有实际帮助。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
