{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Analyze data\n",
    "* Read in the saved DCFC and graph of NC major roads.\n",
    "* Analyze graph to locate all roads within 1/2 of vehicle's range of each DCFC (\"Safe areas\")\n",
    "* Analyze graph to locate all roads within vehicle's full range of each DCFC (\"Full range\")\n",
    "* Spatially subtract \"Safe areas\" from \"full range\" (\"Anxious areas\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Import packages\n",
    "import os\n",
    "import osmnx as ox\n",
    "import networkx as nx\n",
    "import pandas as pd\n",
    "import geopandas as gpd\n",
    "from shapely.geometry import Point\n",
    "import contextily as ctx\n",
    "import matplotlib.pyplot as plt\n",
    "from pyproj import CRS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read in the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Filenames\n",
    "DCFC_filename = 'Data\\\\NREL\\\\NC_DCFC.shp'\n",
    "graph_filename = 'Data\\\\OSM\\\\NC_highways_all.graphml'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Read DCFC points into geodataframe\n",
    "print(f\"Loading DCFC points from {DCFC_filename}\")\n",
    "gdf_DCFC = gpd.read_file(DCFC_filename)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Read graphML file into a graph object\n",
    "print(f\"Loading graph from {graph_filename}\")\n",
    "nc_graph = ox.load_graphml(filename=os.path.basename(graph_filename),\n",
    "                           folder=os.path.dirname(graph_filename))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Convert graph component into geodataframes\n",
    "gdf_nodes, gdf_edges = ox.graph_to_gdfs(nc_graph)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Project to UTM (for analysis)\n",
    "utm17N_prj = CRS.from_epsg(32617)\n",
    "gdf_edges_utm = gdf_edges.to_crs(utm17N_prj)\n",
    "gdf_DCFC_utm = gdf_DCFC.to_crs(utm17N_prj)\n",
    "\n",
    "#Project to Web Mercator (for plotting)\n",
    "wm_prj = CRS.from_epsg(3857)\n",
    "gdf_edges_wm = gdf_edges.to_crs(wm_prj)\n",
    "gdf_DCFC_wm = gdf_DCFC.to_crs(wm_prj)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualize the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Plot the data\n",
    "ax = gdf_edges_wm.plot(figsize=(20,10))\n",
    "gdf_DCFC_wm.plot(color='red',ax=ax)\n",
    "ctx.add_basemap(ax)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Identify \"Safe Areas\"\n",
    "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. \n",
    "\n",
    "Code Source: https://github.com/gboeing/osmnx-examples/blob/master/notebooks/13-isolines-isochrones.ipynb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Set the range\n",
    "range_miles = 100\n",
    "\n",
    "#Convert the distance to meters and halve \n",
    "range_meters = range_miles * 1609.344 / 2\n",
    "\n",
    "#Create a list to hold the subgraphs created\n",
    "subgraphs = []\n",
    "\n",
    "#Iterate through each row and compute a subgraph\n",
    "for i, row in gdf_DCFC.iterrows():\n",
    "    #Status\n",
    "    print('.',end='')\n",
    "    #Get the coordinates\n",
    "    thePoint = (row.geometry.y, row.geometry.x)\n",
    "    #Get the ID\n",
    "    theID = row.id\n",
    "    #Get the nearest node\n",
    "    theNode = ox.get_nearest_node(nc_graph,thePoint)\n",
    "    #Get the subgraph\n",
    "    subgraph = nx.ego_graph(G=nc_graph, n=theNode, radius=range_meters, distance='length')\n",
    "    #Convert to a geodataframe\n",
    "    #subgdf = ox.graph_to_gdfs(G=subgraph,edges=True, nodes=False)\n",
    "    #subgdf['subgraph_id'] = i\n",
    "    #Add to list\n",
    "    subgraphs.append(subgraph)\n",
    "\n",
    "#Status\n",
    "print(f\"\\n{len(subgraphs)} added to subgraphs variable\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Save subgraphs to graphml files\n",
    "if not os.path.exists('./Data/subgraphs'):\n",
    "    os.mkdir('./Data/subgraphs')\n",
    "for i,subgraph in enumerate(subgraphs):\n",
    "    ox.save_graphml(G=subgraph,filename='gml_{}'.format(i),folder='./Data/subgraphs/')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ox.save_graph_shapefile?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Merge the subgraphs gdfs together\n",
    "gdfSubgraphs = pd.concat(subgraphs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Find end nodes in full graph (to separate from subgraph end nodes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Get a list of end nodes\n",
    "nc_endnodes = [n for n in nc_graph.nodes if nc_graph.out_degree(n)==0]\n",
    "#Create a mask\n",
    "nc_end_mask = gdf_nodes.osmid.isin(nc_endnodes)\n",
    "#Subset end nodes in the gdf\n",
    "gdf_nodes.loc[gdf_nodes.osmid.isin(nc_endnodes)].to_file('Data/NC_EndNodes.shp')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Find the end nodes\n",
    "local_endnodes = [n for n in subgraph.nodes if subgraph.out_degree(n)==0]\n",
    "#Create a mask\n",
    "local_end_mask = gdf_nodes.osmid.isin(endnodes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Make mask of only local end nodes\n",
    "set_nc = set(nc_endnodes)\n",
    "set_local = set(local_endnodes)\n",
    "print(len(set_nc),len(set_local))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "set_local_end = set_local.difference(set_nc)\n",
    "len(set_local_end)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gdf_end = gdf_nodes.loc[gdf_nodes.osmid.isin(set_local_end)]\n",
    "gdf_end.to_file(\"Data\\\\EndNodes2.shp\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gdf_nodes, gdf_edges = ox.graph_to_gdfs(subgraph)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gdf_nodes['outdeg'] = gdf_nodes.osmid.apply(lambda x: subgraph.out_degree(x))\n",
    "gdf_nodes.to_file('Data/foo.shp')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gdfSubgraphs['type'] = gdfSubgraphs.highway.apply(lambda x: str(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Save as a shapefile\n",
    "gdfSubgraphs[['subgraph_id','length','type','geometry']].to_file('Data\\\\Subgraphs.shp')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Plot the entire network\n",
    "fig, ax = ox.plot_graph(nc_graph, fig_height=10,show=False, edge_color='k', edge_alpha=0.2, node_color='none',close=False)\n",
    "#Add the subgraphs, in red\n",
    "gdfSubgraphs.plot(ax=ax,color='red')\n",
    "#Add the DCFC sites\n",
    "gdf_DCFC.plot(ax=ax, color='blue')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Next step:\n",
    "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 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Get the coordinates\n",
    "theRow = gdf_DCFC.iloc[6]\n",
    "thePoint = (theRow.geometry.y, theRow.geometry.x)\n",
    "#Get the ID\n",
    "theID = row.id\n",
    "#Get the nearest node\n",
    "theNode = ox.get_nearest_node(nc_graph,thePoint)\n",
    "#Get the subgraph\n",
    "subgraph = nx.ego_graph(G=nc_graph, n=theNode, radius=range_meters, distance='length')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "the_map = ox.plot_graph_folium(nc_graph)\n",
    "ox.plot_graph_folium(subgraph,graph_map=the_map)\n",
    "the_map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = subgraphs[6].to_crs(wm_prj).plot()\n",
    "ctx.add_basemap(ax=ax)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gdfSubgraphs['hwy'] = gdfSubgraphs['highway'].apply(lambda x: x[0])\n",
    "gdfSubgraphs[['length','hwy','geometry']].to_file('Data\\\\Subgraphs.shp')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ox.get_node(2706300546)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nx.edge_boundary(subgraph)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "out_ids = []\n",
    "out_nbrs = []\n",
    "out_geoms = []\n",
    "for n,data in subgraph.nodes(data=True):\n",
    "    if subgraph.out_degree(n)==0 and nc_graph.out_degree(n)==0:\n",
    "        out_ids.append(data['osmid'])\n",
    "        out_geoms.append(Point(data['x'],data['y']))\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gdfTerminals = gpd.GeoDataFrame()\n",
    "gdfTerminals['osmid'] = out_ids\n",
    "gdfTerminals['geometry'] = out_geoms\n",
    "gdfTerminals.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gdfTerminals.to_file('Data/TerminalPts.shp')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gNodes,gEdges = ox.graph_to_gdfs(nc_graph)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fixIt(val):\n",
    "    if type(val) == list:\n",
    "        return val[0]\n",
    "    else:\n",
    "        return val"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gEdges['hwy'] = gEdges.highway.apply(fixIt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gEdges[['hwy','geometry']].to_file('Data/FixedEdges.shp')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}