北京六环附近及往内的可驾驶道路网络（路网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

?

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98

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_USER\SOFTWARE\Microsoft\Windows\CurrentVersion\Internet 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)