北京六环附近及往内的可驾驶道路网络(路网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_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)
个人学习记录