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 import requests
9 from tools.geo_convert import wgs84_to_gcj02, gcj02_to_wgs84
10
11 # default CRS to set when creating graphs
12 default_crs = "epsg:4326"
13
14 def is_projected(crs):
15 """
16 Determine if a coordinate reference system is projected or not.
17
18 This is a convenience wrapper around the pyproj.CRS.is_projected function.
19
20 Parameters
21 ----------
22 crs : string or pyproj.CRS
23 the coordinate reference system
24
25 Returns
26 -------
27 projected : bool
28 True if crs is projected, otherwise False
29 """
30 return CRS.from_user_input(crs).is_projected
31
32 def project_gdf(gdf, to_crs=None, to_latlong=False):
33 """
34 Project a GeoDataFrame from its current CRS to another.
35
36 If to_crs is None, project to the UTM CRS for the UTM zone in which the
37 GeoDataFrame's centroid lies. Otherwise project to the CRS defined by
38 to_crs. The simple UTM zone calculation in this function works well for
39 most latitudes, but may not work for some extreme northern locations like
40 Svalbard or far northern Norway.
41
42 Parameters
43 ----------
44 gdf : geopandas.GeoDataFrame
45 the GeoDataFrame to be projected
46 to_crs : string or pyproj.CRS
47 if None, project to UTM zone in which gdf's centroid lies, otherwise
48 project to this CRS
49 to_latlong : bool
50 if True, project to settings.default_crs and ignore to_crs
51
52 Returns
53 -------
54 gdf_proj : geopandas.GeoDataFrame
55 the projected GeoDataFrame
56 """
57 if gdf.crs is None or len(gdf) < 1: # pragma: no cover
58 raise ValueError("GeoDataFrame must have a valid CRS and cannot be empty")
59
60 # if to_latlong is True, project the gdf to latlong
61 if to_latlong:
62 gdf_proj = gdf.to_crs(default_crs)
63 print(f"Projected GeoDataFrame to {default_crs}")
64
65 # else if to_crs was passed-in, project gdf to this CRS
66 elif to_crs is not None:
67 gdf_proj = gdf.to_crs(to_crs)
68 print(f"Projected GeoDataFrame to {to_crs}")
69
70 # otherwise, automatically project the gdf to UTM
71 else:
72 if is_projected(gdf.crs): # pragma: no cover
73 raise ValueError("Geometry must be unprojected to calculate UTM zone")
74
75 # calculate longitude of centroid of union of all geometries in gdf
76 avg_lng = gdf["geometry"].unary_union.centroid.x
77
78 # calculate UTM zone from avg longitude to define CRS to project to
79 utm_zone = int(np.floor((avg_lng + 180) / 6) + 1)
80 utm_crs = f"+proj=utm +zone={utm_zone} +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
81
82 # project the GeoDataFrame to the UTM CRS
83 gdf_proj = gdf.to_crs(utm_crs)
84 print(f"Projected GeoDataFrame to {gdf_proj.crs}")
85
86 return gdf_proj
87
88 def project_geometry(geometry, crs=None, to_crs=None, to_latlong=False):
89 """
90 Project a shapely geometry from its current CRS to another.
91
92 If to_crs is None, project to the UTM CRS for the UTM zone in which the
93 geometry's centroid lies. Otherwise project to the CRS defined by to_crs.
94
95 Parameters
96 ----------
97 geometry : shapely.geometry.Polygon or shapely.geometry.MultiPolygon
98 the geometry to project
99 crs : string or pyproj.CRS
100 the starting CRS of the passed-in geometry. if None, it will be set to
101 settings.default_crs
102 to_crs : string or pyproj.CRS
103 if None, project to UTM zone in which geometry's centroid lies,
104 otherwise project to this CRS
105 to_latlong : bool
106 if True, project to settings.default_crs and ignore to_crs
107
108 Returns
109 -------
110 geometry_proj, crs : tuple
111 the projected geometry and its new CRS
112 """
113 if crs is None:
114 crs = default_crs
115
116 gdf = gpd.GeoDataFrame(geometry=[geometry], crs=crs)
117 gdf_proj = project_gdf(gdf, to_crs=to_crs, to_latlong=to_latlong)
118 geometry_proj = gdf_proj["geometry"].iloc[0]
119 return geometry_proj, gdf_proj.crs
120
121 def get_grid_from_polygon(polygon_lonlat, threshold = 1000.0, gcj=True):
122 bounds = []
123 ring_plg_proj, crs_utm = project_geometry(polygon_lonlat) # 平面投影
124 ring_plg_proj_buff = ring_plg_proj.buffer(threshold / 2.0) # buffer dis = threshold/2.0
125 ring_plg_buf, ring_plg_buf_epsg = project_geometry(ring_plg_proj_buff, crs=crs_utm, to_latlong=True) # 反平面投影
126 ring_plg_enve = ring_plg_buf.envelope # Polygon
127 ring_plg_enve_proj, crs_utm = project_geometry(ring_plg_enve) # 平面投影
128 (X_l, Y_b, X_r, Y_u) = ring_plg_enve_proj.bounds # tuple
129 delt_X, delt_Y = X_r - X_l, Y_u - Y_b
130 m, n = int(np.ceil(delt_X / threshold)), int(np.ceil(delt_Y / threshold))
131 for i in range(n):
132 y_lb = threshold * i + Y_b
133 y_ur = y_lb + threshold
134 for j in range(m):
135 x_lb = threshold * j + X_l
136 x_ur = x_lb + threshold
137 poly_arr = np.array([[x_lb, y_lb], [x_ur, y_lb], [x_ur, y_ur], [x_lb, y_ur]])
138 grid_plg_tmp = Polygon(poly_arr)
139 grid_plg, grid_plg_epsg = project_geometry(grid_plg_tmp, crs=crs_utm, to_latlong=True) # 反平面投影
140 geom_intersec = polygon_lonlat.intersection(grid_plg) # Polygon
141 if geom_intersec.is_empty is not True:# 存在空间相交
142 (lon_l, lat_b, lon_r, lat_u) = grid_plg.bounds # tuple
143
144 if gcj:
145 lon_l, lat_b = wgs84_to_gcj02(lng=lon_l, lat=lat_b)
146 lon_r, lat_u = wgs84_to_gcj02(lng=lon_r, lat=lat_u)
147
148 bounds.append((lon_l, lat_b, lon_r, lat_u))
149
150 return bounds
151
152 def get_region_txt_file(bounds, region_txt_file = './region.txt'):
153 num = 1
154 file_region = open(region_txt_file, 'w')
155 for (lon_l, lat_b, lon_r, lat_u) in bounds:
156 zs_lon, zs_lat = lon_l, lat_u
157 yx_lon, yx_lat = lon_r, lat_b
158 region_str = "{0}
{1},{2}
{3},{4}
".format(num, zs_lon, zs_lat, yx_lon, yx_lat)
159 print(region_str)
160 file_region.write(region_str)
161 num += 1
162 file_region.close()
163
164 def point_transform_amap(location):
165 """
166 坐标转换，将非高德坐标（GPS坐标、mapbar坐标、baidu坐标）转换为高德坐标
167 """
168 parameters = {
169 'coordsys': 'gps',
170 'locations': location,
171 'key': '***'# gaod key
172 }
173 base = 'http://restapi.amap.com/v3/assistant/coordinate/convert?'
174 response = requests.get(base, parameters)
175 answer = response.json()
176 return answer['locations']
177
178 def get_bj_ring_road(osmfile=r"./beijing_6ring_multiways.osm", buffer_distance=1000, show=True):
179 # # beijing_xring road.osm as the input data
180 root = et.parse(osmfile).getroot()
181
182 nodes = {}
183 for node in root.findall('node'):
184 nodes[int(node.attrib['id'])] = (float(node.attrib['lon']), float(node.attrib['lat']))
185 edges = []
186 way_count = len(root.findall('way'))
187 if way_count > 1:# multiways in osm
188 for way in root.findall('way'):
189 element_nd = way.findall('nd')
190 node_s, node_e = int(element_nd[0].attrib['ref']), int(element_nd[1].attrib['ref'])
191 path = (node_s, node_e)
192 edges.append(path)
193 elif way_count == 1:# one way in osm
194 way_nodes = []
195 for element_nd in root.findall('way')[0].findall('nd'):
196 way_nd_index = element_nd.attrib['ref']
197 way_nodes.append(int(way_nd_index))
198 for i in range(len(way_nodes)-1):# a way must include two points
199 node_s, node_e = way_nodes[i], way_nodes[i+1]
200 path = (node_s, node_e)
201 edges.append(path)
202 else:# error osm
203 print("OSM Error!")
204 return None, None, None, None
205
206 edge = edges[0]
207 E = []
208 E.append(edge)
209 edges.remove(edge)
210 point_s, point_e = nodes[E[0][0]], nodes[E[0][1]]
211 Point_lists = []
212 Point_lists.append(list(point_s))
213 Point_lists.append(list(point_e))
214 while len(edges) > 0:
215 (node_f_s, node_f_e) = E[-1]
216 for (node_n_s, node_n_e) in edges:
217 if node_f_e == node_n_s:
218 E.append((node_n_s, node_n_e))
219 point_f_e = nodes[node_n_e]
220 Point_lists.append(list(point_f_e))
221 # print((node_n_s, node_n_e))
222 edges.remove((node_n_s, node_n_e))
223 break
224 # Points.pop()
225 print(E[0], '...', E[-2], E[-1])
226 print(Point_lists[0], '...', Point_lists[-2], Point_lists[-1])
227 road_line_arr = np.array(Point_lists) # 转换成二维坐标表示
228
229 ring_plg = Polygon(road_line_arr) # Polygon
230 ring_lr = LinearRing(road_line_arr) # LinearRing
231
232 # get_grid_from_polygon and get_region_txt_file
233 bounds = get_grid_from_polygon(ring_plg, threshold=1000.0, gcj=True) # get_grid_from_polygon with threshold
234 get_region_txt_file(bounds, region_txt_file='./region.txt')
235
236 ring_lr_proj, crs_utm = project_geometry(ring_lr)# 平面投影
237 ring_lr_proj_buff =ring_lr_proj.buffer(buffer_distance)# buffer
238 ring_lr_buf, ring_lr_buf_epsg = project_geometry(ring_lr_proj_buff, crs=crs_utm, to_latlong=True) # 反平面投影 LinearRing buffer
239
240 poly_proj, crs_utm = project_geometry(ring_plg)# 平面投影
241 poly_proj_buff = poly_proj.buffer(buffer_distance)# buffer
242 ring_plg_buf, poly_buff_epsg = project_geometry(poly_proj_buff, crs=crs_utm, to_latlong=True)# 反平面投影 Polygon buffer
243
244 if show:
245 # Draw the geo with geopandas
246 geo_lrg = gpd.GeoSeries(ring_lr) # LinearRing
247 geo_lrg.plot()
248 plt.xlabel("$lon$")
249 plt.ylabel("$lat$")
250 plt.title("$LinearRing$")
251 plt.show()
252
253 geo_lrg_buf = gpd.GeoSeries(ring_lr_buf) # LinearRing buffer
254 geo_lrg_buf.plot()
255 plt.xlabel("$lon$")
256 plt.ylabel("$lat$")
257 plt.title("$LinearRing-with-buffer$")
258 plt.show()
259
260 geo_plg = gpd.GeoSeries(ring_plg) # Polygon
261 geo_plg.plot()
262 plt.xlabel("$lon$")
263 plt.ylabel("$lat$")
264 plt.title("$Polygon$")
265 plt.show()
266
267 geo_plg = gpd.GeoSeries(ring_plg_buf) # Polygon
268 geo_plg.plot()
269 plt.xlabel("$lon$")
270 plt.ylabel("$lat$")
271 plt.title("$Polygon-with-buffer$")
272 plt.show()
273
274 return ring_lr, ring_lr_buf, ring_plg, ring_plg_buf
275
276 # 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)
277 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)
278
279 # test
280 # # buffer_dist = 500
281 # # poly_proj, crs_utm = project_geometry(road_2ring_plg)
282 # # poly_proj_buff = poly_proj.buffer(buffer_dist)
283 # # poly_buff, poly_buff_epsg = project_geometry(poly_proj_buff, crs=crs_utm, to_latlong=True)
284 # # ring_lr_proj, crs_utm = project_geometry(road_2ring_lr)
285 # # ring_lr_proj_buff = ring_lr_proj.buffer(buffer_dist)
286 # # ring_lr_buff, ring_lr_buff_epsg = project_geometry(ring_lr_proj_buff, crs=crs_utm, to_latlong=True)
287
288 import osmnx as ox
289 from osmnx import geocoder, graph_from_polygon, graph_from_xml
290 from osmnx import save_graph_shapefile, save_graph_xml, plot_graph
291
292 ######## Create a graph from OSM and save the graph to .osm ########
293 # 下载数据
294
295 # ox.settings
296 utn = ox.settings.useful_tags_node
297 oxna = ox.settings.osm_xml_node_attrs
298 oxnt = ox.settings.osm_xml_node_tags
299 utw = ox.settings.useful_tags_way
300 oxwa = ox.settings.osm_xml_way_attrs
301 oxwt = ox.settings.osm_xml_way_tags
302 utn = list(set(utn + oxna + oxnt))
303 utw = list(set(utw + oxwa + oxwt))
304 ox.config(all_oneway=True, useful_tags_node=utn, useful_tags_way=utw)
305
306 # Create a graph from OSM within the boundaries of some shapely polygon.
307 M = graph_from_polygon(polygon=road_2ring_plg_buf,
308 network_type='drive', # network_type : {"all_private", "all", "bike", "drive", "drive_service", "walk"}
309 simplify=True,# if True, simplify graph topology with the `simplify_graph` function
310 retain_all=False,# if True, return the entire graph even if it is not connected.
311 truncate_by_edge=True,# if True, retain nodes outside boundary polygon if at least one of node's neighbors is within the polygon
312 clean_periphery=True,# if True, buffer 500m to get a graph larger than requested, then simplify, then truncate it to requested spatial boundaries
313 custom_filter=None)# a custom ways filter to be used instead of the network_type presets, e.g., '["power"~"line"]' or '["highway"~"motorway|trunk"]'.
314 # 可能因为代理问题， 则需要 Win+R -- (输入)regedit -- 删除 计算机HKEY_CURRENT_USERSOFTWAREMicrosoftWindowsCurrentVersionInternet Settings下的 ProxyEnable、ProxyOverride、ProxyServer 后重新运行程序
315 print('number_of_nodes:', M.number_of_nodes(), 'number_of_edges:', M.number_of_edges())
316 plot_graph(M, figsize=(32, 32), node_size=15, dpi=600, save=True, filepath='./images/ring_graph.png')
317 # save_graph_xml(M, filepath='./download_osm/graph.osm')
318 save_graph_shapefile(M, 'download_shp')
319
320 # Convert M to DiGraph G
321 G = DiGraph(M, directed=True)
322
323 nodes = G.nodes() # 获取图中的所有节点
324 nd_data = G.nodes.data(data=True) # 获得所有节点的属性
325 print(list(nd_data)[-1])
326 # simplify (9094448576, {'y': 39.8708445, 'x': 116.3840209, 'street_count': 3})
327 # full (9094448578, {'y': 39.8707558, 'x': 116.3840692, 'highway': 'crossing', 'street_count': 2})
328 eg_data = G.edges(data=True) # 获得所有边的属性
329 print(list(eg_data)[-1])
330 # simplify (9094448576, 9063939428, {'osmid': [955313155, 26295235, 955313156], 'oneway': 'yes', 'highway': 'tertiary', 'length': 393.366, 'tunnel': 'yes', 'geometry': <shapely.geometry.linestring.LineString object at 0x0000020FBD457630>})
331 # full (9094448578, 9094448574, {'osmid': 983612462, 'oneway': 'yes', 'highway': 'trunk_link', 'length': 3.189})
332 print('number_of_nodes:', G.number_of_nodes())
333 print('number_of_edges:', G.number_of_edges())
1 # coding: utf-8
2
3 """
4 Geolocataion conversion between WGS84, BD09 and GCJ02.
5 WGS84 / BD09 / GCJ02 / MapBar 经纬度坐标互转。
6
7 - WGS84: GPS coordinates for Google Earth (GPS坐标，谷歌地球使用)
8 - GCJ02: national coordinate system developed by China (国测局坐标，谷歌中国地图、腾讯地图、高德地图使用)
9 - BD09: Baidu coordinates (百度坐标系，百度地图使用)
10 - MapBar: MapBar coordinates (图吧坐标系，图吧地图使用)
11
12 Test website: http://gpsspg.com/maps.htm
13
14 Author: Gaussic
15 Date: 2019-05-09
16 """
17
18 import math
19
20 PI = math.pi
21 PIX = math.pi * 3000.0 / 180.0
22 EE = 0.00669342162296594323
23 A = 6378245.0
24
25
26 def bd09_to_gcj02(lng, lat):
27 """BD09 -> GCJ02"""
28 x, y = lng - 0.0065, lat - 0.006
29 z = math.sqrt(x * x + y * y) - 0.00002 * math.sin(y * PIX)
30 theta = math.atan2(y, x) - 0.000003 * math.cos(x * PIX)
31 lng, lat = z * math.cos(theta), z * math.sin(theta)
32 return lng, lat
33
34
35 def gcj02_to_bd09(lng, lat):
36 """GCJ02 -> BD09"""
37 z = math.sqrt(lng * lng + lat * lat) + 0.00002 * math.sin(lat * PIX)
38 theta = math.atan2(lat, lng) + 0.000003 * math.cos(lng * PIX)
39 lng, lat = z * math.cos(theta) + 0.0065, z * math.sin(theta) + 0.006
40 return lng, lat
41
42
43 def gcj02_to_wgs84(lng, lat):
44 """GCJ02 -> WGS84"""
45 if out_of_china(lng, lat):
46 return lng, lat
47 dlat = transform_lat(lng - 105.0, lat - 35.0)
48 dlng = transform_lng(lng - 105.0, lat - 35.0)
49 radlat = lat / 180.0 * PI
50 magic = math.sin(radlat)
51 magic = 1 - EE * magic * magic
52 sqrtmagic = math.sqrt(magic)
53 dlat = (dlat * 180.0) / ((A * (1 - EE)) / (magic * sqrtmagic) * PI)
54 dlng = (dlng * 180.0) / (A / sqrtmagic * math.cos(radlat) * PI)
55 lng, lat = lng - dlng, lat - dlat
56 return lng, lat
57
58
59 def wgs84_to_gcj02(lng, lat):
60 """WGS84 -> GCJ02"""
61 if out_of_china(lng, lat):
62 return lng, lat
63 dlat = transform_lat(lng - 105.0, lat - 35.0)
64 dlng = transform_lng(lng - 105.0, lat - 35.0)
65 radlat = lat / 180.0 * PI
66 magic = math.sin(radlat)
67 magic = 1 - EE * magic * magic
68 sqrtmagic = math.sqrt(magic)
69 dlat = (dlat * 180.0) / ((A * (1 - EE)) / (magic * sqrtmagic) * PI)
70 dlng = (dlng * 180.0) / (A / sqrtmagic * math.cos(radlat) * PI)
71 lng, lat = lng + dlng, lat + dlat
72 return lng, lat
73
74
75 def mapbar_to_wgs84(lng, lat):
76 """MapBar -> WGS84"""
77 lng = lng * 100000.0 % 36000000
78 lat = lat * 100000.0 % 36000000
79 lng1 = int(lng - math.cos(lat / 100000.0) * lng / 18000.0 - math.sin(lng / 100000.0) * lat / 9000.0)
80 lat1 = int(lat - math.sin(lat / 100000.0) * lng / 18000.0 - math.cos(lng / 100000.0) * lat / 9000.0)
81 lng2 = int(lng - math.cos(lat1 / 100000.0) * lng1 / 18000.0 - math.sin(lng1 / 100000.0) * lat1 / 9000.0 + (1 if lng > 0 else -1))
82 lat2 = int(lat - math.sin(lat1 / 100000.0) * lng1 / 18000.0 - math.cos(lng1 / 100000.0) * lat1 / 9000.0 + (1 if lat > 0 else -1))
83 lng, lat = lng2 / 100000.0, lat2 / 100000.0
84 return lng, lat
85
86
87 def transform_lat(lng, lat):
88 """GCJ02 latitude transformation"""
89 ret = -100 + 2.0 * lng + 3.0 * lat + 0.2 * lat * lat + 0.1 * lng * lat + 0.2 * math.sqrt(math.fabs(lng))
90 ret += (20.0 * math.sin(6.0 * lng * PI) + 20.0 * math.sin(2.0 * lng * PI)) * 2.0 / 3.0
91 ret += (20.0 * math.sin(lat * PI) + 40.0 * math.sin(lat / 3.0 * PI)) * 2.0 / 3.0
92 ret += (160.0 * math.sin(lat / 12.0 * PI) + 320.0 * math.sin(lat * PI / 30.0)) * 2.0 / 3.0
93 return ret
94
95
96 def transform_lng(lng, lat):
97 """GCJ02 longtitude transformation"""
98 ret = 300.0 + lng + 2.0 * lat + 0.1 * lng * lng + 0.1 * lng * lat + 0.1 * math.sqrt(math.fabs(lng))
99 ret += (20.0 * math.sin(6.0 * lng * PI) + 20.0 * math.sin(2.0 * lng * PI)) * 2.0 / 3.0
100 ret += (20.0 * math.sin(lng * PI) + 40.0 * math.sin(lng / 3.0 * PI)) * 2.0 / 3.0
101 ret += (150.0 * math.sin(lng / 12.0 * PI) + 300.0 * math.sin(lng / 30.0 * PI)) * 2.0 / 3.0
102 return ret
103
104
105 def out_of_china(lng, lat):
106 """No offset when coordinate out of China."""
107 if lng < 72.004 or lng > 137.8437:
108 return True
109 if lat < 0.8293 or lat > 55.8271:
110 return True
111 return False
112
113
114 def bd09_to_wgs84(lng, lat):
115 """BD09 -> WGS84"""
116 lng, lat = bd09_to_gcj02(lng, lat)
117 lng, lat = gcj02_to_wgs84(lng, lat)
118 return lng, lat
119
120
121 def wgs84_to_bd09(lng, lat):
122 """WGS84 -> BD09"""
123 lng, lat = wgs84_to_gcj02(lng, lat)
124 lng, lat = gcj02_to_bd09(lng, lat)
125 return lng, lat
126
127
128 def mapbar_to_gcj02(lng, lat):
129 """MapBar -> GCJ02"""
130 lng, lat = mapbar_to_wgs84(lng, lat)
131 lng, lat = wgs84_to_gcj02(lng, lat)
132 return lng, lat
133
134
135 def mapbar_to_bd09(lng, lat):
136 """MapBar -> BD09"""
137 lng, lat = mapbar_to_wgs84(lng, lat)
138 lng, lat = wgs84_to_bd09(lng, lat)
139 return lng, lat
140
141
142 if __name__ == '__main__':
143 blng, blat = 121.4681891220,31.1526609317
144 print('BD09:', (blng, blat))
145 print('BD09 -> GCJ02:', bd09_to_gcj02(blng, blat))
146 print('BD09 -> WGS84:',bd09_to_wgs84(blng, blat))
147 wlng, wlat = 121.45718237717077, 31.14846209914084
148 print('WGS84:', (wlng, wlat))
149 print('WGS84 -> GCJ02:', wgs84_to_gcj02(wlng, wlat))
150 print('WGS84 -> BD09:', wgs84_to_bd09(wlng, wlat))
151 mblng, mblat = 121.4667323772, 31.1450420991
152 print('MapBar:', (mblng, mblat))
153 print('MapBar -> WGS84:', mapbar_to_wgs84(mblng, mblat))
154 print('MapBar -> GCJ02:', mapbar_to_gcj02(mblng, mblat))
155 print('MapBar -> BD09:', mapbar_to_bd09(mblng, mblat))