Shortest path with Tobler function

The following post was made trying to understand the functions withing the gdistance R package by Jacob van Etten.

Here is the documentation of the pacakge

During the process there will by answer the following questions found GIS StackExchange:

  1. Working out travel time from gdistance’s conductance values in R?
  2. Calculating accumulated walking pace from starting location outwards using R?
  3. Integrating use of surface distance when calculating cumulative cost-surface based on walking pace?

Basic functioning

library(dplyr)
library(gdistance)
library(raster)
library(sf)
library(RColorBrewer)
library(mapview)
library(htmlwidgets)
library(rgeos)
library(leaflet)

Usning a randomized raster with height values, we can see the cost of moving from an starting cell. Let’s create the matrix

r <- raster(nrows=6, ncols=7, 
            xmn=0, xmx=7, 
            ymn=0, ymx=6, 
            crs="+proj=utm +units=m")

r[] <- c(2, 2, 1, 1, 5, 5, 5,
         2, 2, 8, 8, 5, 2, 1,
         7, 1, 1, 8, 2, 2, 2,
         8, 7, 8, 8, 8, 8, 5,
         8, 8, 1, 1, 5, 3, 9,
         8, 1, 1, 2, 5, 3, 9)

Calculate CONDUCTANCE (transition) and plot with original raster:

# CONDUCTANCE: 1/mean: reciprocal to get permeability
tr <- transition(r, function(x) 1/mean(x), 8) 
tr <- geoCorrection(tr, scl=FALSE)

# plot raster and conductance
par(mfrow = c(1,2), mar=c(0,0,0,5))
plot(r)
text(r)
x <- raster(tr)
plot(x)
for(i in 1:ncell(x)){
 xycoords <- xyFromCell(x, cell = i)
 value <- extract(x, xycoords)
 text(xycoords, labels = round(value, 3))
}

R-plot

Now lets calculate the accumulated cost and the shortest path between cells

# Calculate accumulated cost
c1 <- c(1,6)
c2 <- c(6,2)
cost <- accCost(tr, c1)

# Shortest path
path1 <- shortestPath(tr, c1, c2, output="SpatialLines")
path2 <- shortestPath(tr, c2, c1, output="SpatialLines")

# Plots
par(mfrow=c(1,2))
# Accumulated cost
plot(cost)
text(cost)
# Shortest path
plot(cost)
plot(path1, add=TRUE, col = 'blue', lwd=5)
plot(path2, add=TRUE, col = 'red')

R-plot

(un) Real example

Now we will calculate the shortest routes between 3 points using a real digital elevation model an a real walking speed function that varies with topology (slope).

First, read raster dtm and generate random points to calculate routes:

# Read raster
r <- raster("../../data/20220720/DTM.tif")

# Select random coords within the raster
set.seed(1234)
A <- coordinates(r)[sample(dim(r)[[1]]*dim(r)[[2]])[[1]],]
B <- coordinates(r)[sample(dim(r)[[1]]*dim(r)[[2]])[[1]],]
C <- coordinates(r)[sample(dim(r)[[1]]*dim(r)[[2]])[[1]],]

# coords as sf
library(data.table)
cs <- data.table(ID=c("A", "B", "C"),
                 x=c(A[[1]], B[[1]], C[[1]]),
                 y=c(A[[2]], B[[2]], C[[2]]))

# st_as_sf() ######
ps <- st_as_sf(cs, coords = c("x", "y"), crs = 25829, agr = "constant")
ps <- ps %>% mutate(x=st_coordinates(ps)[[1]], y = st_coordinates(ps)[[2]])

# plot
plot(r)
points(as_Spatial(ps))
text(st_coordinates(ps), labels = ps$ID, pos = 4)

R-plot

Calculate cell differences between adjacent cells and apply the Toble’s function to correct speed:

# Calculate the height differences between cells (neet to calculate slope)
altDiff <- function(x){x[2] - x[1]}
suppressWarnings({ hd <- gdistance::transition(r, altDiff, 8, symm=FALSE) })
slope <- gdistance::geoCorrection(hd, scl=FALSE) # divide transition by distances between cells

# Calculate speed
adj <- raster::adjacent(x=r, cells=1:ncell(r), direction=8)
speed <- slope
speed[adj] <-  exp(-3.5 * abs(slope[adj] + 0.05)) # m/s

# Geocorrection divides speed by distance (i.e., cells' centers)
conduct <- gdistance::geoCorrection(speed, scl=FALSE) # conductance = speed/distance = 1/travel time

Plot of conductance:

plot(raster(conduct))
points(as_Spatial(ps))
text(st_coordinates(ps), labels = ps$ID, pos = 4)

R-plot

Calculate traveltime.

This is a solution to a QGIS Stack exchange question.

Extracted from “R Package gdistance: Distances and Routes on Geographical Grids” (Jacob van Etten)

What do these values mean? The function geoCorrection divides the values in the matrix between the distance between cell centres. So, with our last command we calculated this: conductance = speed/distance This looks a lot like a measure that we are more familiar with: traveltime = distance/speed In fact, the conductance values we have calculated are the reciprocal of travel time. 1/traveltime = speed/distance = conductance

So, just invert the last formula, so: traveltime = 1/conductance

tv <- raster(conduct)
values(tv) <- 1/values(raster(conduct))

See results of traveltime

plot(tv)
points(as_Spatial(ps))
text(st_coordinates(ps), labels = ps$ID, pos = 4)

R-plot

Calculate cost surface from A to the others

csA <- accCost(conduct, A)
csB <- accCost(conduct, B)
csC <- accCost(conduct, C)

#plot
par(mfrow=c(1,3))
plot(csA)
points(as_Spatial(ps))
text(st_coordinates(ps), labels = ps$ID, pos = 4)
plot(csB)
points(as_Spatial(ps))
text(st_coordinates(ps), labels = ps$ID, pos = 4)
plot(csC)
points(as_Spatial(ps))
text(st_coordinates(ps), labels = ps$ID, pos = 4)

R-plot

Can these results converted into isochrones ?

Here is an example of travel time from A to any other places by isochrones of an hour (3600).

library(rasterVis)
levels <- seq(min(values(csA)), max(values(csA)), 3600)
rasterVis::levelplot(csA, par.settings = BuRdTheme, margin=FALSE,
                     contour=TRUE, at = levels, lty=3,
                     labels = list(cex = 0.7), label.style = 'align')

R-plot

Calculate shortest path:

AB <- shortestPath(conduct, A, B, output="SpatialLines")
AC <- shortestPath(conduct, A, C, output="SpatialLines")
BA <- shortestPath(conduct, B, A, output="SpatialLines")
BC <- shortestPath(conduct, B, C, output="SpatialLines")
CA <- shortestPath(conduct, C, A, output="SpatialLines")
CB <- shortestPath(conduct, C, B, output="SpatialLines")

Check results with mapview:

Check if results make sense

We can check results by calculating the time expended

ACdist <- st_length(st_as_sf(AC))
kmh <- 5 # normal walking speed in km/h
msv <- kmh*(1/60)*(1/60)*1000 # normal walking speed in m/s
traveltimeAC <- as.numeric(ACdist) * msv
traveltimeAC
## [1] 4542.591
tvAC <- raster::extract(csA, ps[3,])[[1]]
tvAC
## [1] 4712.175

So, in a normal pace walking, going from A to C it will take 4543 seconds, while using the tobler function time will increase to 4712 seconds. As expected, topology will increase the travel time.