{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Node elevations and edge grades\n",
    "\n",
    "Author: [Geoff Boeing](https://geoffboeing.com/)\n",
    "\n",
    "  - [Overview of OSMnx](http://geoffboeing.com/2016/11/osmnx-python-street-networks/)\n",
    "  - [GitHub repo](https://github.com/gboeing/osmnx)\n",
    "  - [Examples, demos, tutorials](https://github.com/gboeing/osmnx-examples)\n",
    "  - [Documentation](https://osmnx.readthedocs.io/en/stable/)\n",
    "  - [Journal article/citation](http://geoffboeing.com/publications/osmnx-complex-street-networks/)\n",
    "\n",
    "OSMnx allows you to automatically add elevation attributes to your graph's nodes with the `elevation` module, using either local raster files or the Google Maps Elevation API as the elevation data source. If you use the Google API, you will need an API key. Once your nodes have elevation values, OSMnx can automatically calculate your edges' grades (inclines)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "\n",
    "import numpy as np\n",
    "import osmnx as ox\n",
    "import pandas as pd\n",
    "\n",
    "%matplotlib inline\n",
    "ox.__version__"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Elevation from local raster file(s)\n",
    "\n",
    "OSMnx can attach elevations to graph nodes using either a single raster file or a list of raster files. The latter creates a virtual raster VRT composed of the rasters at those filepaths. By default, it uses all available CPUs but you can configure this with an argument."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "address = \"600 Montgomery St, San Francisco, California, USA\"\n",
    "G = ox.graph_from_address(address=address, dist=500, dist_type=\"bbox\", network_type=\"bike\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# add node elevations from a single raster file\n",
    "# some nodes will be null because the single file does not cover the graph's extents\n",
    "raster_path = \"./input_data/elevation1.tif\"\n",
    "G = ox.elevation.add_node_elevations_raster(G, raster_path, cpus=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# add node elevations from multiple raster files\n",
    "# no nulls should remain\n",
    "raster_paths = [\"./input_data/elevation1.tif\", \"./input_data/elevation2.tif\"]\n",
    "G = ox.elevation.add_node_elevations_raster(G, raster_paths)\n",
    "assert not np.isnan(np.array(G.nodes(data=\"elevation\"))[:, 1]).any()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# add edge grades and their absolute values\n",
    "G = ox.elevation.add_edge_grades(G, add_absolute=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Elevation from Google Maps Elevation API\n",
    "\n",
    "You will need a Google Maps Elevation [API key](https://developers.google.com/maps/documentation/elevation/start). Consider your API usage limits. OSMnx rounds coordinates to 5 decimal places (approx 1 meter) to fit 350 locations in a batch. Note that there is some spatial inaccuracy given Google's dataset's resolution. For example, in San Francisco (where the resolution is 19 meters) a couple of edges in hilly parks have a 50+ percent grade because Google assigns one of their nodes the elevation of a hill adjacent to the street."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# replace this with your own API key!\n",
    "try:\n",
    "    from keys import google_elevation_api_key\n",
    "except ImportError:\n",
    "    sys.exit()  # you need an API key to proceed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get the street network for san francisco\n",
    "place = \"San Francisco\"\n",
    "place_query = {\"city\": \"San Francisco\", \"state\": \"California\", \"country\": \"USA\"}\n",
    "G = ox.graph_from_place(place_query, network_type=\"drive\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# add elevation to each of the nodes, using the google elevation API, then calculate edge grades\n",
    "G = ox.elevation.add_node_elevations_google(G, api_key=google_elevation_api_key)\n",
    "G = ox.elevation.add_edge_grades(G)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate some summary stats\n",
    "\n",
    "Use an undirected representation of the network so we don't overcount two-way streets (because they have reciprocal edges pointing in each direction). We use the absolute value of edge grade because we're interested in steepness, not directionality."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate the edges' absolute grades (and drop any infinite/null values)\n",
    "grades = pd.Series([d[\"grade_abs\"] for _, _, d in ox.get_undirected(G).edges(data=True)])\n",
    "grades = grades.replace([np.inf, -np.inf], np.nan).dropna()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "avg_grade = np.mean(grades)\n",
    "print(\"Average street grade in {} is {:.1f}%\".format(place, avg_grade * 100))\n",
    "\n",
    "med_grade = np.median(grades)\n",
    "print(\"Median street grade in {} is {:.1f}%\".format(place, med_grade * 100))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot the nodes by elevation\n",
    "\n",
    "Plot them colored from low (violet) to high (yellow)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get one color for each node, by elevation, then plot the network\n",
    "nc = ox.plot.get_node_colors_by_attr(G, \"elevation\", cmap=\"plasma\")\n",
    "fig, ax = ox.plot_graph(G, node_color=nc, node_size=5, edge_color=\"#333333\", bgcolor=\"k\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot the edges by grade\n",
    "\n",
    "Grade is the ratio of elevation change to edge length. Plot edges colored from low/flat (violet) to high/steep (yellow)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get a color for each edge, by grade, then plot the network\n",
    "ec = ox.plot.get_edge_colors_by_attr(G, \"grade_abs\", cmap=\"plasma\", num_bins=5, equal_size=True)\n",
    "fig, ax = ox.plot_graph(G, edge_color=ec, edge_linewidth=0.5, node_size=0, bgcolor=\"k\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate shortest paths, considering grade impedance\n",
    "\n",
    "This example approximates the route of \"The Wiggle\" in San Francisco."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# select an origin and destination node and a bounding box around them\n",
    "origin = ox.distance.nearest_nodes(G, -122.426, 37.77)\n",
    "destination = ox.distance.nearest_nodes(G, -122.441, 37.773)\n",
    "bbox = ox.utils_geo.bbox_from_point((37.772, -122.434), dist=1500)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define some edge impedance function here\n",
    "def impedance(length, grade):\n",
    "    penalty = grade**2\n",
    "    return length * penalty\n",
    "\n",
    "\n",
    "# add impedance and elevation rise values to each edge in the projected graph\n",
    "# use absolute value of grade in impedance function if you want to avoid uphill and downhill\n",
    "for _, _, _, data in G.edges(keys=True, data=True):\n",
    "    data[\"impedance\"] = impedance(data[\"length\"], data[\"grade_abs\"])\n",
    "    data[\"rise\"] = data[\"length\"] * data[\"grade\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### First find the shortest path that minimizes *trip distance*:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "route_by_length = ox.shortest_path(G, origin, destination, weight=\"length\")\n",
    "fig, ax = ox.plot_graph_route(G, route_by_length, bbox=bbox, node_size=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Now find the shortest path that avoids slopes by minimizing *impedance* (function of length and grade):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "route_by_impedance = ox.shortest_path(G, origin, destination, weight=\"impedance\")\n",
    "fig, ax = ox.plot_graph_route(G, route_by_impedance, bbox=bbox, node_size=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Print some summary stats about these two routes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def print_route_stats(route):\n",
    "    route_grades = ox.utils_graph.get_route_edge_attributes(G, route, \"grade_abs\")\n",
    "    msg = \"The average grade is {:.1f}% and the max is {:.1f}%\"\n",
    "    print(msg.format(np.mean(route_grades) * 100, np.max(route_grades) * 100))\n",
    "\n",
    "    route_rises = ox.utils_graph.get_route_edge_attributes(G, route, \"rise\")\n",
    "    ascent = np.sum([rise for rise in route_rises if rise >= 0])\n",
    "    descent = np.sum([rise for rise in route_rises if rise < 0])\n",
    "    msg = \"Total elevation change is {:.1f} meters: {:.0f} meter ascent and {:.0f} meter descent\"\n",
    "    print(msg.format(np.sum(route_rises), ascent, abs(descent)))\n",
    "\n",
    "    route_lengths = ox.utils_graph.get_route_edge_attributes(G, route, \"length\")\n",
    "    print(\"Total trip distance: {:,.0f} meters\".format(np.sum(route_lengths)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# stats of route minimizing length\n",
    "print_route_stats(route_by_length)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# stats of route minimizing impedance (function of length and grade)\n",
    "print_route_stats(route_by_impedance)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, we decreased the average slope along the route from a 5% grade to a 2% grade. The total elevation change is obviously (approximately, due to rounding) the same with either route, but using our impedance function we decrease the total ascent from 69 meters to 42 meters (but the trip distance increases from 1.9 km to 3.3 km)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python (ox)",
   "language": "python",
   "name": "ox"
  },
  "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.10.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}