赞
踩
前言:本人是在学习交通数据小旭学长写的Python出租车GPS数据的路网匹配(TransBigData+leuvenmapmatching)过程中发现osmnx库由于受到网络限制很难较为完整的获取一个城市或更大范围的路网数据。为了解决这个问题,本人采用了osm2gmns库,构建网络并进行地图匹配,实现了和交通数据小旭学长一样的效果。
现状osm路网数据获取方法较多,但是对于指定城市或者省份数据较难获取,一般要么太大(国家),要么太小(城区范围),因此本人采用了先采集较大范围数据,再进行裁剪的方法进行处理。具体处理方法如下:
登录OSM网站,选择[Geofabrik下载](定期更新的洲、 国家和特定城市的摘录)选择你要选择的地区,选择.osm.pbf下载。
一般这种数据很容易找到,推荐一个,
输入.shp文件,用如下代码生成深圳市的边界文件.poly格式的,然后保证与osmconvert64-0.8.8p.exe在同一文件夹内。
import geopandas as gpd import os from shapely.geometry import Polygon, MultiPolygon def write_poly(polygons, output_file): with open(output_file, "w") as f: c = 1 for polygon in polygons: if isinstance(polygon, Polygon): poly_ext = list(polygon.exterior.coords) f.write("{}\n".format(c)) f.write("1\n") for coord in poly_ext: f.write("\t{} {}\n".format(coord[0], coord[1])) f.write("END\n") elif isinstance(polygon, MultiPolygon): for poly in polygon.geoms: poly_ext = list(poly.exterior.coords) f.write("{}\n".format(c)) f.write("1\n") for coord in poly_ext: f.write("\t{} {}\n".format(coord[0], coord[1])) f.write("END\n") c += 1 f.write("END\n") shp_file = "深圳市.shp" output_file = "shenzhen_polyfile.poly" gdf = gpd.read_file(shp_file) crs_4326 = {"init": "epsg:4326"} gdf = gdf.to_crs(crs_4326) # 将坐标系统转换为WGS84 (EPSG:4326)如果需要 polygons = gdf.geometry write_poly(polygons, output_file)
首先获取osmconvert软件,获取方法很多,推荐1个然后运行osmconvert,进入如下界面
然后输入a,将之前下载好的.osm.pbf文件名称输入(注意:文件也要与软件在同一文件夹内)
回车输入后,选择4,再将刚才准备好的边界区域.poly文件输入
回车后就可以得到处理好的只包含深圳范围的.osm.pbf文件
在这里选择的是osm2gmns库,它可以将osm路网进行选择性清洗,选择你想要的类型的路并构建为路网,具体参数大家可以去库的官网了解,他会输出一个link.csv和node.csv,这是构建路网的关键
import osm2gmns as og
net = og.getNetFromFile('shenzhen-latest.osm_01.pbf', network_types='auto',link_types={'motorway', 'trunk', 'primary', 'secondary', 'tertiary', 'residential', 'service', 'unclassified', 'connector'},strict_mode=False,combine=True,)
og.consolidateComplexIntersections(net, auto_identify=True)
og.outputNetToCSV(net)
import pandas as pd import networkx as nx node_data = pd.read_csv('node.csv', encoding='gbk') link_data = pd.read_csv('link.csv', encoding='gbk') G = nx.Graph() G.graph["crs"] = "epsg:4326" for idx, row in node_data.iterrows(): node_id = row['node_id'] node_lon = row['x_coord'] node_lat = row['y_coord'] G.add_node(node_id, lon=node_lon, lat=node_lat) for idx, row in link_data.iterrows(): link_id = row['link_id'] link_from = row['from_node_id'] link_to = row['to_node_id'] G.add_edge(link_from, link_to, id=link_id)
这里的数据是从交通数据小旭学长写的Python出租车GPS数据的路网匹配(TransBigData+leuvenmapmatching)里面分享的链接中获得的,大家可以自行下载。
import transbigdata as tbd
import pandas as pd
# 地图匹配包
from leuvenmapmatching.matcher.distance import DistanceMatcher
from leuvenmapmatching.map.inmem import InMemMap
from leuvenmapmatching import visualization as mmviz
#读取数据
data = pd.read_csv('TaxiData-Sample.csv',header = None)
data.columns = ['VehicleNum','Time','Lng','Lat','OpenStatus','Speed']
#从GPS数据提取OD与路径GPS
oddata = tbd.taxigps_to_od(data,col = ['VehicleNum','Time','Lng','Lat','OpenStatus'])
data_deliver,data_idle = tbd.taxigps_traj_point(data,oddata,col=['VehicleNum', 'Time', 'Lng', 'Lat', 'OpenStatus'])
from leuvenmapmatching.map.inmem import InMemMap # 创建一个 InMemMap 对象 map_con = InMemMap(G, use_latlon=True, use_rtree=True, index_edges=True) # 遍历 link_data 添加节点和边 for _, row in link_data.iterrows(): link_id = row['link_id'] link_from = row['from_node_id'] link_to = row['to_node_id'] # 提取节点的经纬度 lat1, lon1 = G.nodes[link_from]['lat'], G.nodes[link_from]['lon'] lat2, lon2 = G.nodes[link_to]['lat'], G.nodes[link_to]['lon'] # 在 map_con 中添加节点 map_con.add_node(link_from, (lat1, lon1)) map_con.add_node(link_to, (lat2, lon2)) # 连接这两个节点 map_con.add_edge(link_from, link_to) # 清除不必要的信息以提高性能 map_con.purge()
#用transbigdata提取出行轨迹 import geopandas as gpd tmp_gdf = data_deliver[data_deliver['ID'] == 22].sort_values(by = 'Time') #轨迹增密 tmp_gdf = tbd.traj_densify(tmp_gdf,col = ['ID', 'Time', 'Lng', 'Lat'],timegap = 15) #转换轨迹的坐标系为地理坐标系 tmp_gdf['geometry'] = gpd.points_from_xy(tmp_gdf['Lng'],tmp_gdf['Lat']) tmp_gdf = gpd.GeoDataFrame(tmp_gdf) tmp_gdf.crs = {'init':'epsg:4326'} #tmp_gdf = tmp_gdf.to_crs(2416) #获得轨迹点 path = list(zip(tmp_gdf.geometry.y, tmp_gdf.geometry.x)) #构建地图匹配工具 matcher = DistanceMatcher(map_con, max_dist=500, max_dist_init=170, min_prob_norm=0.0001, non_emitting_length_factor=0.95, obs_noise=50, obs_noise_ne=50, dist_noise=50, max_lattice_width=20, non_emitting_states=True) #进行地图匹配 states, _ = matcher.match(path, unique=False) #绘制底图匹配结果 mmviz.plot_map(map_con, matcher=matcher, show_labels=False, show_matching=True,#show_graph=False, filename=None)
获得如下结果图
注意:该方法通过matcher.path_pred_onlynodes获得匹配的每个路段节点编号,命名为u,由于路段是前后相接的因此错位合并后得到每个路段的起始和终点路段节点编号
#该方法通过matcher.path_pred_onlynodes获得匹配的每个路段节点编号,命名为u,由于路段是前后相接的因此错位合并后得到每个路段的起始和终点路段节点编号
pathdf = pd.DataFrame(matcher.path_pred_onlynodes,columns = ['u'])
pathdf['v'] = pathdf['u'].shift(-1)
pathdf = pathdf[-pathdf['v'].isnull()]
#使用pd.merge与link表的起始终止节点进行匹配,从而得到匹配路段的地理信息
from shapely import wkt
# 需要确保CSV文件中的"geometry"列已经以合适的格式存储了几何数据(例如,以"Well-Known Text (WKT)格式存储)
link_data['geometry'] = link_data['geometry'].apply(wkt.loads)
pathgdf = pd.merge(pathdf,link_data.reset_index(),left_on=['u','v'],right_on=['from_node_id','to_node_id'])
pathgdf = gpd.GeoDataFrame(pathgdf,geometry='geometry')#生成Geopandas文件
pathgdf.plot()
运行结果如下为匹配的路段
首先要获得所有路网路段的中心点,方便后面筛选要显示的路段
#link_data['geometry'] = link_data['geometry'].apply(wkt.loads) #获得所有线段中心点 from shapely.geometry import Point # 创建空的列表以存储经度和纬度值 lons = [] lats = [] # 对于表中的每一行 for geometry in link_data['geometry']: # 计算中心点 centroid = geometry.centroid # 将经度和纬度值添加到相应的列表中 lons.append(centroid.x) lats.append(centroid.y) # 将经度和纬度值添加到表中 link_data['lon'] = lons link_data['lat'] = lats
然后进行轨迹匹配可视化
#与路网一起可视化 import matplotlib as mpl import matplotlib.pyplot as plt fig = plt.figure(1,(8,8),dpi = 100) ax = plt.subplot(111) plt.sca(ax) fig.tight_layout(rect = (0.05,0.1,1,0.9)) #设定可视化边界 bounds = pathgdf.unary_union.bounds gap = 0.003 bounds = [bounds[0]-gap,bounds[1]-gap,bounds[2]+gap,bounds[3]+gap] #绘制匹配的路径 pathgdf.plot(ax = ax,zorder = 1) #绘制底图路网 link_data = gpd.GeoDataFrame(link_data,geometry='geometry') tbd.clean_outofbounds(link_data,bounds,col = ['lon','lat']).plot(ax = ax,color = '#333',lw = 0.1) #绘制GPS点 tmp_gdf.to_crs(4326).plot(ax = ax,color = 'r',markersize = 5,zorder = 2) plt.axis('off') plt.xlim(bounds[0],bounds[2]) plt.ylim(bounds[1],bounds[3]) plt.show()
最后,对前人的无私分享表示感谢!
Copyright © 2003-2013 www.wpsshop.cn 版权所有,并保留所有权利。