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 consisting simply of lines is imported into a GIS such as GRASS and converted directly into a network some features of a real road network will be overlooked. A straightforward conversion of interescting lines into 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 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.
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
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)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)
}So now we can find nodes 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
mapNow we need a function to find the routes. The sql to use the function is a little more complex than you might expect. The pgr_dijkstra function returns only the identities of the edges, so a join on the table 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 looks rather complex but it is easy to adapt to other situations. 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 the two suggested routes would be reasonable and 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.09539That looks right and coincides with the value provided by Google maps .
How long does it take?
sum(fastest$cost_s)/60## [1] 20.28495This 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.
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.
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) The code shown here could be adapted to provide some realistic adjustments to the travel times to suit the Dorset traffic.
https://geospatialwandering.wordpress.com/2015/09/09/pg-routing-shortest-path-algorithms/#more-220
For example ..
query<-"alter table pgrouting.ways
add column real_speed double precision,
add column time_car double precision,
add column time_walking double precision,
add column len_miles double precision,
add column time_cycling double precision;
update pgrouting.ways set len_miles = length*0.6214;
"
dbSendQuery(conn,query)query<-"update pgrouting.ways set real_speed = case
when class_id = 101 then 50
when class_id = 102 then 50
when class_id = 103 then 50
when class_id = 106 then 40
when class_id = 107 then 40
when class_id = 108 or class_id = 109 then 30
when class_id = 110 or class_id = 111 then 20
else 10 end;
update pgrouting.ways set time_car = ((60/real_speed::double precision)*len_miles);
update pgrouting.ways set time_walking = ((60/3::double precision)*len_miles);
update pgrouting.ways set time_cycling = ((60/9.32::double precision)*len_miles);"
dbSendQuery(conn,query)nodes<-function(x="Arne,Dorset",y=0,dist=0){
site<-mksite(x,y,dist)
query<- "select p.*
from pgrouting.ways_vertices_pgr p, site t
where st_intersects(p.the_geom,st_transform(t.geometry,4326))"
nodes<-st_read(conn,query=query)
dbSendQuery(conn, "drop table site")
nodes
}
ways<-function(x="Arne,Dorset",y=0,dist=0){
site<-mksite(x,y,dist)
query<- "select class_id,the_geom
from pgrouting.ways p, site t
where st_intersects(p.the_geom,st_transform(t.geometry,4326))"
ways<-st_read(conn,query=query)
dbSendQuery(conn, "drop table site")
ways
}