Calculating accumulated walking pace from starting location outwards using R?

Question

Goal

In R, I would like to calculate an accumulated cost surface around a starting location, with the cost representing the time (say, hours) to walk from the location outwards in all directions (say, 8 directions in raster's terms). The walking time (pace) is devised via the Tobler's hiking function (https://en.wikipedia.org/wiki/Tobler%27s_hiking_function).

I have indeed used ArcGIS to achieve that in the past, using the 'Path Distance' tool (http://desktop.arcgis.com/en/arcmap/10.3/tools/spatial-analyst-toolbox/how-the-path-distance-tools-work.htm), which implements the method devised by Cixiang Zhan, Sudhakar Menon, Peng Gao (https://link.springer.com/chapter/10.1007/3-540-57207-4_29). I would like to reproduce something similar in R.

Earlier threads

I have made a search in this Forum, but at the best of my knowledge the issue has been always tackled from an ArcGIS perspective.

Issue

In theory, the implementation of the methodology is simple, but I have hard time to devise the correct workflow in R. Also, I do believe that the gdistance package is the key to the calculation of an accumulated cost surface, but I cannot wrap my head around how to properly use the transition() function.

What I did so far

# Define the Tober's function for speed in kmh
tobler.kmh <- function(x){ 6 * exp(-3.5 * abs(tan(x*pi/180) + 0.05)) }

# calculate the slope from an input DTM
SLOPE <- raster::terrain(dtm, opt='slope', unit='degrees', neighbors=8)

# calculate the speed in Kmh
v.kmh <- calc(SLOPE, tobler.kmh)

# time in h to cross 1 km (p stands for 'pace')
p.hkm <- 1/v.kmh    

# express the speed in meters per hour
v.mh <- v.kmh * 1000

# time in hours to cross 1 m
p.hm <- 1 / v.mh

Where I am stuck

Now, with the ultimate goal to devise isochrones, I think I would need to accumulate the pace (in hours per meters = p.hm) outwards from a given location. I guess that the gdistance's accCost() function would make the trick, but I cannot understand how to tackle the intermediate step that entails creating a transition matrix via the transition() function. The latter (at the best of my understanding) would allow to define a traversing cost other than distance in which one incurrs when moving from one cell to an adjacent cell.

I do not know: (a) if the above formula to calculate the pace (in hours per meters) has to be somehow implemented withing the transition() function; (b) if I am approacing the whole issue from the wrong perspective.

Answer 1

I think that this main documentation of the gdistance contains a practical example with that.

If you want to create isochrones, you will need the travel times from one point. You can get that with accCost function and then you can make the time cuts wherever you need.

My approach, following the main documentation is the next:

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

# Calculate speed
adj <- 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 <- geoCorrection(speed, scl=FALSE) # conductance = speed/distance = 1/travel time

See results of traveltime

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

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

This tv object can help you to test if the result of time are correct or almost they are close to the expected.

You can see the conduction from each point with:

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)

And to see in isochrones (just from point A)(I followed you own post to do so :-)):

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

Here I wrote a post about the topic with other examples and layouts. There are a final test of the outputs in order to see if travel times are as expected.

Answer 2

Look into the walkalytics or osrm packages. At the end of the day, what you want is an isochrone. Both the packages use the OpenStreetMap vector layer as a default. But it seems feasible to develop a custom 'profile' for OSRM, using something like the inverse of your walking velocity as the impedance measure of the links. I'm uncertain if raster inputs are permitted, or if you'd need to make it into a vector network first.