ARTICLE · INTELLIGENCE

战地情报 · 详情页

来自尧图项目组的一线实战观察与深度解析

郑州OSM道路矢量数据:从拓扑修复到路网分析实战

郑州OSM道路矢量数据:从拓扑修复到路网分析实战 简介这份资源是面向GIS从业者、城市规划与交通研究人员以及相关专业学生的郑州市道路矢量数据集基于OpenStreetMap志愿者维护的开源地理信息整理加工而成可直接用于空间分析与专题制图。压缩包共9个文件约4.49MB以Shapefile格式为核心包含存储道路线几何的shp、记录道路类型与名称等属性的dbf、定义坐标参考系的prj以及加速读取的shx、sbn、sbx索引和指定中文编码的cpg文件另附一份OSM道路类别对照表图片便于正确识别道路等级与标签含义。目前已有395人学习下载。借助QGIS、ArcGIS等软件加载后可开展交通热点识别、最短路径计算并与人口、建筑、公交线路等数据叠加服务于城市规划、交通管理与应急响应等场景是一份上手即用的基础地理数据。1. 郑州市 OSM 道路矢量数据一份能直接进 GIS 的成品数据拿到一份城市路网数据最怕的不是数据量不够而是拿到手发现字段乱、拓扑断、坐标系对不上光清洗就得搭进去两天。这份郑州市 OSM 道路矢量数据已处理解决的就是这个环节——它把 OpenStreetMap 原始道路数据里那些让人头疼的毛刺提前处理掉了拿到之后可以直接拖进 QGIS、ArcGIS 或者用 GeoPandas 读出来做空间分析。适合做城市交通研究、路网密度计算、可达性分析、路径规划验证以及需要底图路网做可视化叠加的从业者。如果你之前下过 OSM 的 shp 包打开发现属性表里全是英文 tag、道路断成一段一段、投影还是 WGS84 经纬度那这份数据的价值就很具体了省掉从 raw data 到可用图层之间的那一段脏活。2. 先搞清楚 OSM 道路数据的原始结构为什么不能直接拿来用2.1 OSM 的数据模型和道路要素的对应关系OpenStreetMap 的数据组织方式和传统 GIS 矢量文件有本质区别。它底层是 XML 格式由 node节点、way路径、relation关系三种元素构成。一条道路在 OSM 里是一个 way由若干 node 串联而成每个 node 带经纬度坐标。道路的属性不放在固定字段里而是以 key-value 形式的 tag 挂在 way 上比如highwayprimary、name中原路、onewayyes、maxspeed60。这种灵活的结构对社区编辑友好但对做空间分析的人不友好。因为 tag 是开放的同一种道路可能被标成highwayprimary也可能被标成highwaytrunk加上trunk:category1。更麻烦的是OSM 的 way 只表示一段几何两条道路交叉处如果没有共享 node拓扑上就是断开的。直接拿原始数据算路网连通性结果会严重偏低。常见做法是先用 osmium 或 osm2pgsql 把.osm.pbf转成 GIS 可读格式再做 tag 映射和拓扑修复。这份已处理数据相当于把这几步前置完成了。2.2 已处理到底处理了什么从文件名和常见处理流程推断这份数据的处理链路大致包含以下几个环节坐标系转换。OSM 原生是 EPSG:4326WGS84 经纬度做距离计算和面积统计需要投影坐标系。郑州位于东经 112°42′ 到 114°14′北纬 34°16′ 到 34°58′对应的 CGCS2000 3 度带投影是 EPSG:4547中央经线 114°E。如果这份数据已经转到了投影坐标系那算路网密度时就不用自己再转一步。道路等级归一化。OSM 的 highway 值有二十多种从 motorway 到 footway 到 steps。做机动车路网分析时通常只保留 motorway、trunk、primary、secondary、tertiary 及其 link外加 residential、service 等。已处理数据一般会把原始 tag 映射成一个简洁的等级字段比如road_class或fclass。拓扑修复。在交叉口处打断道路确保共享节点这样网络数据集才能正确构建连通关系。QGIS 里的「打断相交线」、PostGIS 的ST_Node函数都是干这个的。字段精简。原始 OSM 的 tag 可能有几十个处理后会保留 name、road_class、oneway、maxspeed、bridge、tunnel 等对分析真正有用的字段。注意不同来源的「已处理」标准不一样。拿到数据后第一件事是用 QGIS 打开属性表确认字段名和坐标系别假设它一定符合你的预期。3. 把数据用起来从加载到路网分析的完整操作链3.1 在 QGIS 中加载与坐标系确认假设你拿到的是 shapefile 或 GeoPackage 格式第一步是确认坐标系。在 QGIS 里右键图层 → 属性 → 信息看 CRS 是什么。如果是 EPSG:4326做缓冲区分析前需要重投影。# 用 ogr2ogr 查看数据的基本信息和坐标系 ogrinfo -so -al zhengzhou_roads.gpkg # 如果坐标系是 4326重投影到 CGCS2000 3度带EPSG:4547 ogr2ogr -f GPKG zhengzhou_roads_projected.gpkg zhengzhou_roads.gpkg \ -t_srs EPSG:4547 -nln roads_projogrinfo -so -al会输出图层名、几何类型、要素数量、字段列表和坐标系。ogr2ogr的-t_srs指定目标坐标系-nln指定新图层名。重投影后长度和面积计算的单位就从度变成了米结果才有物理意义。如果数据已经是投影坐标系跳过这步。判断方法很简单看坐标值。经纬度的经度范围在 112 左右投影坐标的 X 值通常在 500000 上下3 度带带号 38false easting 500000。3.2 用 GeoPandas 做路网密度统计路网密度是城市交通分析的基础指标单位是 km/km²。计算逻辑是在目标区域内统计道路总长度除以区域面积。下面用 GeoPandas 走一遍完整流程。import geopandas as gpd import pandas as pd # 读取道路数据和区域边界 roads gpd.read_file(zhengzhou_roads_projected.gpkg, layerroads_proj) districts gpd.read_file(zhengzhou_districts.gpkg, layerdistricts) # 确认坐标系一致 assert roads.crs districts.crs, 坐标系不一致先统一 # 按道路等级筛选机动车道路 motor_roads roads[roads[road_class].isin( [motorway, trunk, primary, secondary, tertiary] )].copy() # 计算每条道路的长度投影坐标系下单位为米 motor_roads[length_m] motor_roads.geometry.length # 空间连接把道路按所在区域分组 joined gpd.sjoin(motor_roads, districts[[district_name, geometry]], howinner, predicateintersects) # 按区域汇总道路长度公里 length_by_district joined.groupby(district_name)[length_m].sum() / 1000 # 计算各区域面积平方公里 districts[area_km2] districts.geometry.area / 1e6 # 合并计算路网密度 result districts[[district_name, area_km2]].merge( length_by_district.rename(road_km), ondistrict_name) result[density_km_per_km2] result[road_km] / result[area_km2] print(result.sort_values(density_km_per_km2, ascendingFalse))这段代码的关键点有三个。第一geometry.length在投影坐标系下返回米在经纬度下返回度所以前面的重投影不是可选项。第二gpd.sjoin用intersects谓词做空间连接一条道路跨两个区时会被分配到两个区长度会重复计算。如果需要精确到不重复得先按区域边界裁剪道路再统计。第三road_class字段名要按实际数据的字段名替换不同处理流程可能叫fclass、highway或type。3.3 构建网络数据集做路径分析路网分析的核心是构建拓扑网络。在 QGIS 里可以用 GRASS 的v.net模块在 Python 里用 networkx 或 igraph。下面用 networkx 演示从矢量路网到最短路径的完整过程。import networkx as nx import geopandas as gpd from shapely.geometry import Point roads gpd.read_file(zhengzhou_roads_projected.gpkg, layerroads_proj) # 构建无向图有单行道需求时改用 DiGraph G nx.Graph() for idx, row in roads.iterrows(): coords list(row.geometry.coords) # 用首尾节点坐标作为图节点标识 start coords[0] end coords[-1] length row.geometry.length # 添加边权重为长度 G.add_edge(start, end, weightlength, road_namerow.get(name, )) # 找两个最近节点作为起终点 def nearest_node(G, point): return min(G.nodes, keylambda n: (n[0]-point.x)**2 (n[1]-point.y)**2) origin nearest_node(G, Point(760000, 3850000)) dest nearest_node(G, Point(780000, 3860000)) # Dijkstra 最短路径 path nx.dijkstra_path(G, origin, dest, weightweight) total_length nx.dijkstra_path_length(G, origin, dest, weightweight) print(f路径经过 {len(path)} 个节点总长度 {total_length/1000:.2f} 公里)这里有个容易翻车的地方直接用道路的起终点作为图节点只在道路首尾相接时才能正确连通。如果数据没有在交叉口打断两条交叉道路各走各的图就是断的。判断方法是在 QGIS 里放大看交叉口如果两条路在交叉处没有共享顶点说明拓扑没修。修复方式是用 QGIS 的「处理工具箱 → 矢量几何 → 打断相交线」或者 PostGIS 里ST_Node(ST_Collect(geom))。另一个参数是weight。用几何长度做权重得到的是最短距离路径改成1/maxspeed或通行时间才是最快路径。如果数据里有maxspeed字段可以换算成秒做权重。4. 避坑与排查处理 OSM 道路数据时最容易翻车的五个地方4.1 坐标系搞混导致长度算出来差几十倍现象路网密度算出来是 0.05 km/km²明显偏低或者高到几百明显离谱。原因数据是 EPSG:4326 经纬度但代码里直接用了geometry.length返回的单位是度。一度纬度约 111 公里一度经度在郑州纬度约 92 公里和米的差距是五个数量级。解决计算长度和面积前必须确认 CRS。用gdf.crs查看如果是EPSG:4326先to_crs(epsg4547)再算。养成习惯任何涉及距离的运算前打印一下gdf.crs和第一条几何的长度值心里有个数。4.2 道路等级字段映射不全导致筛选后丢数据现象筛选road_class为 primary、secondary 等之后要素数量比预期少很多有些明显的主干道不见了。原因OSM 原始 tag 里同一条路可能被标为highwayprimary_link匝道或highwaytrunk但带trunk:category1。如果处理脚本只映射了标准值这些变体就被归到 other 或直接丢了。解决先看属性表里road_class有哪些唯一值用roads[road_class].value_counts()打印分布。如果发现大量要素归在 other 或 unknown回溯原始 OSM tag 看是什么值。必要时手动补充映射规则把primary_link归到 primarytrunk_link归到 trunk。4.3 单行道方向信息丢失现象做路径分析时明明有单行道算法却允许逆行结果不合理。原因OSM 里单行道用onewayyes或oneway-1表示-1表示几何方向与通行方向相反。如果处理时只保留了oneway字段但没做方向翻转或者构建图时用了无向图单行道就失效了。解决构建有向图时对onewayyes的边只加正向对oneway-1的边只加反向对onewayno或空值的边加双向。用 networkx 的DiGraph加边时判断方向。4.4 交叉口未打断导致网络不连通现象路径分析报「目标节点不可达」或者连通分量数量远大于预期。原因OSM 原始数据里两条道路交叉但不在交叉处共享 node几何上是两条独立的线。直接拿首尾节点建图交叉口处没有连接。解决在 QGIS 里用「打断相交线」工具或者 PostGIS 里用ST_Node处理。打断后重新构建图连通分量数量应该大幅下降。验证方法nx.number_connected_components(G)如果结果是个位数或几十说明连通性正常如果是几百上千拓扑有问题。4.5 中文路名字段编码问题现象属性表里name字段显示乱码或者导出 CSV 后中文变成问号。原因shapefile 对中文支持差DBF 文件的编码默认是 Latin-1 或系统编码。如果数据是 shapefile 格式且包含中文字段很容易出问题。解决优先用 GeoPackage 格式它对 UTF-8 支持完善。如果必须用 shapefile在 QGIS 里导出时指定编码为 UTF-8或者在 GeoPandas 里读的时候加encodingutf-8参数。已经乱码的数据很难恢复只能重新从原始数据转一遍。5. 进阶用法用这份数据做等时圈分析和可视化等时圈isochrone是交通分析里很实用的一个产出——给定一个点算出 15 分钟、30 分钟能到达的范围。用这份路网数据配合 networkx 就能做不需要调外部 API。思路是这样的从起点出发用 Dijkstra 算法计算到所有节点的最短通行时间然后筛选出时间小于阈值的节点对这些节点做凸包或缓冲区合并得到等时圈多边形。import networkx as nx import geopandas as gpd from shapely.geometry import Point from shapely.ops import unary_union roads gpd.read_file(zhengzhou_roads_projected.gpkg, layerroads_proj) # 构建图权重用通行时间秒 G nx.Graph() for idx, row in roads.iterrows(): coords list(row.geometry.coords) start, end coords[0], coords[-1] length row.geometry.length # 默认速度 40 km/h有 maxspeed 字段则用实际值 speed_kmh float(row[maxspeed]) if row.get(maxspeed) else 40 travel_time length / (speed_kmh * 1000 / 3600) # 秒 G.add_edge(start, end, weighttravel_time) # 起点郑州市中心某点投影坐标 origin min(G.nodes, keylambda n: (n[0]-762000)**2 (n[1]-3852000)**2) # 计算到所有节点的最短时间 times nx.single_source_dijkstra_path_length(G, origin, weightweight) # 提取 15 分钟900秒内可达的节点 reachable_15min [Point(n) for n, t in times.items() if t 900] # 生成等时圈多边形缓冲区半径 100 米模拟道路宽度 buffers [pt.buffer(100) for pt in reachable_15min] isochrone unary_union(buffers) # 导出为 GeoJSON gdf_iso gpd.GeoDataFrame({name: [15min], geometry: [isochrone]}, crsroads.crs) gdf_iso.to_file(isochrone_15min.geojson, driverGeoJSON) print(f15分钟等时圈面积{isochrone.area/1e6:.2f} 平方公里)这段代码有几个可以调的地方。速度参数直接影响等时圈大小40 km/h 是城市道路的保守估计快速路可以设 60-80支路设 20-30。缓冲区半径 100 米是经验值太小会导致多边形碎片化太大则边界模糊。如果数据里有oneway字段把G换成DiGraph并处理方向结果更准确。等时圈做出来之后叠加到 QGIS 里和底图一起看能直观判断某个位置的交通便利程度。做选址分析时把多个候选点的等时圈叠在一起覆盖范围重叠少的那个通常更优。从那以后我每次拿到新的路网数据都会先跑一遍连通分量检查和坐标系确认这两个指标正常了再往下做分析。希望帮到你。本文还有配套的精品资源点击获取
RELATED READING

延伸阅读

更多一线实战笔记与深度复盘,助您持续精进