1 from shapely.geometry import Polygon, LinearRing
2 import xml.etree.cElementTree as et
3 import matplotlib.pyplot as plt
4 from networkx import DiGraph
5 import geopandas as gpd
6 import numpy as np
7 from pyproj import CRS
8
9 # default CRS to set when creating graphs
10 default_crs = "epsg:4326"
11
12 def is_projected(crs):
13 """
14 Determine if a coordinate reference system is projected or not.
15
16 This is a convenience wrapper around the pyproj.CRS.is_projected function.
17
18 Parameters
19 ----------
20 crs : string or pyproj.CRS
21 the coordinate reference system
22
23 Returns
24 -------
25 projected : bool
26 True if crs is projected, otherwise False
27 """
28 return CRS.from_user_input(crs).is_projected
29
30 def project_gdf(gdf, to_crs=None, to_latlong=False):
31 """
32 Project a GeoDataFrame from its current CRS to another.
33
34 If to_crs is None, project to the UTM CRS for the UTM zone in which the
35 GeoDataFrame's centroid lies. Otherwise project to the CRS defined by
36 to_crs. The simple UTM zone calculation in this function works well for
37 most latitudes, but may not work for some extreme northern locations like
38 Svalbard or far northern Norway.
39
40 Parameters
41 ----------
42 gdf : geopandas.GeoDataFrame
43 the GeoDataFrame to be projected
44 to_crs : string or pyproj.CRS
45 if None, project to UTM zone in which gdf's centroid lies, otherwise
46 project to this CRS
47 to_latlong : bool
48 if True, project to settings.default_crs and ignore to_crs
49
50 Returns
51 -------
52 gdf_proj : geopandas.GeoDataFrame
53 the projected GeoDataFrame
54 """
55 if gdf.crs is None or len(gdf) < 1: # pragma: no cover
56 raise ValueError("GeoDataFrame must have a valid CRS and cannot be empty")
57
58 # if to_latlong is True, project the gdf to latlong
59 if to_latlong:
60 gdf_proj = gdf.to_crs(default_crs)
61 print(f"Projected GeoDataFrame to {default_crs}")
62
63 # else if to_crs was passed-in, project gdf to this CRS
64 elif to_crs is not None:
65 gdf_proj = gdf.to_crs(to_crs)
66 print(f"Projected GeoDataFrame to {to_crs}")
67
68 # otherwise, automatically project the gdf to UTM
69 else:
70 if is_projected(gdf.crs): # pragma: no cover
71 raise ValueError("Geometry must be unprojected to calculate UTM zone")
72
73 # calculate longitude of centroid of union of all geometries in gdf
74 avg_lng = gdf["geometry"].unary_union.centroid.x
75
76 # calculate UTM zone from avg longitude to define CRS to project to
77 utm_zone = int(np.floor((avg_lng + 180) / 6) + 1)
78 utm_crs = f"+proj=utm +zone={utm_zone} +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
79
80 # project the GeoDataFrame to the UTM CRS
81 gdf_proj = gdf.to_crs(utm_crs)
82 print(f"Projected GeoDataFrame to {gdf_proj.crs}")
83
84 return gdf_proj
85
86 def project_geometry(geometry, crs=None, to_crs=None, to_latlong=False):
87 """
88 Project a shapely geometry from its current CRS to another.
89
90 If to_crs is None, project to the UTM CRS for the UTM zone in which the
91 geometry's centroid lies. Otherwise project to the CRS defined by to_crs.
92
93 Parameters
94 ----------
95 geometry : shapely.geometry.Polygon or shapely.geometry.MultiPolygon
96 the geometry to project
97 crs : string or pyproj.CRS
98 the starting CRS of the passed-in geometry. if None, it will be set to
99 settings.default_crs
100 to_crs : string or pyproj.CRS
101 if None, project to UTM zone in which geometry's centroid lies,
102 otherwise project to this CRS
103 to_latlong : bool
104 if True, project to settings.default_crs and ignore to_crs
105
106 Returns
107 -------
108 geometry_proj, crs : tuple
109 the projected geometry and its new CRS
110 """
111 if crs is None:
112 crs = default_crs
113
114 gdf = gpd.GeoDataFrame(geometry=[geometry], crs=crs)
115 gdf_proj = project_gdf(gdf, to_crs=to_crs, to_latlong=to_latlong)
116 geometry_proj = gdf_proj["geometry"].iloc[0]
117 return geometry_proj, gdf_proj.crs
118
119 def get_bj_ring_road(osmfile=r"./beijing_6ring_multiways.osm", buffer_distance=1000, show=True):
120 # # beijing_xring road.osm as the input data
121 root = et.parse(osmfile).getroot()
122
123 nodes = {}
124 for node in root.findall('node'):
125 nodes[int(node.attrib['id'])] = (float(node.attrib['lon']), float(node.attrib['lat']))
126 edges = []
127 way_count = len(root.findall('way'))
128 if way_count > 1:# multiways in osm
129 for way in root.findall('way'):
130 element_nd = way.findall('nd')
131 node_s, node_e = int(element_nd[0].attrib['ref']), int(element_nd[1].attrib['ref'])
132 path = (node_s, node_e)
133 edges.append(path)
134 elif way_count == 1:# one way in osm
135 way_nodes = []
136 for element_nd in root.findall('way')[0].findall('nd'):
137 way_nd_index = element_nd.attrib['ref']
138 way_nodes.append(int(way_nd_index))
139 for i in range(len(way_nodes)-1):# a way must include two points
140 node_s, node_e = way_nodes[i], way_nodes[i+1]
141 path = (node_s, node_e)
142 edges.append(path)
143 else:# error osm
144 print("OSM Error!")
145 return None, None, None, None
146
147 edge = edges[0]
148 E = []
149 E.append(edge)
150 edges.remove(edge)
151 point_s, point_e = nodes[E[0][0]], nodes[E[0][1]]
152 Point_lists = []
153 Point_lists.append(list(point_s))
154 Point_lists.append(list(point_e))
155 while len(edges) > 0:
156 (node_f_s, node_f_e) = E[-1]
157 for (node_n_s, node_n_e) in edges:
158 if node_f_e == node_n_s:
159 E.append((node_n_s, node_n_e))
160 point_f_e = nodes[node_n_e]
161 Point_lists.append(list(point_f_e))
162 # print((node_n_s, node_n_e))
163 edges.remove((node_n_s, node_n_e))
164 break
165 # Points.pop()
166 print(E[0], '...', E[-2], E[-1])
167 print(Point_lists[0], '...', Point_lists[-2], Point_lists[-1])
168 road_line_arr = np.array(Point_lists) # 转换成二维坐标表示
169
170 ring_plg = Polygon(road_line_arr) # Polygon
171 ring_lr = LinearRing(road_line_arr) # LinearRing
172
173 ring_lr_proj, crs_utm = project_geometry(ring_lr)# 平面投影
174 ring_lr_proj_buff =ring_lr_proj.buffer(buffer_distance)# buffer
175 ring_lr_buf, ring_lr_buf_epsg = project_geometry(ring_lr_proj_buff, crs=crs_utm, to_latlong=True) # 反平面投影 LinearRing buffer
176
177 poly_proj, crs_utm = project_geometry(ring_plg)# 平面投影
178 poly_proj_buff = poly_proj.buffer(buffer_distance)# buffer
179 ring_plg_buf, poly_buff_epsg = project_geometry(poly_proj_buff, crs=crs_utm, to_latlong=True)# 反平面投影 Polygon buffer
180
181 if show:
182 # Draw the geo with geopandas
183 geo_lrg = gpd.GeoSeries(ring_lr) # LinearRing
184 geo_lrg.plot()
185 plt.xlabel("$lon$")
186 plt.ylabel("$lat$")
187 plt.title("$LinearRing$")
188 plt.show()
189
190 geo_lrg_buf = gpd.GeoSeries(ring_lr_buf) # LinearRing buffer
191 geo_lrg_buf.plot()
192 plt.xlabel("$lon$")
193 plt.ylabel("$lat$")
194 plt.title("$LinearRing-with-buffer$")
195 plt.show()
196
197 geo_plg = gpd.GeoSeries(ring_plg) # Polygon
198 geo_plg.plot()
199 plt.xlabel("$lon$")
200 plt.ylabel("$lat$")
201 plt.title("$Polygon$")
202 plt.show()
203
204 geo_plg = gpd.GeoSeries(ring_plg_buf) # Polygon
205 geo_plg.plot()
206 plt.xlabel("$lon$")
207 plt.ylabel("$lat$")
208 plt.title("$Polygon-with-buffer$")
209 plt.show()
210
211 return ring_lr, ring_lr_buf, ring_plg, ring_plg_buf
212
213 # road_6ring_lr, road_6ring_lr_buf, road_6ring_plg, road_6ring_plg_buf = get_bj_ring_road(osmfile=r"./beijing_6ring_multiways.osm", buffer_distance=1000, show=True)
214 road_2ring_lr, road_2ring_lr_buf, road_2ring_plg, road_2ring_plg_buf = get_bj_ring_road(osmfile=r"./beijing_2ring.osm", buffer_distance=1000, show=True)
215
216 # test
217 # # buffer_dist = 500
218 # # poly_proj, crs_utm = project_geometry(road_2ring_plg)
219 # # poly_proj_buff = poly_proj.buffer(buffer_dist)
220 # # poly_buff, poly_buff_epsg = project_geometry(poly_proj_buff, crs=crs_utm, to_latlong=True)
221 # # ring_lr_proj, crs_utm = project_geometry(road_2ring_lr)
222 # # ring_lr_proj_buff = ring_lr_proj.buffer(buffer_dist)
223 # # ring_lr_buff, ring_lr_buff_epsg = project_geometry(ring_lr_proj_buff, crs=crs_utm, to_latlong=True)
224
225 import osmnx as ox
226 from osmnx import geocoder, graph_from_polygon, graph_from_xml
227 from osmnx import save_graph_shapefile, save_graph_xml, plot_graph
228
229 ######## Create a graph from OSM and save the graph to .osm ########
230 # 下载数据
231
232 # ox.settings
233 utn = ox.settings.useful_tags_node
234 oxna = ox.settings.osm_xml_node_attrs
235 oxnt = ox.settings.osm_xml_node_tags
236 utw = ox.settings.useful_tags_way
237 oxwa = ox.settings.osm_xml_way_attrs
238 oxwt = ox.settings.osm_xml_way_tags
239 utn = list(set(utn + oxna + oxnt))
240 utw = list(set(utw + oxwa + oxwt))
241 ox.config(all_oneway=True, useful_tags_node=utn, useful_tags_way=utw)
242
243 # Create a graph from OSM within the boundaries of some shapely polygon.
244 M = graph_from_polygon(polygon=road_2ring_plg_buf,
245 network_type='drive', # network_type : {"all_private", "all", "bike", "drive", "drive_service", "walk"}
246 simplify=True,# if True, simplify graph topology with the `simplify_graph` function
247 retain_all=False,# if True, return the entire graph even if it is not connected.
248 truncate_by_edge=True,# if True, retain nodes outside boundary polygon if at least one of node's neighbors is within the polygon
249 clean_periphery=True,# if True, buffer 500m to get a graph larger than requested, then simplify, then truncate it to requested spatial boundaries
250 custom_filter=None)# a custom ways filter to be used instead of the network_type presets, e.g., '["power"~"line"]' or '["highway"~"motorway|trunk"]'.
251 # 可能因为代理问题， 则需要 Win+R -- (输入)regedit -- 删除 计算机\HKEY_CURRENT_USER\SOFTWARE\Microsoft\Windows\CurrentVersion\Internet Settings下的 ProxyEnable、ProxyOverride、ProxyServer 后重新运行程序
252 print('number_of_nodes:', M.number_of_nodes(), 'number_of_edges:', M.number_of_edges())
253 plot_graph(M, figsize=(32, 32), node_size=15, dpi=600, save=True, filepath='./images/ring_graph.png')
254 # save_graph_xml(M, filepath='./download_osm/graph.osm')
255 save_graph_shapefile(M, 'download_shp')
256
257 # Convert M to DiGraph G
258 G = DiGraph(M, directed=True)
259
260 nodes = G.nodes() # 获取图中的所有节点
261 nd_data = G.nodes.data(data=True) # 获得所有节点的属性
262 print(list(nd_data)[-1])
263 # simplify (9094448576, {'y': 39.8708445, 'x': 116.3840209, 'street_count': 3})
264 # full (9094448578, {'y': 39.8707558, 'x': 116.3840692, 'highway': 'crossing', 'street_count': 2})
265 eg_data = G.edges(data=True) # 获得所有边的属性
266 print(list(eg_data)[-1])
267 # simplify (9094448576, 9063939428, {'osmid': [955313155, 26295235, 955313156], 'oneway': 'yes', 'highway': 'tertiary', 'length': 393.366, 'tunnel': 'yes', 'geometry': <shapely.geometry.linestring.LineString object at 0x0000020FBD457630>})
268 # full (9094448578, 9094448574, {'osmid': 983612462, 'oneway': 'yes', 'highway': 'trunk_link', 'length': 3.189})
269 print('number_of_nodes:', G.number_of_nodes())
270 print('number_of_edges:', G.number_of_edges())
北京23456环路OSM数据见百度网盘：链接：https://pan.baidu.com/s/1Y3H0s9natD77rgDP-kMvxQ
提取码：8866