Tutorial 2.1 - Shortest path analysis#

Attention

Finnish university students are encouraged to use the CSC Notebooks platform.
CSC badge

Others can follow the lesson interactively using Binder. Check the rocket icon on the top of this page.

Lesson objectives

This tutorial focuses on spatial networks and learn how to construct a routable directed graph for Networkx and find shortest paths along the given street network based on travel times or distance by car. In addition, we will learn how to calculate travel times from a single source into all nodes in the graph.

Tutorial#

In this tutorial we will focus on a network analysis methods that relate to way-finding. Finding a shortest path from A to B using a specific street network is a very common spatial analytics problem that has many practical applications.

Python provides easy to use tools for conducting spatial network analysis. One of the easiest ways to start is to use a library called Networkx which is a Python module that provides a lot tools that can be used to analyze networks on various different ways. It also contains algorithms such as Dijkstra’s algorithm or A* algoritm that are commonly used to find shortest paths along transportation network.

Next, we will learn how to do spatial network analysis in practice.

Typical workflow for routing#

If you want to conduct network analysis (in any programming language) there are a few basic steps that typically needs to be done before you can start routing. These steps are:

  1. Retrieve data (such as street network from OSM or Digiroad + possibly transit data if routing with PT).

  2. Modify the network by adding/calculating edge weights (such as travel times based on speed limit and length of the road segment).

  3. Build a routable graph for the routing tool that you are using (e.g. for NetworkX, igraph or OpenTripPlanner).

  4. Conduct network analysis (such as shortest path analysis) with the routing tool of your choice.

1. Retrieve data#

As a first step, we need to obtain data for routing. Pyrosm library makes it really easy to retrieve routable networks from OpenStreetMap (OSM) with different transport modes (walking, cycling and driving).

  • Let’s first extract OSM data for Helsinki that are walkable. In pyrosm, we can use a function called osm.get_network() which retrieves data from OpenStreetMap. It is possible to specify what kind of roads should be retrieved from OSM with network_type -parameter (supports walking, cycling, driving).

from pyrosm import OSM, get_data
import geopandas as gpd
import pandas as pd
import networkx as nx

# We will use test data for Helsinki that comes with pyrosm
osm = OSM(get_data("helsinki_pbf"))

# Parse roads that can be driven by car
roads = osm.get_network(network_type="driving")
roads.plot(figsize=(10,10))
<AxesSubplot: >
../../_images/e8f456f84bae89dbb4ec3e037732c70f5b6afa95f0bdbdd202720e085f7bc4fd.png
roads.head(2)
access area bicycle bridge cycleway foot footway highway int_ref lanes ... surface tunnel width id timestamp version tags osm_type geometry length
0 None None None None None None None unclassified None 2 ... paved None None 4236349 1380031970 21 {"name:fi":"Erottajankatu","name:sv":"Skillnad... way MULTILINESTRING ((24.94327 60.16651, 24.94337 ... 14.0
1 None None None None None None None unclassified None 2 ... paved None None 4243035 1543430213 12 {"name:fi":"Korkeavuorenkatu","name:sv":"H\u00... way MULTILINESTRING ((24.94567 60.16767, 24.94567 ... 51.0

2 rows × 30 columns

Okay, now we have drivable roads as a GeoDataFrame for the city center of Helsinki. If you look at the GeoDataFrame (scroll to the right), we can see that pyrosm has also calculated us the length of each road segment (presented in meters). The geometries are presented here as MultiLineString objects. From the map above we can see that the data also includes short pieces of roads that do not lead to anywhere (i.e. they are isolated). This is a typical issue when working with real-world data such as roads. Hence, at some point we need to take care of those in someway (remove them (typical solution), or connect them to other parts of the network).

In OSM, the information about the allowed direction of movement is stored in column oneway. Let’s take a look what kind of values we have in that column:

roads["oneway"].unique()
array(['yes', None, 'no'], dtype=object)

As we can see the unique values in that column are "yes", "no" or None. We can use this information to construct a directed graph for routing by car. For walking and cycling, you typically want create a bidirectional graph, because the travel is typically allowed in both directions at least in Finland. Notice, that the rules vary by country, e.g. in Copenhagen you have oneway rules also for bikes but typically each road have the possibility to travel both directions (you just need to change the side of the road if you want to make a U-turn). Column maxspeed contains information about the speed limit for given road:

roads["maxspeed"].unique()
array(['30', '40', None, '20', '10', '5', '50'], dtype=object)

As we can see, there are also None values in the data, meaning that the speed limit has not been tagged for some roads. This is typical, and often you need to fill the non existing speed limits yourself. This can be done by taking advantage of the road class that is always present in column highway:

roads["highway"].unique()
array(['unclassified', 'residential', 'secondary', 'service', 'tertiary',
       'primary', 'primary_link', 'cycleway', 'footway', 'tertiary_link',
       'pedestrian', 'trail', 'crossing'], dtype=object)

Based on these values, we can make assumptions that e.g. residential roads in Helsinki have a speed limit of 30 kmph. Hence, this information can be used to fill the missing values in maxspeed. As we can see, the current version of the pyrosm tool seem to have a bug because some non-drivable roads were also leaked to our network (e.g. footway, cycleway). If you notice these kind of issues with any of the libraries that you use, please notify the developers by raising an Issue in GitHub. This way, you can help improving the software. For this given problem, an issue has already been raised so you don’t need to do it again (it’s always good to check if a related issue exists in GitHub before adding a new one).

Okay, but how can we make a routable graph out of this data of ours? Let’s remind us about the basic elements of a graph that we went through in the lecture slides:

Basic elements of a graph

So to be able to create a graph we need to have nodes and edges. Now we have a GeoDataFrame of edges, but where are those nodes? Well they are not yet anywhere, but with pyrosm we can easily retrieve the nodes as well by specifying nodes=True, when parsing the streets:

# Parse nodes and edges
nodes, edges = osm.get_network(network_type="driving", nodes=True)

# Plot the data
ax = edges.plot(figsize=(10,10), color="gray", lw=1.0)
ax = nodes.plot(ax=ax, color="red", markersize=2)

# Zoom in to take a closer look
#ax.set_xlim([24.9375, 24.945])
ax.set_ylim([60.17, 60.173])
(60.17, 60.173)
../../_images/27361d1a39eb2347e314578e3a2d32f56d0225faf1850afa0e4e920324f35d2d.png

Okay, as we can see now we have both the roads (i.e. edges) and the nodes that connect the street elements together (in red) that are typically intersections. However, we can see that many of the nodes are in locations that are clearly not intersections. This is intented behavior to ensure that we have full connectivity in our network. We can at later stage clean and simplify this network by merging all roads that belong to the same link (i.e. street elements that are between two intersections) which also reduces the size of the network.

Note

In OSM, the street topology is typically not directly suitable for graph traversal due to missing nodes at intersections which means that the roads are not splitted at those locations. The consequence of this, is that it is not possible to make a turn if there is no intersection present in the data structure. Hence, pyrosm will separate all road segments/geometries into individual rows in the data.

Let’s take a look what our nodes data look like:

nodes.head()
lon lat tags timestamp version changeset id geometry
0 24.943271 60.166514 None 1390926206 2 0 1372477605 POINT (24.94327 60.16651)
1 24.943365 60.166444 {'highway': 'crossing', 'crossing': 'traffic_s... 1383915357 6 0 292727220 POINT (24.94337 60.16644)
2 24.943403 60.166408 None 1374595731 1 0 2394117042 POINT (24.94340 60.16641)
3 24.945668 60.167668 {'highway': 'crossing', 'crossing': 'uncontrol... 1290714658 5 0 296250563 POINT (24.94567 60.16767)
4 24.945671 60.167630 {'traffic_calming': 'divider'} 1354578076 1 0 2049084195 POINT (24.94567 60.16763)

As we can see, the nodes GeoDataFrame contains information about the coordinates of each node as well as a unique id for each node. These id values are used to determine the connectivity in our network. Hence, pyrosm has also added two columns to the edges GeoDataFrame that specify from and to ids for each edge. Column u contains information about the from-id and column v about the to-id accordingly:

# Check last four columns
edges.iloc[:5,-4:]
geometry u v length
0 LINESTRING (24.94327 60.16651, 24.94337 60.16644) 1372477605 292727220 9.370
1 LINESTRING (24.94337 60.16644, 24.94340 60.16641) 292727220 2394117042 4.499
2 LINESTRING (24.94567 60.16767, 24.94567 60.16763) 296250563 2049084195 4.174
3 LINESTRING (24.94567 60.16763, 24.94569 60.16744) 2049084195 60072359 21.692
4 LINESTRING (24.94569 60.16744, 24.94571 60.16726) 60072359 6100704327 19.083

We can see that the geometries are now stored as LineString instead of MultiLineString. At this point, we can fix the issue related to having some pedestrian roads in our network. We can do this by removing all edges from out GeoDataFrame that have highway value in 'cycleway', 'footway', 'pedestrian', 'trail', 'crossing':

edges = edges.loc[~edges["highway"].isin(['cycleway', 'footway', 'pedestrian', 'trail', 'crossing'])].copy()
edges.plot()
<AxesSubplot: >
../../_images/0aa396788145f658d7b659a6bb8de99c75dc650f318dc31fb9b566346809eec6.png

Now we can see, that some of the isolated edges were removed from the data. The character ~ (tilde) in the command above is a negation operator that is handy if you want to e.g. remove some rows from your GeoDataFrame based on criteria such as we used here.

2. Modify the data#

At this stage, we have the necessary components to build a routable graph (nodes and edges) based on distance. However, in real life the network distance is not the best cost metric to use, because the shortest path (based on distance) is not necessarily always the optimal route in terms of travel time. Time is typically the measure that people value more (plus it is easier to comprehend), so at this stage we want to add a new cost attribute to our edges GeoDataFrame that converts the metric distance information to travel time (in seconds) based on following formula:

  • <distance-in-meters> / (<speed-limit-kmph> / 3.6)

Before we can do this calculation, we need to ensure that all rows in maxspeed column have information about the speed limit. Let’s check the value counts of the column and also include information about the NaN values with dropna parameter:

# Count values
edges["maxspeed"].value_counts(dropna=False)
30      1110
None     628
40       422
10        50
20        36
5         21
50         2
Name: maxspeed, dtype: int64

As we can see, the rows which do not contain information about the speed limit is the second largest group in our data. Hence, we need to apply a criteria to fill these gaps. We can do this based on following “rule of thumb” criteria in Finland (notice that these vary country by country):

Road class

Speed limit within urban region

Speed limit outside urban region

motorway

100

120

motorway_link

80

80

trunk

60

100

trunk_link

60

60

primary

50

80

primary_link

50

50

secondary

50

50

secondary_link

50

50

tertiary

50

60

tertiary_link

50

50

unclassified

50

80

unclassified_link

50

50

residential

50

80

living_street

20

NA

service

30

NA

other

50

80

For simplicity, we can consider that all the roads in Helsinki Region follows the within urban region speed limits, although this is not exactly true (the higher speed limits start somewhere at the outer parts of the city region). For making the speed limit values more robust / correct, you could use data about urban/rural classification which is available in Finland from Finnish Environment Institute. Let’s first convert our maxspeed values to integers using astype() method:

edges["maxspeed"] = edges["maxspeed"].astype(float).astype(pd.Int64Dtype())
edges["maxspeed"].unique()
<IntegerArray>
[30, 40, <NA>, 20, 10, 5, 50]
Length: 7, dtype: Int64

As we can see, now the maxspeed values are stored in integer format inside an IntegerArray, and the None values were converted into pandas.NA objects that are assigned with <NA>. Now we can create a function that returns a numeric value for different road classes based on the criteria in the table above:

def road_class_to_kmph(road_class):
    """
    Returns a speed limit value based on road class, 
    using typical Finnish speed limit values within urban regions.
    """
    if road_class == "motorway":
        return 100
    elif road_class == "motorway_link":
        return 80
    elif road_class in ["trunk", "trunk_link"]:
        return 60
    elif road_class == "service":
        return 30
    elif road_class == "living_street":
        return 20
    else:
        return 50

Now we can apply this function to all rows that do not have speed limit information:

# Separate rows with / without speed limit information 
mask = edges["maxspeed"].isnull()
edges_without_maxspeed = edges.loc[mask].copy()
edges_with_maxspeed = edges.loc[~mask].copy()

# Apply the function and update the maxspeed
edges_without_maxspeed["maxspeed"] = edges_without_maxspeed["highway"].apply(road_class_to_kmph)
edges_without_maxspeed.head(5).loc[:, ["maxspeed", "highway"]]
maxspeed highway
30 30 service
47 30 service
48 30 service
49 30 service
50 30 service

Okay, as we can see now the maxspeed value have been updated according our criteria, and e.g. the service road class have been given the speed limit 30 kmph. Now we can recreate the edges GeoDataFrame by combining the two frames:

edges = edges_with_maxspeed.append(edges_without_maxspeed)
edges["maxspeed"].unique()
/tmp/ipykernel_385531/2154535046.py:1: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
  edges = edges_with_maxspeed.append(edges_without_maxspeed)
<IntegerArray>
[30, 40, 20, 10, 5, 50]
Length: 6, dtype: Int64

Great, now all of our edges have information about the speed limit. We can also visualize them:

# Convert the value into regular integer Series (the plotting requires having Series instead of IntegerArray) 
edges["maxspeed"] = edges["maxspeed"].astype(int)
ax = edges.plot(column="maxspeed", figsize=(10,10), legend=True)
../../_images/c38433408562e780e045c9a932cde54e84507f9d66a2db32d9dae4434f3b2cb4.png

Finally, we can calculate the travel time in seconds using the formula we saw earlier and add that as a new cost attribute for our network:

edges["travel_time_seconds"] = edges["length"] / (edges["maxspeed"]/3.6)
edges.iloc[0:10, -4:]
u v length travel_time_seconds
0 1372477605 292727220 9.370 1.12440
1 292727220 2394117042 4.499 0.53988
2 296250563 2049084195 4.174 0.50088
3 2049084195 60072359 21.692 2.60304
4 60072359 6100704327 19.083 2.28996
5 6100704327 296250223 6.027 0.72324
6 264015226 25345665 9.644 1.15728
7 25345665 296248024 7.016 0.84192
8 296248024 426911766 4.137 0.49644
9 426911766 60072364 21.132 2.53584

Excellent! Now our GeoDataFrame has all the information we need for creating a graph that can be used to conduct shortest path analysis based on length or travel time. Notice that here we assume that the cars can drive with the same speed as what the speed limit is. Considering the urban dynamics and traffic congestion, this assumption might not hold, but for simplicity, we assume so in this tutorial.

3. Build a directed graph for routing using pyrosm#

Now as we have calculated the travel time for our edges. We still need to convert our nodes and edges into a directed graph, so that we can start using it for routing. There are easy-to-use functionalities for doing this in pyrosm and osmnx.

As a background information, it is good to understand that our edges represents a directed network. This means that the information stored in oneway column will be used to determine the network structure for the edges based on the rules in that column. If the oneway is 'yes', it means that the street can be driven only to one direction, and if it is None or has a value "no", then that road can be driven to both directions. This means that the tool will make new duplicate edge and reversing the from-id and to-id values in the u and v columns in the data. In addition, value -1 in the oneway column means that the road can only be driven to one direction but against the digitization direction. In such cases, the edge is flipped: the to-id (u) and from-id (u) is reversed so that the directionality in the graph is correctly specified.

  • Let’s see how we can create a routable NetworkX graph using pyrosm with one command:

G = osm.to_graph(nodes, edges, graph_type="networkx")
G
<networkx.classes.multidigraph.MultiDiGraph at 0x7fd04cb39eb0>

Now we have a routable graph. pyrosm actually does some additional steps in the background. By default, pyrosm cleans all unconnected edges from the graph and only keeps edges that can be reached from every part of the network. In addition, pyrosm automatically modifies the graph attribute information in a way that they are compatible with OSMnx that provides many handy functionalities to work with graphs. Such as plotting an interactive map based on the graph:

import osmnx as ox 
ox.plot_graph_folium(G)
Make this Notebook Trusted to load map: File -> Trust Notebook