#Import packages
import os
import osmnx as ox
import networkx as nx
import pandas as pd
import geopandas as gpd
from shapely.geometry import Point
import contextily as ctx
import matplotlib.pyplot as plt
from pyproj import CRS
#Filenames
DCFC_filename = 'Data\\NREL\\NC_DCFC.shp'
graph_filename = 'Data\\OSM\\NC_highways_all.graphml'
#Read DCFC points into geodataframe
print(f"Loading DCFC points from {DCFC_filename}")
gdf_DCFC = gpd.read_file(DCFC_filename)
#Read graphML file into a graph object
print(f"Loading graph from {graph_filename}")
nc_graph = ox.load_graphml(filename=os.path.basename(graph_filename),
folder=os.path.dirname(graph_filename))
#Convert graph component into geodataframes
gdf_nodes, gdf_edges = ox.graph_to_gdfs(nc_graph)
#Project to UTM (for analysis)
utm17N_prj = CRS.from_epsg(32617)
gdf_edges_utm = gdf_edges.to_crs(utm17N_prj)
gdf_DCFC_utm = gdf_DCFC.to_crs(utm17N_prj)
#Project to Web Mercator (for plotting)
wm_prj = CRS.from_epsg(3857)
gdf_edges_wm = gdf_edges.to_crs(wm_prj)
gdf_DCFC_wm = gdf_DCFC.to_crs(wm_prj)
#Plot the data
ax = gdf_edges_wm.plot(figsize=(20,10))
gdf_DCFC_wm.plot(color='red',ax=ax)
ctx.add_basemap(ax)
First, identify all "safe" areas, i.e., all alreas within a 1/2 the range of a typical EV from a DCFC charging location. We chose 1/2 the range because beyond that, the car is effectively beyond its ability to get to (i.e. return to) a charger. In the code below, we assume the car has a range of 100 miles.
Code Source: https://github.com/gboeing/osmnx-examples/blob/master/notebooks/13-isolines-isochrones.ipynb
#Set the range
range_miles = 100
#Convert the distance to meters and halve
range_meters = range_miles * 1609.344 / 2
#Create a list to hold the subgraphs created
subgraphs = []
#Iterate through each row and compute a subgraph
for i, row in gdf_DCFC.iterrows():
#Status
print('.',end='')
#Get the coordinates
thePoint = (row.geometry.y, row.geometry.x)
#Get the ID
theID = row.id
#Get the nearest node
theNode = ox.get_nearest_node(nc_graph,thePoint)
#Get the subgraph
subgraph = nx.ego_graph(G=nc_graph, n=theNode, radius=range_meters, distance='length')
#Convert to a geodataframe
#subgdf = ox.graph_to_gdfs(G=subgraph,edges=True, nodes=False)
#subgdf['subgraph_id'] = i
#Add to list
subgraphs.append(subgraph)
#Status
print(f"\n{len(subgraphs)} added to subgraphs variable")
#Save subgraphs to graphml files
if not os.path.exists('./Data/subgraphs'):
os.mkdir('./Data/subgraphs')
for i,subgraph in enumerate(subgraphs):
ox.save_graphml(G=subgraph,filename='gml_{}'.format(i),folder='./Data/subgraphs/')
ox.save_graph_shapefile?
#Merge the subgraphs gdfs together
gdfSubgraphs = pd.concat(subgraphs)
#Get a list of end nodes
nc_endnodes = [n for n in nc_graph.nodes if nc_graph.out_degree(n)==0]
#Create a mask
nc_end_mask = gdf_nodes.osmid.isin(nc_endnodes)
#Subset end nodes in the gdf
gdf_nodes.loc[gdf_nodes.osmid.isin(nc_endnodes)].to_file('Data/NC_EndNodes.shp')
#Find the end nodes
local_endnodes = [n for n in subgraph.nodes if subgraph.out_degree(n)==0]
#Create a mask
local_end_mask = gdf_nodes.osmid.isin(endnodes)
#Make mask of only local end nodes
set_nc = set(nc_endnodes)
set_local = set(local_endnodes)
print(len(set_nc),len(set_local))
set_local_end = set_local.difference(set_nc)
len(set_local_end)
gdf_end = gdf_nodes.loc[gdf_nodes.osmid.isin(set_local_end)]
gdf_end.to_file("Data\\EndNodes2.shp")
gdf_nodes, gdf_edges = ox.graph_to_gdfs(subgraph)
gdf_nodes['outdeg'] = gdf_nodes.osmid.apply(lambda x: subgraph.out_degree(x))
gdf_nodes.to_file('Data/foo.shp')
gdfSubgraphs['type'] = gdfSubgraphs.highway.apply(lambda x: str(x))
#Save as a shapefile
gdfSubgraphs[['subgraph_id','length','type','geometry']].to_file('Data\\Subgraphs.shp')
#Plot the entire network
fig, ax = ox.plot_graph(nc_graph, fig_height=10,show=False, edge_color='k', edge_alpha=0.2, node_color='none',close=False)
#Add the subgraphs, in red
gdfSubgraphs.plot(ax=ax,color='red')
#Add the DCFC sites
gdf_DCFC.plot(ax=ax, color='blue')
The above figure reveals where a car with 100 mile range could drive. To increase this range, we'd need to add a charger anywhere within 50 miles of the existing "safe zone". To do that, we need to find all the terminal nodes
#Get the coordinates
theRow = gdf_DCFC.iloc[6]
thePoint = (theRow.geometry.y, theRow.geometry.x)
#Get the ID
theID = row.id
#Get the nearest node
theNode = ox.get_nearest_node(nc_graph,thePoint)
#Get the subgraph
subgraph = nx.ego_graph(G=nc_graph, n=theNode, radius=range_meters, distance='length')
the_map = ox.plot_graph_folium(nc_graph)
ox.plot_graph_folium(subgraph,graph_map=the_map)
the_map
ax = subgraphs[6].to_crs(wm_prj).plot()
ctx.add_basemap(ax=ax)
gdfSubgraphs['hwy'] = gdfSubgraphs['highway'].apply(lambda x: x[0])
gdfSubgraphs[['length','hwy','geometry']].to_file('Data\\Subgraphs.shp')
ox.get_node(2706300546)
nx.edge_boundary(subgraph)
out_ids = []
out_nbrs = []
out_geoms = []
for n,data in subgraph.nodes(data=True):
if subgraph.out_degree(n)==0 and nc_graph.out_degree(n)==0:
out_ids.append(data['osmid'])
out_geoms.append(Point(data['x'],data['y']))
gdfTerminals = gpd.GeoDataFrame()
gdfTerminals['osmid'] = out_ids
gdfTerminals['geometry'] = out_geoms
gdfTerminals.shape
gdfTerminals.to_file('Data/TerminalPts.shp')
gNodes,gEdges = ox.graph_to_gdfs(nc_graph)
def fixIt(val):
if type(val) == list:
return val[0]
else:
return val
gEdges['hwy'] = gEdges.highway.apply(fixIt)
gEdges[['hwy','geometry']].to_file('Data/FixedEdges.shp')