离线OSM数据条件下，出租车轨迹数据的地图匹配

离线OSM数据条件下，出租车轨迹数据的地图匹配

前言：本人是在学习交通数据小旭学长写的Python出租车GPS数据的路网匹配（TransBigData+leuvenmapmatching）过程中发现osmnx库由于受到网络限制很难较为完整的获取一个城市或更大范围的路网数据。为了解决这个问题，本人采用了osm2gmns库，构建网络并进行地图匹配，实现了和交通数据小旭学长一样的效果。

01 OSM路网数据的获取

现状osm路网数据获取方法较多，但是对于指定城市或者省份数据较难获取，一般要么太大（国家），要么太小（城区范围），因此本人采用了先采集较大范围数据，再进行裁剪的方法进行处理。具体处理方法如下：
登录OSM网站，选择[Geofabrik下载]（定期更新的洲、 国家和特定城市的摘录）选择你要选择的地区，选择.osm.pbf下载。

02 使用osmconvert软件，对数据进行切割

在使用之前先获取所需边界的shp文件，比如深圳市的区域范围。确保是WGS84坐标系的

一般这种数据很容易找到，推荐一个，

然后编写如下代码：构建.poly文件

输入.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对深圳市路网进行提取

首先获取osmconvert软件，获取方法很多，推荐1个然后运行osmconvert，进入如下界面
然后输入a,将之前下载好的.osm.pbf文件名称输入（注意：文件也要与软件在同一文件夹内）
回车输入后，选择4，再将刚才准备好的边界区域.poly文件输入
回车后就可以得到处理好的只包含深圳范围的.osm.pbf文件

03 路网数据的处理

001 首先将.pbf格式的osm路网数据，转化为link和节点

在这里选择的是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)

002 构建图

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)

003 GPS数据读取

这里的数据是从交通数据小旭学长写的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'])

004 构建leuvenmapmatching所需的地图网络（用于地图匹配）

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() 

005 地图匹配

#用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)

获得如下结果图

006 获取地图匹配路径

注意：该方法通过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()

运行结果如下为匹配的路段

007叠加路网可视化

首先要获得所有路网路段的中心点，方便后面筛选要显示的路段

#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()

最后，对前人的无私分享表示感谢！