get_beijing_roadnetwork（工作太忙，仅仅作为记录）
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}\n{1},{2}\n{3},{4}\n".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_USER\SOFTWARE\Microsoft\Windows\CurrentVersion\Internet 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())
geo_convert.py
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))
个人学习记录
