Data simulation

In order to simulate the sort of data that would be collected by manually placing points on a map I digitised 16 points in the Bournemouth area using the giscourse emap function. They were saved in the data base as a table called “test_points”

library(giscourse)
conn<-sconnect11()
pts<-st_read(conn,"test_points")
mapview(pts)

Converting points to pgrouting nodes

The open street map netowrk held in the pgrouting schema consists of nodes, with id numbers, and edges (lines) between them that represent the road network. The nodes are junctions and the end of cul de sacs. The nearest node can be found to any point whether or not it is actually on a street.

The following function takes the points and returns the nearest nodes on the network.

pnts_to_nodes<-function(pnts=pts)
{
cds<-st_coordinates(pnts)
f<-function(x){
query<-sprintf("SELECT id, the_geom from pgrouting.ways_vertices_pgr
ORDER BY the_geom <-> 'SRID=4326;POINT(%s %s)'::geometry
LIMIT 1",x[1],x[2])
st_read(conn,query=query) 
}
nds<-do.call("rbind",apply(cds,1,f))
nds}

nds<- pnts_to_nodes(pts)

Zooming in on the map below will show the way the points in red have been linked to nodes in green. The distances between the two are generally very small.

mapview(pts,color="red", col.regions="red", legend=FALSE) + mapview(nds, col.regions="darkgreen",legend=FALSE)

Finding the routes

The first node is the wntrance to Hengistbury head. So this is the end point.

end<-nds$id[1]

To obtain the routes apply a function to a list of all the other nodes.

list<-nds$id[-1]
f<-function(x)fastest_route(c(x,end))

routes<-lapply(list,f)
routes<-do.call("rbind",routes)
mapview(routes)

Road usage by section

In order to calculate how many cars use each section of road we need to convert the multistring routes into a set of lines representing the sections of road. This can be done by casting the geometry from multilinestring to linestring.

sections<-st_cast(routes,"LINESTRING")

Now to calculate the density of usage we group the unique line sections together and sumarise the number of times each one appears in the sf data object.

library(dplyr)
sections$id<-as.character(sections$the_geom)
sections %>% group_by(id) %>% summarise(n=n()) ->sections

Plotting

A simple mapview plot will shade the sections by the default colour. Hovering over the section gives a number showing usage.

mapview(sections,zcol="n")

Using the lwd parameter in mapview turns the lines invisible for some strange reason that I have not been able to fix.

The line widths can be shown using the tmap package. A little more work would be needed to add to this map in order to show other features and scale.

dor<-st_read(conn,query ="select * from Dorset")
library(tmap)
tm_shape(sections) +tm_lines(lwd="n",scale=20,col="red") +tm_scale_bar() +tm_compass() + tm_shape(dor) + tm_borders()

Example with more points

conn<-sconnect11()
set.seed(10)
library(concaveman)
pol<-concaveman(pts,100)
pnts2<-st_sample(pol,50)
nds2<-pnts_to_nodes(pnts2)
list<-nds2$id[-1]
f<-function(x)fastest_route(c(x,end))

routes<-lapply(list,f)
routes<-do.call("rbind",routes)
sections<-st_cast(routes,"LINESTRING")
sections$id<-as.character(sections$the_geom)
sections %>% group_by(id) %>% summarise(n=n()) ->sections

tm_shape(sections) +tm_lines(lwd="n",scale=20,col="red") +tm_scale_bar() +tm_compass() + tm_shape(dor) + tm_borders()