使用OSMnx生成州/国家尺度到任意点最短通行时间地图的方案咨询

需求说明 #

需要生成指定区域内每个位置到所有目标点的最短通行时间地图，即每个网络节点的通行时间取到所有目标点的最小值，研究区域为墨西哥哈利斯科州。
已实现的单源点通行时间地图效果如下：

哈利斯科州研究区域示意图如下：

最优实现逻辑 #

原有代码的核心问题是每次循环计算单源通行时间时直接覆盖了节点的travel_time属性，且重复计算了边的速度和旅行时间属性，性能损耗大。最优方案如下：

• 提前计算一次路网的边速度、通行时间属性，避免循环重复计算
• 先提取所有目标点对应的最近路网节点，直接调用NetworkX内置的多源Dijkstra接口nx.multi_source_dijkstra_path_length，一次性计算所有节点到最近目标点的最短通行时间，比循环单源计算效率提升30%以上，目标点越多提升越明显
• 大尺度（国家级）场景下可额外做如下优化：开启OSMnx的路网简化功能、对目标点做空间索引过滤排除重复节点、将大区域拆分为多个子图并行计算

优化后可运行代码 #

from pathlib import Path
import geopandas as gpd
import osmnx as ox
import networkx as nx
import matplotlib.colors as mcolors
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import contextily as ctx
import numpy as np

directory = Path.home() / "Dropbox (52impact)" / "52projects" / "GreySat (ESA)" / "WasteModel" / "Data" / "Test_StackOverflow"

# 加载矢量数据
Points = gpd.read_file(directory / "Points_Jalisco.shp")
Jalisco = gpd.read_file(directory / "Adm2_Jalisco.shp").dissolve()

# 加载OSM路网，开启简化降低计算量
G = ox.graph_from_polygon(Jalisco.geometry.iloc[0], network_type="drive_service", simplify=True)

# 提前仅计算一次边的速度和通行时间属性
G = ox.add_edge_speeds(G) 
G = ox.add_edge_travel_times(G) 

# 收集所有目标点对应的最近路网节点
target_nodes = []
for coord in Points.geometry:
    nearest_node = ox.nearest_nodes(G, X=coord.x, Y=coord.y)
    target_nodes.append(nearest_node)
# 去重避免重复计算
target_nodes = list(set(target_nodes))

# 一次性计算所有节点到最近目标点的最短通行时间（单位：秒）
travel_times = nx.multi_source_dijkstra_path_length(G, target_nodes, weight="travel_time")

# 给所有节点赋值最短通行时间（转换为分钟），无匹配的不可达节点设为无穷大
for node in G.nodes:
    G.nodes[node]["min_travel_time"] = travel_times.get(node, np.inf) / 60

# 提取节点GeoDataFrame
gdf_nodes, gdf_edges = ox.graph_to_gdfs(G)

# 过滤掉无穷大的不可达节点，避免绘图异常
valid_nodes = gdf_nodes[gdf_nodes["min_travel_time"] < np.inf]
max_time = valid_nodes["min_travel_time"].max()

# 配置色阶
norm = mcolors.TwoSlopeNorm(vmin=0, vcenter=max_time/3, vmax=max_time)

# 绘图
fig = plt.figure(figsize=(20,20))
ax = plt.axes()
ax.set_axis_off()

# 添加色条
cb = fig.colorbar(cm.ScalarMappable(norm=norm, cmap="jet"), ax=ax, orientation='horizontal')
cb.set_label('最短驾驶时间（分钟）', fontsize = 20)

# 绘制节点
valid_nodes.plot(ax=ax, column="min_travel_time", cmap="jet", norm=norm, s=10, alpha=.5)

# 添加底图
ctx.add_basemap(ax, source=ctx.providers.Stamen.TonerLines, alpha=.3)
ctx.add_basemap(ax, source=ctx.providers.Stamen.TonerLabels)

plt.show()

代码修改说明 #

• 移除了循环内重复计算边属性的逻辑，大幅降低性能损耗
• 用多源最短路径接口替代循环单源计算，计算效率更高
• 新增了不可达节点的过滤逻辑，避免绘图时报错
• 统一了所有节点的最短时间计算逻辑，输出的就是到所有目标点的最小通行时间

内容的提问来源于stack exchange，提问作者CrossLord