## Finding the Shortest Path using the 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:

- Working out travel time from gdistance’s conductance values in R?
- Calculating accumulated walking pace from starting location outwards using R?
- 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))
}
```

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')
```

# (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)
```

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)
```

## 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)
```

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)
```

## 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')
```

## 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.