in libs/apls/osmnx_funcs.py [0:0]
def project_graph(G, to_crs=None, verbose=False):
"""
https://github.com/gboeing/osmnx/blob/v0.9/osmnx/projection.py#L126
Project a graph from lat-long to the UTM zone appropriate for its geographic
location.
Parameters
----------
G : networkx multidigraph
the networkx graph to be projected
to_crs : dict
if not None, just project to this CRS instead of to UTM
Returns
-------
networkx multidigraph
"""
G_proj = G.copy()
start_time = time.time()
# create a GeoDataFrame of the nodes, name it, convert osmid to str
nodes, data = zip(*G_proj.nodes(data=True))
gdf_nodes = gpd.GeoDataFrame(list(data), index=nodes)
# gdf_nodes.crs = G_proj.graph['crs']
gdf_nodes.gdf_name = '{}_nodes'.format(G_proj.name)
# create new lat/lon columns just to save that data for later, and create a
# geometry column from x/y
gdf_nodes['lon'] = gdf_nodes['x']
gdf_nodes['lat'] = gdf_nodes['y']
gdf_nodes['geometry'] = gdf_nodes.apply(lambda row: Point(row['x'], row['y']), axis=1)
if verbose:
print('Created a GeoDataFrame from graph in {:,.2f} seconds'.format(time.time()-start_time))
gdf_nodes.crs = default_crs
# project the nodes GeoDataFrame to UTM
gdf_nodes_utm = project_gdf(gdf_nodes, to_crs=to_crs)
# extract data for all edges that have geometry attribute
edges_with_geom = []
for u, v, key, data in G_proj.edges(keys=True, data=True):
if 'geometry' in data:
edges_with_geom.append({'u':u, 'v':v, 'key':key, 'geometry':data['geometry']})
# create an edges GeoDataFrame and project to UTM, if there were any edges
# with a geometry attribute. geom attr only exists if graph has been
# simplified, otherwise you don't have to project anything for the edges
# because the nodes still contain all spatial data
if len(edges_with_geom) > 0:
gdf_edges = gpd.GeoDataFrame(edges_with_geom)
gdf_edges.crs = G_proj.graph['crs']
gdf_edges.gdf_name = '{}_edges'.format(G_proj.name)
gdf_edges_utm = project_gdf(gdf_edges, to_crs=to_crs)
# extract projected x and y values from the nodes' geometry column
start_time = time.time()
gdf_nodes_utm['x'] = gdf_nodes_utm['geometry'].map(lambda point: point.x)
gdf_nodes_utm['y'] = gdf_nodes_utm['geometry'].map(lambda point: point.y)
gdf_nodes_utm = gdf_nodes_utm.drop('geometry', axis=1)
if verbose:
print('Extracted projected node geometries from GeoDataFrame in {:,.2f} seconds'.format(time.time()-start_time))
# clear the graph to make it a blank slate for the projected data
start_time = time.time()
edges = list(G_proj.edges(keys=True, data=True))
graph_name = G_proj.graph['name']
G_proj.clear()
# add the projected nodes and all their attributes to the graph
G_proj.add_nodes_from(gdf_nodes_utm.index)
attributes = gdf_nodes_utm.to_dict()
for label in gdf_nodes_utm.columns:
nx.set_node_attributes(G_proj, name=label, values=attributes[label])
# add the edges and all their attributes (including reconstructed geometry,
# when it exists) to the graph
for u, v, key, attributes in edges:
if 'geometry' in attributes:
row = gdf_edges_utm[(gdf_edges_utm['u']==u) & (gdf_edges_utm['v']==v) & (gdf_edges_utm['key']==key)]
attributes['geometry'] = row['geometry'].iloc[0]
# attributes dict contains key, so we don't need to explicitly pass it here
G_proj.add_edge(u, v, **attributes)
# set the graph's CRS attribute to the new, projected CRS and return the
# projected graph
G_proj.graph['crs'] = gdf_nodes_utm.crs
G_proj.graph['name'] = '{}_UTM'.format(graph_name)
if 'streets_per_node' in G.graph:
G_proj.graph['streets_per_node'] = G.graph['streets_per_node']
if verbose:
print('Rebuilt projected graph in {:,.2f} seconds'.format(time.time()-start_time))
return G_proj