Section 16 Analysing road networks with pgrouting

16.1 Introduction

PgRouting is an extension to postgis that provides functions for working with networks. It is not particularly easy to use. For small site based projects it may be simpler to use the network functions included in GRASS. pgRouting is designed for large networks held in the database with a particular focus on road networks.

A road network consists of nodes that are connected by edges, i.e sections of roads between intersections. If a road network depicted simply of a set of lines is imported into a GIS such as GRASS and then converted directly into a network some features of a real road network will be overlooked. The conversion of interecting lines into network nodes will not take into account features such as bridges, flyovers and one way streets. The open street map project aimed to collect relevant data on the transport network to allow road routing to be conducted correctly.

To run pgrouting using open street map data the first step is to import the latest osm data and convert the data into a postgis network.

16.2 Downloading the data

Running the following commands in a Linux terminal downloads the data and converts the compressed file into the format needed for importing. This assumes that osmtools are installed.

wget http://download.geofabrik.de/europe/great-britain/england/dorset-latest.osm.pbf osmconvert dorset-latest.osm.pbf > dorset.osm

The data for Hampshire can be added to that for Dorset to build up a network for a larger area.

wget http://download.geofabrik.de/europe/great-britain/england/hampshire-latest.osm.pbf osmconvert hampshire-latest.osm.pbf > hampshire.osm

16.3 Importing

Now we need to form a schema and import the data. On the BU server these commands set up the connection. If you are not on the server you would have to set up your own pgrouting enabled postgis data base and connect to that.

library(giscourse)
conn<-connect11()

The choice of mapconfig file is important when setting up a network using osm2pgrouting. The default is to load all data including tracks and paths. However for a satnav type application it makes more sense to just include roads. The mapconfig_for_cars.xml file is set up for this.

dbSendQuery(conn,"drop schema if exists pgrouting cascade; create schema pgrouting")

command<-"osm2pgrouting --clean -c mapconfig_for_cars.xml -f dorset.osm --host 172.16.49.31 --port 5433 --username gis_course --password gis_course123 --schema pgrouting --dbname gis_course"

system(command)

command<-"osm2pgrouting -c mapconfig_for_cars.xml -f hampshire.osm --host 172.16.49.31 --port 5433 --username gis_course --password gis_course123 --schema pgrouting --dbname gis_course"

system(command)

16.4 Finding the nodes

Osm2pgrouting sets up two spatial tables in the schema. One contains the edges (i.e. roads) of the network. This is named “ways”. The other table contains the vertices, i.e. the nodes that are at the ends of roads and at intersections. This is named ways_vertices_pgr. In order to navigate through the network we need to identify the nodes at the start and end of a route.

As we will usually have either an address or some coordinates on the map as input, rather than the nodes themsleves, we need a nearest neighbour query that takes coordinates as input and finds the nearest point on the network itself.. If we use the geocoding function in tmaptools to find coordinates then we can get the node for the centre of a town or other geocodable place by typing the name into a function.

require(tmaptools)
get_node<-function(place="Wareham"){
 g<- geocode_OSM(place)
 query<-sprintf("SELECT id, the_geom from pgrouting.ways_vertices_pgr
ORDER BY the_geom <-> 'SRID=4326;POINT(%s %s)'::geometry
LIMIT 1",g$coords[1],g$coords[2])
st_read(conn,query=query) 
}

Now we can find nodes just by typing a name, let’s plan a journey from Wareham to Bournemouth University.

get_node("Streche road, Wareham") -> start
get_node("Bournemouth University") -> end
mapview(start) + mapview(end) -> map
map

We will need a function to find the routes . There are several in pgrouting. Wrting sql to use these function is a little more complex than you might expect. The pgr_dijkstravia function returns only the identities of the edges. A join onto the table of edges is needed in order to return the roads.

The pgr_dijstraVia function allows any number of points to be passed to the function in order to form a trip and takes an array as input. The array can be formed pasting a comma between values in R.

shortest_route<-function(rt=c(33200,15389) ){
  query<-sprintf(" SELECT  seq, length_m, cost_s, the_geom
        FROM pgr_dijkstraVia('
                SELECT gid as id, source, target,
                        length as cost FROM pgrouting.ways',
                ARRAY[%s], false
        ) a INNER JOIN pgrouting.ways b ON (a.edge = b.gid) ORDER BY seq;",paste(rt,collapse=","))
  st_read(conn,query=query)
}

fastest_route<-function(rt=c(33200,15389) ){
  query<-sprintf(" SELECT  seq, length_m, cost_s, the_geom
        FROM pgr_dijkstraVia('
                SELECT gid as id, source, target,
                        cost_s as cost FROM pgrouting.ways',
                ARRAY[%s], false
        ) a INNER JOIN pgrouting.ways b ON (a.edge = b.gid) ORDER BY seq;",paste(rt,collapse=","))
  st_read(conn,query=query)
}

OK, the function does look rather complex but once written it is easy to adapt to other situations. It has been added to the giscourse package for convenience. Now, let’s plan a trip from Wareham to Bournemouth University.

fastest<-fastest_route(c(start$id,end$id) )
shortest<-shortest_route(c(start$id,end$id) )
map + mapview(fastest,color="red") + mapview(shortest,color="black")

Either one of these two suggested routes are reasonable options They conincide with the choices made by a habitual driver on the route.

How far was the journey?

sum(shortest$length_m)/1000
#> [1] 21.09539

That looks right and it coincides exactly with the value provided by Google maps.

How long does the journey take?

sum(fastest$cost_s)/60
#> [1] 20.28495

This is not right. The algorithm simply takes the maximum speeds allowed for the class of road. It does not take into account additional speed limits nor traffic. So the time is far too optimistic. Even without traffic the drive time is around 30 minutes and can increase to over an hour. So although pgrouting using OSM data provides reliable shortest distances, some further work to adjust the network would be needed is the intention was to provide more realistic drive times.

16.5 Round trip via Ferndown

Another example.

via<-get_node("Ferndown")
shortest<-shortest_route(c(start$id,via$id, end$id,start$id) )
fastest<-fastest_route(c(start$id,via$id, end$id,start$id) )
map + mapview(shortest,color="black") + mapview(fastest,color="red")

Again, additional work on the timings would be needed to improve the results in a satnav type application.

16.6 Driving distances

Pgrouting can reproduce more or less the same results that would be provided by Google maps using purely open source data and software. This has limited research value in the case of single trips, as Google provides accurate information without charge. However packages in R such as gmapsdistance now require a Google maps api, and if many requests are made the service has a charge.

A very useful feature of pgrouting in this context is the pgr_drivingDistance function. The function does not just provide driving distance. It finds all the distances in a network between one node and all others within a set distance. This is rather like v.net.iso in GRASS.

catchment_distance<-function(id= 15389,distance= 10) {
query<-sprintf("SELECT * FROM pgrouting.ways_vertices_pgr ver 
JOIN ( SELECT * FROM 
pgr_drivingDistance( 'SELECT gid id, source, target, 
length_m/1000 as cost, length_m/1000 as reverse_cost, 
the_geom FROM pgrouting.ways', %s, %s ) ) r 
on ver.id = r.node",id,distance)
  st_read(conn,query=query)
}


catchment_time<-function(id= 15389,time= 20) {
query<-sprintf("SELECT * FROM pgrouting.ways_vertices_pgr ver 
JOIN ( SELECT * FROM 
pgr_drivingDistance( 'SELECT gid id, source, target, 
cost_s/60 as cost, cost_s/60 as reverse_cost, 
the_geom FROM pgrouting.ways', %s, %s ) ) r 
on ver.id = r.node",id,time)
  st_read(conn,query=query)
}

bu_catch<-catchment_time()

The driving time may be rather unreliable, although as the road speeds are taken as the maximum speeds some calibration could be performed to improve the results.

The distances are as accurate as Google maps. Open street map is constantly updated so these results are useful as a means of defining catchment areas by distance along road networks, rather than euclidean distances.

The function returns all the nodes within the desired distance from a central point together with the aggregate travel cost from the source to each target.

bu_catch<-catchment_distance(distance=1)
mapview(bu_catch, zcol="agg_cost",legend=FALSE)

The data can be fed into an algorith to calculate concave hulls.

library(concaveman)
library(dplyr)
bu_catch<-catchment_distance(distance=12)
bu_catch %>% filter(agg_cost < 1)  %>% concaveman() ->km_1
bu_catch %>% filter(agg_cost < 3)  %>% concaveman() ->km_3
bu_catch %>% filter(agg_cost < 5)  %>% concaveman() ->km_5
bu_catch %>% filter(agg_cost < 10)  %>% concaveman() ->km_10

mapview(km_1,alpha.regions=0) + 
  mapview(km_3,alpha.regions=0) +
  mapview(km_5,alpha.regions=0) +
  mapview(km_10,alpha.regions=0)