 Estudios de paisaje, cartografía (R y QGIS), diseño ambiental 3D (Blender) y análisis de datos geográficos en R

##### Contact Info
Ferrol, Santiago de Compostela, Galicia, Spain, The World

+00(34) not yet...

## Finding the Shortest Path using the Tobler Function

20 Jul 2022

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:

## 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)
`````` # (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)[]*dim(r)[])[],]
B <- coordinates(r)[sample(dim(r)[]*dim(r)[])[],]
C <- coordinates(r)[sample(dim(r)[]*dim(r)[])[],]

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

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

# 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 - x}
suppressWarnings({ hd <- gdistance::transition(r, altDiff, 8, symm=FALSE) })
slope <- gdistance::geoCorrection(hd, scl=FALSE) # divide transition by distances between cells

# Calculate speed
speed <- slope

# 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
``````
``````##  4542.591
``````
``````tvAC <- raster::extract(csA, ps[3,])[]
tvAC
``````
``````##  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.