北京六环附近及往内的可驾驶道路网络（路网graph为连通图）

一、北京六环的道路.osm文件（beijing_car_seg_6ring road.osm，相关文件存网盘，网盘地址链接：https://pan.baidu.com/s/1ODT1BsN1HhyS1rTJW2v4og  提取码：szeq），如下图（JOSM软件查看）

 

二、解析北京六环的道路.osm文件（beijing_car_seg_6ring road.osm）生成6 ring Polygon geometry.

三、osmnx: Create a graph from OSM within the boundaries of some shapely polygon geometry.

四、完整python

from shapely.geometry import Polygon, LinearRing
import xml.etree.cElementTree as et
import matplotlib.pyplot as plt
from networkx import DiGraph
import geopandas as gpd
import numpy as np
# beijing_car_seg_6ring road.osm as the input data
root = et.parse(r"./beijing_car_seg_6ring road.osm").getroot()
nodes = {}
for node in root.findall('node'):
    nodes[int(node.attrib['id'])] = (float(node.attrib['lon']), float(node.attrib['lat']))
edges = []
for way in root.findall('way'):
    element_nd = way.findall('nd')
    node_s, node_e = int(element_nd[0].attrib['ref']), int(element_nd[1].attrib['ref'])
    path = (node_s, node_e)
    edges.append(path)
edge = edges[0]
E = []
E.append(edge)
edges.remove(edge)
point_s, point_e = nodes[E[0][0]], nodes[E[0][1]]
Point_lists = []
Point_lists.append(list(point_s))
Point_lists.append(list(point_e))
while len(edges) > 0:
    (node_f_s, node_f_e) = E[-1]
    for (node_n_s, node_n_e) in edges:
        if node_f_e == node_n_s:
            E.append((node_n_s, node_n_e))
            point_f_e = nodes[node_n_e]
            Point_lists.append(list(point_f_e))
            # print((node_n_s, node_n_e))
            edges.remove((node_n_s, node_n_e))
            break
# Points.pop()
print(E[0], '...', E[-2], E[-1])
print(Point_lists[0], '...', Point_lists[-2], Point_lists[-1])

road_line_arr = np.array(Point_lists)# 转换成二维坐标表示
six_ring_plg = Polygon(road_line_arr)# Polygon
six_ring_lr = LinearRing(road_line_arr)# LinearRing
six_ring_lr_buf = six_ring_lr.buffer(0.0225, 16)# LinearRing buffer
six_ring_lr_off = Polygon(six_ring_lr.parallel_offset(0.02)[-1])# LinearRing parallel offset

# Draw the geo with geopandas
geo_plg = gpd.GeoSeries(six_ring_plg)# Polygon
geo_plg.plot()
plt.show()

geo_lrg = gpd.GeoSeries(six_ring_lr)# LinearRing
geo_lrg.plot()
plt.show()

geo_lrg_buf = gpd.GeoSeries(six_ring_lr_buf)# LinearRing buffer
geo_lrg_buf.plot()
plt.show()

geo_lrg_off = gpd.GeoSeries(six_ring_lr_off)# LinearRing parallel offset
geo_lrg_off.plot()
plt.show()

from osmnx import geocoder, graph_from_polygon, graph_from_xml
from osmnx import save_graph_shapefile, save_graph_xml, plot_graph

######## Create a graph from OSM and save the graph to .osm ########
# 下载数据
# Create a graph from OSM within the boundaries of some shapely polygon.
G = graph_from_polygon(polygon=six_ring_lr_off,
                       network_type='drive',# network_type : {"all_private", "all", "bike", "drive", "drive_service", "walk"}
                       simplify=False,
                       truncate_by_edge=True)
# 可能因为代理问题， 则需要 Win+R -- (输入)regedit -- 删除 计算机HKEY_CURRENT_USERSOFTWAREMicrosoftWindowsCurrentVersionInternet Settings下的 ProxyEnable、ProxyOverride、ProxyServer 后重新运行程序
plot_graph(G, figsize=(32, 32), node_size=15, dpi=600, save=True, filepath='./images/image.png')
save_graph_xml(G, filepath='./osm/six_ring_graph.osm')
save_graph_shapefile(G, 'six_ring_shp')
# Read graph from xml(.osm)
M = graph_from_xml('./osm/six_ring_graph.osm', simplify=False)
print('number_of_nodes:', M.number_of_nodes(), 'number_of_edges:', M.number_of_edges())
plot_graph(M, figsize=(32, 32), node_size=15, dpi=600, save=True, filepath='./images/six_ring_graph.png')

# Convert M to DiGraph G
G = DiGraph(M, directed=True, n_res=5)

nodes = G.nodes()   # 获取图中的所有节点
nd_data = G.nodes.data(data=True)  # 获得所有节点的属性
nd_attr = nd_data[25248785]   # 获得某节点的属性 osmid=25248785
eg_data = G.edges(data=True)  # 获得所有边的属性
eg_attr = list(eg_data)[0]   # 获得某边-[0]的所有属性
eg_at_dic = eg_attr[2]   # 获得某边-[0]的第[2]组属性
eg_res_cost = eg_at_dic['res_cost']

# 根据地址查询经纬度信息
# Geocode the address string to a (lat, lng) point
address = 'Ceramstraat, Delft, Netherlands'# Adjust 根据需要修改 地址
address = 'Ceramstraat'# Adjust 根据需要修改 地址
point = geocoder.geocode(query=address)
print('(lat, lng) -->', point)