Support Our Open Source Project

cppRouting v4 : Algorithms for Routing and Solving the Traffic Assignment Problem

Package presentation

Installation

Stable version from CRAN

install.packages("cppRouting")

or from github

library(remotes)
remotes::install_github("vlarmet/cppRouting")

What are we talking about ?

Readme Data

library(cppRouting)
library(dplyr)
library(sf)
library(ggplot2)
library(concaveman)
library(ggmap)
library(tmap)
library(microbenchmark)
library(reshape2)
library(kableExtra)

#Reading french road data
roads  <-  read.csv("roads.csv",colClasses = c("character","character","numeric"))
#Shapefile data of communes (polygons)
com  <-  read_sf("com_simplified_geom.shp")
#Correspondance file between communes and nodes in the graph (nearest node to each commune centroid)
ndcom  <-  read.csv("node_commune.csv",colClasses = c("character","character","numeric"))
#General practitioners locations
med  <-  read.csv("doctor.csv",colClasses = c("character","numeric","character","numeric"))
#Import materinty ward locations
maternity  <-  read.csv("maternity.csv",colClasses = c("character","numeric"))
#Commuting data from national census
load("commuting.Rds")
#Import nodes coordinates (projected in EPSG : 2154)
coord  <-  read.csv("coordinates.csv",colClasses = c("character","numeric","numeric"))

Head of road network data

head(roads)

##   from to    weight
## 1    0  1 2.1200000
## 2    1  0 2.1200000
## 3    1  2 0.8600000
## 4    2  1 0.8600000
## 5    3  4 0.1745455
## 6    4  3 0.1745455

Head of coordinates data

head(coord)

##   ID        X       Y
## 1  0 894833.5 6265744
## 2  1 894051.2 6262806
## 3  2 893427.9 6261751
## 4  3 340196.8 6276241
## 5  4 340108.9 6275932
## 6  5 433102.7 6715832

Set the number of threads used by cppRouting

RcppParallel::setThreadOptions(numThreads = 1)

Instantiate the graph

#Instantiate a graph with coordinates
graph  <-  makegraph(roads, directed = T, coords = coord)

Main functions

Network extraction

# path to ROUTE500 LINESTRING shapefile
shp <- ".../RESEAU_ROUTIER/TRONCON_ROUTE.shp"

# extract network, nodes coordinates, auxiliary edges weights (length) in csv files
cppRouting::cpp_extract(
  shp, 
  output = "csv", 
  out_prefix = "route500", 
  coords = TRUE, 
  keep_length = TRUE,
  geom_params = list(
    speed_field = "vitesse",
    length_field = "LONGUEUR",
    direction_field = "SENS",
    direction_values = list(forward = "Sens direct", backward = "Sens inverse", twoway = "Double sens"),
    keep_fields = "LONGUEUR"
  )
)

Routing

Algorithms

#Generate 2000 random origin and destination nodes
origin  <-  sample(graph$dict$ref,  2000)
destination  <-  sample(graph$dict$ref,  2000)
microbenchmark(dijkstra=pair_dijkstra  <-  get_distance_pair(graph, origin, destination, algorithm = "Dijkstra"), 
               bidir=pair_bidijkstra  <-  get_distance_pair(graph, origin, destination, algorithm = "bi"), 
               astar=pair_astar  <-  get_distance_pair(graph, origin, destination, algorithm = "A*", constant = 110/0.06), 
               nba=pair_nba  <-  get_distance_pair(graph, origin, destination, algorithm = "NBA", constant = 110/0.06), 
               times=1)

## Unit: seconds
##      expr       min        lq      mean    median        uq       max neval
##  dijkstra 122.87958 122.87958 122.87958 122.87958 122.87958 122.87958     1
##     bidir  83.18095  83.18095  83.18095  83.18095  83.18095  83.18095     1
##     astar  74.30699  74.30699  74.30699  74.30699  74.30699  74.30699     1
##       nba  34.97340  34.97340  34.97340  34.97340  34.97340  34.97340     1

Output

head(cbind(pair_dijkstra,pair_bidijkstra,pair_astar,pair_nba))

##      pair_dijkstra pair_bidijkstra pair_astar  pair_nba
## [1,]     191.74489       191.74489  191.74489 191.74489
## [2,]     519.59685       519.59685  519.59685 519.59685
## [3,]     308.14500       308.14500  308.14500 308.14500
## [4,]      77.99867        77.99867   77.99867  77.99867
## [5,]     194.39985       194.39985  194.39985 194.39985
## [6,]     314.71795       314.71795  314.71795 314.71795

Compute isochrones

#Node nearest to the Dijon commune centroid (from the node-to-commune correspondence file)
dijon_id <- as.character(ndcom$id_noeud[ndcom$com == "21231"])
#Compute isochrones
iso  <-  get_isochrone(graph, from = dijon_id, lim = c(15, 25, 45, 60, 90, 120))
#Convert nodes to concave polygons with concaveman package
poly  <-  lapply(iso[[1]], function(x){
  x  <-  data.frame(noeuds=x, stringsAsFactors = F)
  x  <-  left_join(x, coord, by=c("noeuds"="ID"))
  return(concaveman(st_as_sf(x, coords=c("X", "Y"), crs=2154)))
})

poly  <-  do.call(rbind, poly)
poly$time  <-  as.factor(names(iso[[1]]))
#Multipolygon
poly2  <-  st_cast(poly, "MULTIPOLYGON")
poly2$time  <-  reorder(poly2$time, c(120, 90, 60, 45, 25, 15))
#Draw the largest isochrone first, so the smaller ones are not hidden
poly2  <-  poly2[order(as.integer(poly2$time)), ]
#Reproject for plotting
poly2  <-  st_transform(poly2, "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
#Plot the map over the commune polygons, reprojected in longlat like the isochrones
com_ll  <-  st_transform(com, "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
p  <-  ggplot() +
  geom_sf(data = com_ll, fill = "white", color = "grey85", linewidth = 0.1) +
  geom_sf(data=poly2, aes(fill=time), alpha=.8, inherit.aes = FALSE)+
  scale_fill_brewer(palette = "YlOrRd")+
  labs(fill="Minutes")+
  ggtitle("Isochrones around Dijon")+
  coord_sf(xlim = c(2.5, 8), ylim = c(45.126, 49)) +
  theme(axis.text.x = element_blank(), 
        axis.text.y = element_blank(), 
        axis.ticks = element_blank(), 
        axis.title.y=element_blank(), axis.title.x=element_blank())
p

Compute possible detours within a fixed additional cost

#Node nearest to the Lyon city centroid (Lyon is not in the node-to-commune correspondence file)
lyon_pt <- st_transform(st_sfc(st_point(c(4.8706, 45.777)), crs = 4326), 2154)
lyon_xy <- st_coordinates(lyon_pt)[1, ]
lyon_id <- coord$ID[which.min((coord$X - lyon_xy[1])^2 + (coord$Y - lyon_xy[2])^2)]
#Compute shortest path
trajet <- get_path_pair(graph,from=dijon_id,to=lyon_id)

#Compute shortest path
distance <- get_distance_pair(graph,from=dijon_id,to=lyon_id)

#Compute detour time of 25 and 45 minutes
det25 <- get_detour(graph,from=dijon_id,to=lyon_id,extra=25)
det45 <- get_detour(graph,from=dijon_id,to=lyon_id,extra=45)

#Create sf object of nodes
pts <- st_as_sf(coord,coords=c("X","Y"),crs=2154)
pts <- st_transform(pts,crs=4326)
pts$time <- ifelse(pts$ID %in% unlist(det45),"45","0")
pts$time <- ifelse(pts$ID %in% unlist(det25),"25",pts$time)
pts$time <- ifelse(pts$ID %in% unlist(trajet),"Shortest Path",pts$time)
pts$time <- factor(pts$time,levels = c("25","45","Shortest Path","0"))

#Plot
#Plot over the commune polygons, reprojected in longlat like the points
com_ll <- st_transform(com, "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")

p <- ggplot() +
  geom_sf(data = com_ll, fill = "white", color = "grey85", linewidth = 0.1) +
  geom_sf(data=pts[pts$time %in% c("25","45"),],aes(color=time),inherit.aes = FALSE)+
  geom_sf(data=pts[pts$time =="Shortest Path",],aes(color=time),inherit.aes = FALSE)+
  scale_color_brewer(palette = "YlOrRd") +
  ggtitle(paste0("Detours around Dijon-lyon path - ",round(distance,digits = 2)," minutes"))+
  labs(color="Minutes") +
  coord_sf(xlim = c(3.2, 6.5), ylim = c(45.126, 47.8)) +
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        axis.title.y=element_blank(),axis.title.x=element_blank())
p

Contraction hierarchies (CH)

#Contraction of input graph using CH
graph3 <- cpp_contract(graph, silent=TRUE)

#Calculate distances on the contracted graph
system.time(
  pair_ch <- get_distance_pair(graph3, origin, destination)
)

## utilisateur     système      écoulé 
##       0.701       0.014       0.715

Compare outputs

summary(pair_ch-pair_dijkstra)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##       0       0       0       0       0       0      14

Customizable contraction hierarchies (CCH)

#Contraction of input graph using CCH
cch_graph <- cpp_contract(graph, customizable = TRUE, silent=TRUE)

# Assign original weights
graph4 <- cpp_customize(cch_graph)

#Calculate distances on the contracted graph
system.time(
  pair_cch <- get_distance_pair(graph4, origin, destination)
)

## utilisateur     système      écoulé 
##       2.610       0.000       2.609

Compare outputs

summary(pair_cch-pair_dijkstra)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##       0       0       0       0       0       0      14

Performance comparison

Distance by pair

Distance matrix

Work with multiple weights network

# The weight to be minimized is set in 'df'; the additional weights are set in 'aux'.
# 'aux' is a data.frame with one column per auxiliary weight. Here we set two:
#  - n_edges : a constant of 1, so summing it counts the edges in each shortest path
#  - length  : the length of each edge (km), read from roads_aux.csv (aligned with roads)
auxdf <- read.csv("roads_aux.csv", header = TRUE, stringsAsFactors = FALSE)
aux_graph <- makegraph(df = roads, directed = TRUE, coords = coord,
                 aux = data.frame(n_edges = 1, length = auxdf$length))
# 2000 random origin-destination pairs
set.seed(1)
from <- sample(aux_graph$dict$ref, 2000)
to <- sample(aux_graph$dict$ref, 2000)
# Travel time, and both auxiliary weights aggregated along each shortest path
# (aggregate_aux = TRUE returns one data.frame with the aggregated aux columns)
time_pair <- get_distance_pair(aux_graph, from = from, to = to)
aux_pair <- get_distance_pair(aux_graph, from = from, to = to, aggregate_aux = TRUE)
res <- data.frame(travel_time = time_pair, aux_pair)
res$length[is.na(res$length)] <- 0
res$n_edges[is.na(res$n_edges)] <- 0

# Scatter plot 1: travel time vs number of edges
p1 <- ggplot(data = res, aes(x = n_edges, y = travel_time)) +
  geom_point() +
  labs(x = "Number of edges", y = "Travel time (minutes)")
p1

# Scatter plot 2: travel time vs path length
p2 <- ggplot(data = res, aes(x = length, y = travel_time)) +
  geom_point() +
  labs(x = "Path length (km)", y = "Travel time (minutes)")
p2

Network simplification

library(igraph)
#Create directed graph
edges<-data.frame(from=c("a","b","c","d",
                         "d","e","e","e",
                         "f","f","g","h","h","h",
                         "i","j","k","k","k",
                         "l","l","l","m","m","m",
                         "n","n","o","p","q","r"),
                  to=c("b","c","d","e","k","f","d",
                       "h","g","e","f","e","i","k",
                       "j","i","h","d","l","k",
                       "m","n","n","o","l","l","m","m",
                       "r","p","q"),
                  dist=rep(1,31))

#Plotting with igraph
par(mfrow=c(1,2),mar=c(3,0,3,0))

igr1<-graph_from_data_frame(edges)
set.seed(2)
plot.igraph(igr1,edge.arrow.size=.3,main="Original graph")
box(col="black")

#Instantiate cppRouting graph, then simplify without iterations
graph_ex<-makegraph(edges,directed = TRUE)
simp<-cpp_simplify(graph_ex,rm_loop = FALSE)
#Convert graph to df
edges2<-to_df(simp)

#Plotting simplified graph
igr2<-graph_from_data_frame(edges2)
set.seed(20)
plot(igr2,edge.arrow.size=.3,edge.label=E(igr2)$dist,main="One iteration - keeping loop")
box(col="black")

par(mfrow=c(1,2),mar=c(3,0,3,0))
#Simplify with iterations
simp2<-cpp_simplify(graph_ex,rm_loop = FALSE,iterate = TRUE)
edges3<-to_df(simp2)
igr3<-graph_from_data_frame(edges3)
set.seed(2)
plot(igr3,edge.arrow.size=.3,edge.label=E(igr3)$dist,main="Second iteration - keeping loop")
box(col="black")

#The same but removing loops
simp3<-cpp_simplify(graph_ex,rm_loop = TRUE,iterate = TRUE)
edges4<-to_df(simp3)

igr4<-graph_from_data_frame(edges4)
set.seed(2)
plot(igr4,edge.arrow.size=.3,edge.label=E(igr4)$dist,main="Second iteration - removing loop")
box(col="black")

French road network simplification

#Simplify original graph by keeping nodes of interest
graph2<-cpp_simplify(graph,
                     iterate = TRUE,
                     keep = unique(c(origin,destination))) #we don't want to remove origin and destination nodes

#Running NBA*
system.time(
  pair_nba_2<-get_distance_pair(graph2,origin,destination,algorithm = "NBA",constant = 110/0.06)
)

## utilisateur     système      écoulé 
##      25.151       0.007      25.159

Compare outputs

summary(pair_nba-pair_nba_2)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##       0       0       0       0       0       0      14

Running time

Traffic assignment

net <- read.csv("siouxfalls_net.csv")
head(net)

##   Init.node Term.node  Capacity Length Free.Flow.Time    B Power Speed.limit
## 1         1         2 25900.201      6              6 0.15     4           0
## 2         1         3 23403.473      4              4 0.15     4           0
## 3         2         1 25900.201      6              6 0.15     4           0
## 4         2         6  4958.181      5              5 0.15     4           0
## 5         3         1 23403.473      4              4 0.15     4           0
## 6         3         4 17110.524      4              4 0.15     4           0
##   Toll Type
## 1    0    1
## 2    0    1
## 3    0    1
## 4    0    1
## 5    0    1
## 6    0    1

trips <- read.csv("siouxfalls_trips.csv")
head(trips)

##   from to demand
## 1    1  2    100
## 2    1  3    100
## 3    1  4    500
## 4    1  5    200
## 5    1  6    300
## 6    1  7    500

All-or-Nothing (AON)

sgr <- makegraph(df = net[,c("Init.node", "Term.node", "Free.Flow.Time")], directed = TRUE) 

# Compute AON assignment using OD trips
aon <- get_aon(Graph = sgr, from = trips$from, to = trips$to, demand = trips$demand)
head(aon)

##   from to cost  flow
## 1    1  2    6  3800
## 2    1  3    4  6000
## 3    2  1    6  3800
## 4    2  6    5  6600
## 5    3  1    4  6000
## 6    3  4    4 10200

g <- graph_from_data_frame(aon[,c("from", "to", "flow")], directed = TRUE)
c_scale <- colorRamp(c('#60AF66', 'yellow', 'red'))
maxflow <- max(c(aon$flow))

E(g)$color = apply(c_scale(aon$flow/maxflow), 1, function(x) rgb(x[1]/255,x[2]/255,x[3]/255) )
set.seed(5)
plot(g, edge.width = net$Capacity/2000, edge.arrow.size = 0, main = "AON assignment")
box(col = "black")

Choosing the best routing algorithm

User Equilibrium (UE)

Link-based algorithms

sgr <- makegraph(df = net[,c("Init.node", "Term.node", "Free.Flow.Time")], 
                 directed = TRUE,
                 capacity = net$Capacity,
                 alpha = net$B,
                 beta = net$Power)


traffic <- assign_traffic(Graph = sgr,  from = trips$from, to = trips$to, demand = trips$demand, 
                          max_gap = 1e-6, algorithm = "bfw", verbose = FALSE)
  
traffic$gap

## [1] 0.0000009230227

head(traffic$data)

##   from to ftt     cost      flow  capacity alpha beta
## 1    1  2   6 6.000779  4442.890 25900.201  0.15    4
## 2    1  3   4 4.008658  8111.525 23403.473  0.15    4
## 3    2  1   6 6.000829  4511.531 25900.201  0.15    4
## 4    2  6   5 6.545270  5940.297  4958.181  0.15    4
## 5    3  1   4 4.008369  8042.884 23403.473  0.15    4
## 6    3  4   4 4.262114 13910.673 17110.524  0.15    4

par(mfrow=c(1,2),mar=c(3,0,3,0))
g <- graph_from_data_frame(aon[,c("from", "to", "flow")], directed = TRUE)
c_scale <- colorRamp(c('#60AF66', 'yellow', 'red'))
maxflow <- max(c(aon$flow, traffic$data$flow))

#Applying the color scale to edge weights.
#rgb method is to convert colors to a character vector.
E(g)$color = apply(c_scale(aon$flow/maxflow), 1, function(x) rgb(x[1]/255,x[2]/255,x[3]/255) )
set.seed(5)
plot(g, edge.width = traffic$data$capacity/2000, edge.arrow.size = 0, main = "AON assignment")
box(col = "black")


g <- graph_from_data_frame(traffic$data[,c("from", "to", "flow")], directed = TRUE)
E(g)$color = apply(c_scale(traffic$data$flow/maxflow), 1, function(x) rgb(x[1]/255,x[2]/255,x[3]/255) )
set.seed(5)
plot(g, edge.width = traffic$data$capacity/2000, edge.arrow.size = 0, main = "Flows at equilibrium")
box(col = "black")

Bush-based algorithms

traffic <- assign_traffic(Graph = sgr,  from = trips$from, to = trips$to, demand = trips$demand, 
                          max_gap = 1e-6, algorithm = "dial", verbose = TRUE)

## Bushes initialization...
## Iterating...
## iteration 1 : 0.05019
## iteration 2 : 0.0152742
## iteration 3 : 0.00363109
## iteration 4 : 0.000592523
## iteration 5 : 0.000106401
## iteration 6 : 3.00525e-05
## iteration 7 : 3.5758e-05
## iteration 8 : 1.21886e-05
## iteration 9 : 6.67853e-06
## iteration 10 : 1.659e-06
## iteration 11 : 2.34866e-07

traffic$gap

## [1] 0.000000234866

Performance comparison

net <- read.csv("chicagoregional_net.csv")
head(net)

##   from    to capacity length ftime    B power speed toll type
## 1    1 10293   100000   0.45     0 0.15     4    25    0    3
## 2    2 10294   100000   0.50     0 0.15     4    25    0    3
## 3    3 10295   100000   0.39     0 0.15     4    25    0    3
## 4    4 10296   100000   0.58     0 0.15     4    25    0    3
## 5    5 10297   100000   0.50     0 0.15     4    25    0    3
## 6    6 10298   100000   0.58     0 0.15     4    25    0    3

trips <- read.csv("chicagoregional_trips.csv")
head(trips)

##   from to demand
## 1    1  1  28.94
## 2    1  2  16.58
## 3    1  3   2.01
## 4    1  4  13.99
## 5    1  5   4.90
## 6    1  6   5.10

# Construct graph with link attributes
sgr <- makegraph(df = net[,c("from", "to", "ftime")], directed = TRUE,
                 capacity = net$capacity, alpha = net$B, beta = net$power)

# benchmark using all cores. We don't test NBA* because we don't have node coordinates
RcppParallel::setThreadOptions(numThreads = parallel::detectCores())

cch <- cpp_contract(sgr, customizable = TRUE, silent = TRUE)

microbenchmark(
  pairwise = pw <- get_aon(Graph = sgr, from = trips$from, to = trips$to, demand = trips$demand, algorithm = "bi"),
  matrixlike = ml <- get_aon(Graph = sgr, from = trips$from, to = trips$to, demand = trips$demand, algorithm = "d"),
  cch = {
    cch_metric <- cpp_customize(cch)
    caon <- get_aon(Graph = cch_metric, from = trips$from, to = trips$to, demand = trips$demand)
  },
  times = 1
)

## Unit: seconds
##        expr       min        lq      mean    median        uq       max neval
##    pairwise 55.368205 55.368205 55.368205 55.368205 55.368205 55.368205     1
##  matrixlike  1.092179  1.092179  1.092179  1.092179  1.092179  1.092179     1
##         cch  1.909323  1.909323  1.909323  1.909323  1.909323  1.909323     1

sgr <- makegraph(df = net[,c("from", "to", "ftime")], directed = TRUE,
                 capacity = net$capacity, alpha = net$B, beta = net$power)

# we set 5 minutes for each algorithm (about 600 iterations for link-based, and 20 for algorithm-b)
## capture console output for plotting
outs <- list()
times <- list()
for (i in c("msa", "fw", "cfw", "bfw", "dial")){
  max_it <- 600
  if (i == "dial"){
    max_it <- 20
  }
  
  stdout <- vector(mode = "character")
  con    <- textConnection('stdout', 'wr', local = TRUE)
  sink(con)
  st <- system.time(
  traffic <- assign_traffic(Graph = sgr,  from = trips$from, to = trips$to, demand = trips$demand, 
                          max_gap = 1e-6, algorithm = i, aon_method = "d", verbose = TRUE, max_it = max_it)
  )
  sink()
  outs[[i]] <- stdout
  times[[i]] <- st
  rm(stdout)
}

# Plotting gap vs time

dfp <- list()
for (i in c("msa", "fw", "cfw", "bfw", "dial")){
  out <- outs[[i]]
  out <- out[grepl("iteration", out)]
  out <- sapply(strsplit(out, " : "), function(x) as.numeric(x[2]))
  
  df <- data.frame(gap = out, time = 1:length(out)  * (times[[i]][1]/length(out)), algo = i)
  dfp[[i]] <- df
}
dfp <- do.call(rbind, dfp)

options(scipen=999)
p <- ggplot(dfp, aes(x = time, y = gap, colour = algo))+
  geom_line(size = 1)+
  scale_y_log10()+
  theme_bw()+
  labs(x = "Time(second)", y = "Relative gap (log10 scale)", colour = "Algorithm")+
  scale_color_hue(labels = c("Bi-Conjugate FW", "Conjugate FW", "Algorithm B", "Frank-Wole", "MSA"))
  
p

Algorithm compatibility

Parallel implementation

Applications

Application 1 : Calculate Two Step Floating Catchment Areas (2SFCA) of general practitioners in France

#Isochrone around doctor locations with time limit of 15 minutes
iso<-get_isochrone(graph,
                   from = ndcom[ndcom$com %in% med$CODGEO,"id_noeud"],
                   lim = 15,
                   keep=ndcom$id_noeud, #return only commune nodes
                   long = TRUE) #data.frame output
#Joining and summing population located in each isochrone
iso<-left_join(iso,ndcom[,c("id_noeud","POPULATION")],by=c("node"="id_noeud"))
df<-iso %>% group_by(origin) %>%
  summarise(pop=sum(POPULATION))
#Joining number of doctors 
df<-left_join(df,med[,c("id_noeud","NB_D201")],by=c("origin"="id_noeud"))
#Calculate ratios
df$ratio<-df$NB_D201/df$pop

#Isochrone around each commune with time limit of 15 minutes (few seconds to compute)
iso2<-get_isochrone(graph,
                    from=ndcom$id_noeud,
                    lim = 15,
                    keep=ndcom$id_noeud,
                    long=TRUE)
#Joining and summing ratios calculated in first step
df2<-left_join(iso2,df[,c("origin","ratio")],by=c("node"="origin"))
df2<-df2 %>% group_by(origin) %>%
  summarise(sfca=sum(ratio,na.rm=T))

#Joining commune IDs to nodes
df2<-left_join(df2,ndcom[,c("id_noeud","com")],by=c("origin"="id_noeud"))
#Joining 2SFCA to shapefile
com<-left_join(com,df2[,c("com","sfca")],by=c("INSEE_COM"="com"))
#Plot for one region
p<-tm_shape(com[com$NOM_REG=="BOURGOGNE-FRANCHE-COMTE",]) + 
  tm_fill(col = "sfca",style="cont", lwd=0)+
  tm_layout(main.title="2SFCA applied to general practitioners",legend.outside=TRUE)
p

Application 2 : Calculate the minimum travel time to the closest maternity ward in France

#Distance matrix on contracted graph
dists<-get_distance_matrix(graph3,
                           from=ndcom$id_noeud,
                           to=ndcom$id_noeud[ndcom$com %in% maternity$CODGEO],
                           algorithm = "phast")#because of the rectangular shape of the matrix
#We extract each minimum travel time for all the communes
dists2<-data.frame(node=ndcom$id_noeud,mindist=apply(dists,1,min,na.rm=T))
#Joining commune IDs to nodes
dists2<-left_join(dists2,ndcom[,c("id_noeud","com")],by=c("node"="id_noeud"))
#Joining minimum travel time to the shapefile
com<-left_join(com,dists2[,c("com","mindist")],by=c("INSEE_COM"="com"))

p<-tm_shape(com[com$NOM_REG=="BOURGOGNE-FRANCHE-COMTE",]) + 
  tm_fill(col = "mindist",style="cont",palette="Greens",title="Minutes", lwd=0)+
  tm_layout(main.title="Travel time to the closest maternity ward",legend.outside=T)
p

Application 3 : Calculate average commuting time to go to job in France

#Merge node to communes
ndcom$id_noeud<-as.character(ndcom$id_noeud)
cmt$node1<-ndcom$id_noeud[match(cmt$CODGEO,ndcom$com)]
cmt$node2<-ndcom$id_noeud[match(cmt$DCLT,ndcom$com)]
cmt<-cmt[!is.na(cmt$node1) & !is.na(cmt$node2),]

#Calculate distance for each pair using contracted graph
dist<-get_distance_pair(graph3,from=cmt$node1,to=cmt$node2)

#Mean weighted by the flow of each pair
traveltime<-cbind(cmt,dist)
traveltime<-traveltime %>% group_by(CODGEO) %>% summarise(time=weighted.mean(dist,flux,na.rm=T))

#Merge to shapefile
com_time <- left_join(com, traveltime, by = c("INSEE_COM"="CODGEO"))
p <- tm_shape(com_time[com_time$NOM_REG=="BOURGOGNE-FRANCHE-COMTE", ]) +
  tm_fill(col = "time",style="quantile",n=5,title="Minutes",palette="Reds", lwd=0)+
  tm_layout(main.title="Average travel time",legend.outside=TRUE)
p

Application 4 : Calculate the flow of people crossing each municipality in the context of commuting in Bourgogne-Franche-Comte region

#Convert nodes coordinates to spatial points
pts<-st_as_sf(coord,coords = c("X","Y"),crs=2154)
st_crs(com)<-st_crs(pts)
#One row per commune: the app 1/2 joins may have duplicated communes sharing
#the same nearest node (node_commune.csv is not unique by node)
com<-com[!duplicated(com$INSEE_COM),]
#Spatial join commune ID to each node
spjoin<-st_join(pts,com[,"INSEE_COM"])
st_geometry(spjoin)<-NULL
spjoin$ID<-as.character(spjoin$ID)
#Vector of commune IDs of Bourgogne-Franche-Comte region
comBourg<-com$INSEE_COM[com$NOM_REG=="BOURGOGNE-FRANCHE-COMTE"]
#Calculate shortest paths for all commuters on the contracted graph
shortPath<-get_path_pair(graph3,
                         from=cmt$node1,
                         to=cmt$node2,
                         keep = spjoin$ID[spjoin$INSEE_COM %in% comBourg],#return only the nodes in Bourgogne Franche-Comte because of memory usage
                         long = TRUE)#return a long data.frame instead of a list

#Joining commune ID to each traveled node
shortPath<-left_join(shortPath,spjoin,by=c("node"="ID"))

#Remove duplicated 
shortPath<-shortPath[!duplicated(shortPath[,c("from","to","INSEE_COM")]),]

#Joining flow per each commuter
shortPath<-left_join(shortPath,
                     cmt[,c("node1","node2","flux")],
                     by=c("from"="node1","to"="node2"))

#Aggregating flows by each commune
shortPath<-shortPath %>% group_by(INSEE_COM) %>% summarise(flow=sum(flux))

#Merge to shapefile
com<-left_join(com,shortPath,by="INSEE_COM")

p<-tm_shape(com[com$NOM_REG=="BOURGOGNE-FRANCHE-COMTE",]) + 
  tm_fill(col = "flow",style="quantile",n=10,palette="Blues", lwd=0)+
  tm_layout(main.title="Commuters flow",legend.outside=TRUE)
p

Application 5 : Estimate the commuter traffic at user equilibrium in Bourgogne-Franche-Comte

#Nodes located in the region (node ids are 0-based, aligned with the road file)
nodecom  <- spjoin$INSEE_COM
bfc_nodes <- which(nodecom %in% comBourg) - 1L
#Sub-network: keep only arcs whose two endpoints are in the region
keep_e <- (as.integer(roads$from) %in% bfc_nodes) & (as.integer(roads$to) %in% bfc_nodes)
sub_road <- roads[keep_e, ]
old_ids <- sort(unique(c(as.integer(sub_road$from), as.integer(sub_road$to))))
#Renumber the nodes from 0
idmap <- setNames(as.character(seq_along(old_ids) - 1L), as.character(old_ids))
sub_road$from <- idmap[as.character(sub_road$from)]
sub_road$to   <- idmap[as.character(sub_road$to)]
sub_coord <- coord[match(old_ids, as.integer(coord$ID)), ]
sub_coord$ID <- as.character(seq_along(old_ids) - 1L)
#Keep only origin-destination pairs whose two endpoints are in the sub-network
od <- cmt[as.character(cmt$node1) %in% old_ids & as.character(cmt$node2) %in% old_ids, ]
od$node1 <- idmap[as.character(od$node1)]
od$node2 <- idmap[as.character(od$node2)]
cat(sprintf("Sub-network: %d nodes, %d arcs - %d origin-destination pairsn",
            length(old_ids), nrow(sub_road), nrow(od)))

## Sub-network: 64310 nodes, 179340 arcs - 47628 origin-destination pairs

#Traffic assignment needs the BPR parameters, defined at graph construction
RcppParallel::setThreadOptions(numThreads = parallel::detectCores() - 1) #the AON step is multithreaded
g_ue <- makegraph(sub_road, coords = sub_coord, directed = TRUE,
                  capacity = 1800, alpha = 0.15, beta = 4)
cch_ue <- cpp_contract(g_ue, customizable = TRUE) #topology prepared once

## Computing elimination order...
## Eliminating nodes...

ue <- assign_traffic(g_ue, from = od$node1, to = od$node2, demand = od$flux,
                     algorithm = "bfw", aon_method = "cch", cch = cch_ue,
                     max_gap = 0.001, verbose = FALSE)
cat(sprintf("Gap achieved: %.5f (%d iterations)n", ue$gap, ue$it))

## Gap achieved: 0.00100 (635 iterations)

#Coordinates of the two endpoints of each arc
cof <- sub_coord; colnames(cof) <- c("ID_f","X_f","Y_f")
cot <- sub_coord; colnames(cot) <- c("ID_t","X_t","Y_t")
ue_df <- left_join(ue$data, cof, by = c("from" = "ID_f"))
ue_df <- left_join(ue_df, cot, by = c("to" = "ID_t"))
ue_df <- ue_df[!is.na(ue_df$X_f) & !is.na(ue_df$X_t) & ue_df$flow > 0, ]
ue_df$geom <- mapply(function(a,b,c,dd) st_linestring(rbind(c(a,b), c(c,dd))),
                     ue_df$X_f, ue_df$Y_f, ue_df$X_t, ue_df$Y_t, SIMPLIFY = FALSE)
ue_df <- st_sf(ue_df, crs = 2154)
com_bfc <- com[com$NOM_REG=="BOURGOGNE-FRANCHE-COMTE", ]
bbox <- st_bbox(com_bfc$geometry)
p <- ggplot() +
  geom_sf(data = com_bfc, fill = "white", color = "grey85", linewidth = 0.2) +
  geom_sf(data = ue_df, aes(size = flow, color = flow), alpha = 0.9) +
  scale_size_continuous(range = c(0.2, 2.5), trans = "sqrt", name = "Commuter flow (people)") +
  scale_color_gradient(low="yellow", high = "red", trans = "log10", name = NULL, guide = "none") +
  coord_sf(xlim = bbox[c("xmin","xmax")], ylim = bbox[c("ymin","ymax")]) +
  theme_void()
p

Citation

citation("cppRouting")