提示：文章写完后，目录可以自动生成，如何生成可参考右边的帮助文档


文章目录

  • 问题描述
  • 一、核心步骤
  • 1. 生成图
  • 2. 最短路径
  • 二、使用中的问题
  • 三、利用momepy和osmnx生成带编号的MultiGraph
  • 1. 利用momepy生成带编号的edges和nodes
  • 2. 利用nodes和edges生成图
  • 3. 最短路径
  • 总结


问题描述

只有线要素矢量地图，想求任意两点之间的最短路径，尝试用arcgis实现，发现arcgis输出的路径只有长度信息，没有详细的线路信息，于是利用python地理处理包实现最短路径。

一、核心步骤

核心步骤就是两步，首先生成图，再基于图实现最短路径。

1. 生成图

只有线要素矢量地图可以直接用networkx生成图，其中read_shp函数只在2.x版本中有，3.x版本中删除了：

import networkx as nx
shapefile = 'C:/liaoninggaosu/RoadSegment.shp'
graph = nx.read_shp(shapefile)

图的基本元素就是边和节点：

graph.edges

Python打开图片_python

graph.nodes

Python打开图片_结点_02


graph.edges中的((120.02891072, 41.80614275), (120.03331112, 41.80518226))代表一条边，它由两个结点(120.02891072, 41.80614275)和(120.03331112, 41.80518226)组成，结点也记录在graph.nodes中。

2. 最短路径

调用networkx的shortest_path方法，即可实现最短路径，默认使用dijkstra算法，起点和终点需要是graph.nodes中的元素：

nx.shortest_path(G=graph, source=(120.02891072, 41.80614275), target=(120.03544861, 41.80373951), weight='length')

输出如下，输出的是路径经过的节点，可以通过节点再找到边。

[(120.02891072, 41.80614275),
(120.03331112, 41.80518226),
(120.03544861, 41.80373951)]

二、使用中的问题

问题一：networkx.shortest_path的输入要求是图的结点，如果我们的点不在结点处，则需要利用osmnx找到最近的结点，nearest_nodes要求graph为networkx.MultiDiGraph类型，而用networkx.read_shp生成的graph是networkx.DiGraph类型，使用会出报错，KeyError: ‘crs’。

import osmnx as ox
sNode = ox.distance.nearest_nodes(graph, sPoint.x, sPoint.y)

问题二：即使利用networkx.to_pandas_edgelist将图的边转为dataframe后，边与结点的对应关系仍是通过坐标，不方便查找。

df_edge = nx.to_pandas_edgelist(graph)
df_edge.head(3)

Python打开图片_结点_03

三、利用momepy和osmnx生成带编号的MultiGraph

1. 利用momepy生成带编号的edges和nodes

import pandas as pd
import geopandas as gpd
import networkx as nx
import osmnx as ox
import momepy
from shapely.geometry import Point

shapefile = 'C:/liaoninggaosu/RoadSegment.shp'
gdf_map = gpd.read_file(shapefile)
graph = momepy.gdf_to_nx(gdf_map, approach="primal") # graph为MultiGraph类型
nodes, edges, sw = momepy.nx_to_gdf(graph, points=True, lines=True, spatial_weights=True)

生成的nodes和edges如下，edges会自动生成node_start与node_end字段，其与nodes的nodeID字段对应。

Python打开图片_图论_04

2. 利用nodes和edges生成图

然后利用nodes和edges生成graph，需要改成osmnx要求的格式，即节点表索引为osmid、包含x和y字段，边表的多级索引名为u, v, key，详见osmnx - graph_from_gdfs：

nodes['x'] = nodes.geometry.x
nodes['y'] = nodes.geometry.y
nodes = nodes.set_index('nodeID')
nodes = nodes.rename_axis(index={'nodeID':'osmid'})
edges['key'] = 0
edges = edges.set_index(['node_start','node_end','key'])
edges = edges.rename_axis(index={'node_start':'u','node_end':'v'})
graph_ox = ox.graph_from_gdfs(nodes, edges)

生成的graph_ox是MultiDiGraph，是有向图，而通过momepy根据坐标信息自动生成的图实际是无向的，即从node_start到node_end的流向并不正确，需要改成无向图：

graph_oxu = ox.get_undirected(graph_ox)

3. 最短路径

def getPath(route, edges):
    gdf_path = pd.DataFrame()
    path_list = []
    for i in range(len(route)-1):
        node1 = route[i]
        node2 = route[i+1]
        e1 = edges.loc[(edges.index.get_level_values('u')==node1)&(edges.index.get_level_values('v')==node2)]
        if len(e1)==0:
            e1 = edges.loc[(edges.index.get_level_values('u')==node2)&(edges.index.get_level_values('v')==node1)]
        assert len(e1)==1, f'无边或边不唯一'
        path_list.append(e1)
    if len(path_list):
        gdf_path = pd.concat(path_list)
    return gdf_path

sPoint = Point(124.34460452690661, 40.20653233525042) # 起点坐标
ePoint = Point(124.31843430480932, 40.175612247996426) # 终点坐标
sNode = ox.distance.nearest_nodes(graph_oxu, sPoint.x, sPoint.y) # 起始结点
eNode = ox.distance.nearest_nodes(graph_oxu, ePoint.x, ePoint.y) # 终止结点
route = nx.shortest_path(G=graph_oxu, source=sNode, target=eNode, weight='length') # 最短路径沿线结点
gdf_path = getPath(route, edges) # 根据结点找到边

起点1和终点2的位置，以及找到的最短路径效果如图：

Python打开图片_最短路径_05

找到的路径与输入点的位置有差异，是因为最短路径以图的结点为起终点，如果希望精确到输入点位，则可以求输入点到线的投影点，在投影点处新建结点，重构路网图，效果如下：

Python打开图片_python_06

总结

展示的是只有线要素的情况下，如何生成方便操作的图，如果边和结点数据都有，则可以直接用osmnx.graph_from_gdfs(nodes, edges)生成有向图，不需要这么麻烦。
生成的最短路径是以图的结点为起终点，要精确到任意输入点的位置，需要重构路网图。