exercise_test_network.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise for OR-INF transport module 2021\n",
    "Introduction to the basics of transport modelling. If you want to learn more by yourself, have a read in [Modelling Transport by Ortúzar and Willumsen](https://www.wiley.com/en-us/Modelling+Transport%2C+4th+Edition-p-9780470760390) or have a look into a more complex [transport model for Germany](https://github.com/marlinarnz/quetzal_germany)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Libraries we need\n",
    "import pandas as pd\n",
    "import geopandas as gpd\n",
    "import numpy as np\n",
    "from shapely import geometry\n",
    "from geopy.distance import geodesic"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `pandas` and `numpy` libraries are core tools for all kinds of data handling and analysis in Python. `pandas` allows easy and quick handling of data in so-called DataFrames ([pandas tutorials](https://pandas.pydata.org/pandas-docs/stable/getting_started/intro_tutorials/index.html)). `geopandas` is an extension of `pandas` which adds functionalities for geometric operations. `numpy` is the essential library for mathematical operations. `shapely.geometry` provides classes for geometric operations, while `geopy.distance.geodesic` provides a method for calculation geometric distances in meters."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create a test network\n",
    "Before we get startet we need nodes and links to carry out transport analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create our test network nodes\n",
    "nodes = pd.DataFrame(\n",
    "    data={'name': ['Berlin', 'Hamburg', 'Frankfurt', 'Stuttgart'],\n",
    "          'lat': [52.518611, 53.550556, 50.110556, 48.775556],\n",
    "          'lon': [13.408333, 9.993333, 8.682222, 9.182778]})\n",
    "nodes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a GeoDataFrame. It looks similar to the normal\n",
    "# DataFrame but has additional attributes like plotting\n",
    "nodes = gpd.GeoDataFrame(\n",
    "    nodes, geometry=gpd.points_from_xy(nodes.lon, nodes.lat))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load basemap shapes\n",
    "countries = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))\n",
    "ger = countries[countries['name']=='Germany']\n",
    "# Plot test network nodes\n",
    "nodes.plot(c='r', ax=ger.plot())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create links from every node to every node for each mode\n",
    "links = pd.DataFrame(\n",
    "    {'a': [0,0,0,1,1,1,2,2,2,3,3,3]*2,\n",
    "     'b': [1,2,3,0,2,3,0,1,3,0,1,2]*2,\n",
    "     'mode': ['rail']*12+['car']*12})\n",
    "# Add the geometry attributes\n",
    "links = gpd.GeoDataFrame(\n",
    "    links, geometry=[geometry.LineString([\n",
    "        nodes.loc[row.a, 'geometry'],nodes.loc[row.b, 'geometry']\n",
    "    ]) for _, row in links.iterrows()])\n",
    "# Add the distance in km\n",
    "links['dist'] = [geodesic(l.coords[0], l.coords[-1]).km\n",
    "                 for l in list(links['geometry'])]\n",
    "# Add level of service (LoS) attribute time in h\n",
    "links['time'] = links['dist'] / 110 #km/h is about average\n",
    "# Adjust time of rail services\n",
    "links.loc[links['mode']=='rail', 'time'] = [\n",
    "    1.8,4.2,6.2, 1.8,3.5,5.5, 4.2,3.5,1.5, 6.2,5.5,1.5]\n",
    "# Sort values by OD-pairs\n",
    "links = links.sort_values(['a', 'b']).reset_index(drop=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# In our links table, every two rows are one OD pair now\n",
    "links"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot\n",
    "nodes.plot(c='r', ax=links.plot(color='y', ax=ger.plot()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 1: Calculate mode choice probabilities\n",
    "Now we can calculate the probabilities of using one or the other mode in a logit formulation. First, we use the level of service (LoS) attribute time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Mathematical formulation of the MNL model\n",
    "def mnl(od_pair, los_col_names=['time'], beta=[1.5]):\n",
    "    '''Calculates the probabilities for each mode on the given\n",
    "        OD pair. One OD pair are all links between the same\n",
    "        origin and destination (all modes of analysis)\n",
    "    @param od_pair: pandas DataFrame slice of links from one \n",
    "        origin to one destination for all modes of analysis\n",
    "    @param los_col_names: list of DataFrame column names for all\n",
    "        LoS attributes of analysis (in the same order as betas)\n",
    "    @param beta: list of numerical calibration factors for \n",
    "        every LoS attribute given\n",
    "    @return: list of probabilities for each mode in od_pair\n",
    "        in order of appearence'''\n",
    "    v = [- sum([row[los_col_names[i]] * beta[i]\n",
    "                for i in range(len(los_col_names))]\n",
    "              ) for _, row in od_pair.iterrows()]\n",
    "    p = [np.exp(val) / sum([np.exp(m) for m in v])\n",
    "         for val in v]\n",
    "    return p"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use the above function to calculate a new column `probs_1` in links. Hint: Iterate over OD-pairs, give them to the MNL function and save the returned values to the new column of this OD pair\n",
    "\n",
    "With `.loc[]` you can access rows, columns and single values in tables based on their index and column names ([docs](https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.DataFrame.loc.html#pandas.DataFrame.loc)).\n",
    "With `.iloc[]` you can access rows, columns and single values in tables based on their index and column numbers ([docs](https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.DataFrame.iloc.html#pandas.DataFrame.iloc))."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# There are many ways to get probabilities as a new column\n",
    "# This might be easily understandable\n",
    "links['probs_1'] = "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "links"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate the mean of the results\n",
    "links.groupby('mode').agg({'probs_1': 'mean'})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 2: New LoS attribute travel cost\n",
    "After adding a new column for travel cost, you can apply the above formula again. First, add a column `cost` to links. Cost in EUR have the following dependencies with trip distance:\n",
    "\n",
    "- car (average direct cost in Germany without vehicle maintenance or purchase cost): 0.12 * distance\n",
    "- rail (ticket prices without Bahncard or Sparpreis): 0.7 * distance^0.8\n",
    "\n",
    "Hint: Use `np.power(base, exponent)` for calculation of rail fares"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "links['cost'] = "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "links['probs_2'] = "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "links"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate the mean of all results\n",
    "links.groupby('mode').agg({'probs_1': 'mean', 'probs_2': 'mean'})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 3: Political action for sustainable transport\n",
    "Now assume: Due to tax reductions for rail operators rail services cost 10% less. Moreover, a speed limit of 120km/h on motorways is enforced, which reduces average car travel speed to 100km/h. What are the new mode choice probabilities?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "links['probs_3'] = "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "links"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate the mean of all results\n",
    "links.groupby('mode').agg({'probs_1': 'mean',\n",
    "                           'probs_2': 'mean',\n",
    "                           'probs_3': 'mean'})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 4: Trip volumes\n",
    "Given an OD-matrix, you can calculate trip columes on each link. How many people are travelling where using which mode?\n",
    "\n",
    "What is the summed up travel time?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Number of trips per year\n",
    "od_matrix = pd.read_csv('od_matrix.csv', header=None)\n",
    "od_matrix.columns = ['a', 'b', 'volumes']\n",
    "od_matrix.set_index(['a', 'b'], inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "od_matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use a dictionary to translate numbers in columns \n",
    "# a and b of links table to names in the OD matrix\n",
    "name_dict = nodes['name'].to_dict()\n",
    "# Access items as follows\n",
    "name_dict[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate trip volumes\n",
    "links['volumes'] = "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "links.head(4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate passenger kilometer (pkm)\n",
    "# You can use all mathematical operations on whole columns\n",
    "links['pkm'] = "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "links.groupby('mode').agg({'pkm': 'mean'})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the volumes on the network as width of links\n",
    "# for road transport\n",
    "nodes.plot(c='r', ax=links.loc[links['mode']=='car'].plot(\n",
    "    color='y', linewidth=[v / links['volumes'].mean()\n",
    "                          for v in links['volumes']], ax=ger.plot()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the volumes on the network as width of links\n",
    "# for rail transport\n",
    "nodes.plot(c='r', ax=links.loc[links['mode']=='rail'].plot(\n",
    "    color='y', linewidth=[v / links['volumes'].mean()\n",
    "                          for v in links['volumes']], ax=ger.plot()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Aggregated travel time\n",
    "t = \n",
    "t.sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 5: Capacities and congestion\n",
    "So far, travel on our test network was not restricted. However, in reality roads become congested and trains crowded. Given the capacities of links and a volume-dependent congestion formula, how do people travel now? What is the new aggregated travel time?\n",
    "\n",
    "If you want, you can try multiple iterations to reach an equilibrium between volumes and capacities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Average occupancy in ICE trains is low in Germany while\n",
    "# main motorways are quite packed\n",
    "links['cap'] = [od_matrix.loc[\n",
    "    (name_dict[row['a']], name_dict[row['b']]), 'volumes'] / 2\n",
    "#                if row['mode']=='car' else row['volumes']*3\n",
    "                for _, row in links.iterrows()]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def travel_time(volumes, capacity, time):\n",
    "    '''Calculates the travel time of a link with capacity\n",
    "    restriction applying the formula used by Smock (1962)\n",
    "    @param volumes: volume of trips on this link\n",
    "    @param capacity: capacity for trips on this link\n",
    "    @param time: travel time without capacity restriction\n",
    "    @return: new travel time'''\n",
    "    return 0.5 * time * np.exp(volumes / capacity)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "links['time_2'] = [travel_time(row['volumes'], row['cap'], row['time'])\n",
    "                   if row['mode']=='car' else row['time']\n",
    "                   for _, row in links.iterrows()]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "links.head(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "links['volumes_2'] = np.nan\n",
    "links['probs_4'] = np.nan"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for j in range(10):\n",
    "    # Probabilities\n",
    "    \n",
    "    # New volumes\n",
    "    links['volumes_2'] = \n",
    "    # New travel times\n",
    "    links['time_2'] = \n",
    "    print('Iteration {} av. mode shares: car {}; rail {}'.format(\n",
    "        j, links.loc[links['mode']=='car', 'probs_4'].mean(),\n",
    "        links.loc[links['mode']=='rail', 'probs_4'].mean()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "links.sample(4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate the mean of all results\n",
    "links.groupby('mode').agg({'probs_1': 'mean',\n",
    "                           'probs_2': 'mean',\n",
    "                           'probs_3': 'mean',\n",
    "                           'probs_4': 'mean'})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Aggregated travel time\n",
    "# Compare to step 5\n",
    "t = \n",
    "t.sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 6: Impact of lorries\n",
    "Looking into the future, roads will become even more congested by lorries due to the rising demand for goods transport, if there is no political ambition for a shift to rail in freight transport. Given these new capacities, how does the mode share in passenger transport change? And the aggregated travel time?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "links['cap_2'] = [row['cap']*0.9\n",
    "                  if row['mode']=='car' else row['cap']\n",
    "                  for _, row in links.iterrows()]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# End\n",
    "Congrats! You've learned the basics of open source transport modelling with Python. If you want you can go on with altering capacities, adding new nodes or modes or even adding a new LoS attribute."
   ]
  }
 ],
 "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.8.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}