7

Background and Aim

I would like to create a cumulative raster representing walking time across a surface.

For my purposes, the cost to be accumulated has to be the walking time as influenced by terrain slope, according to the Tobler's hiking function below (which has to be be modified to express time instead of speed):

enter image description here

The corner stone for what I am after is a raster representing terrain elevation (let's call it dtm).

Where I am stuck

I seek some elucidations on the results I got in the very preliminary steps I took to move toward my goal. In particular, I am at loss of understanding the values of what are called transition layers produced by the gdistance package in R.

I am going to provide a reproducible example of what I did (working with a small raster representing terrain elevation), and I will finally point out what I am having difficulties to grasp.

Reproducible example

(1) Load a sample dataset, which will be subsetted to have a smaller raster to work with:

library(gdistance)
dtm <- raster(system.file("external/maungawhau.grd", package="gdistance"))
subR <- crop(dtm, extent(dtm, 30, 36, 50, 56))

enter image description here

(2) Calculate the slope using gdistance; first, the function altDiff defines the difference in elevation between adjacent cells; secondly, the function is applied to the subR raster (elevation) via the transition() command, considering moves in 8 directions; finally, geoCorrection divides the differences by the inter-cells distance, eventually arriving at (if I am not mistaken) slope expressed as rise/run:

altDiff <- function(x){x[2] - x[1]}
hd <- transition(subR, altDiff, 8, symm=FALSE)
slope <- geoCorrection(hd,scl=FALSE)
plot(raster(slope), sub="slope (rise over run)", cex.sub=0.8)

enter image description here

(3) The code below calculates the pace (hours per meter) by reworking the Tobler's function and applying the latter to the slope raster; multipling by 1000 and taking the reciprocal helps turning speed (kmh) into pace (hours per meter) (please, disregard the error in the image subtitle, which shows kmh instead of hours per meter).

tobler.hm <- function(x){ 1 / ((6 * exp(-3.5 * abs(x + 0.05))) * 1000) }
adj <- adjacent(x=subR, cells=1:ncell(subR), direction=8)
pace <- slope
pace[adj] <- tobler.kmh(slope[adj]) 
plot(raster(pace), sub="m per hour", cex.sub=0.8)
plotvalues(raster(pace))

enter image description here

The above would be the cost surface (expressing pace, i.e. meters per hour) that should be accumulated.

Question

My issue at this point is to understand the values stored in the two above rasters (slope and pace) devised using the transion() function. While I understand that, by using the transition() function,

"users define a function f(i,j) to calculate the transition value for each pair of adjacent cells i and j" (quoting from https://cran.r-project.org/web/packages/gdistance/vignettes/gdistance1.pdf),

I am at loss of understanding how those values (in each respective raster) related to each other? For instance, in the slope raster, how does the 0.291 (upper-left cell) value relate to the 0.102 immediately at its right, and to the 0.249 laying immediately below? The same applies to the values in the pace raster. I need to grasp this (possibly very basic) step in order for me to perform additional calculation to eventually accumulate the pace values from a starting location.

3

The results of these functions are not trivial to mentally visualise as they are based on transition matrices, but consider this toy e.g.:

N.B. i'm using a projected coordinate system, British National Grid, to simplify things a bit as geoCorrection will begin to take into account north/south & east/west distortions of latlong rasters by using great circle calculations etc.

# 1) make this raster
library("gdistance")
BNG <- CRS("+init=epsg:27700")

r <- raster(xmn=0,xmx=4,ymn=0,ymx=4,res=1,crs=BNG)
r[] <- c(2,2,3,3,4,3,2,2,2,4,1,2,1,3,3,1)
plot(r)

enter image description here

# 2) define altitude difference function and create a transition matrix from the original
altDiff <- function(x){x[2] - x[1]}
hd <- transition(r, altDiff, 8, symm=FALSE)
plot(raster(hd), sub="non geo-corr", cex=0.8)
text(raster(hd),digits=3)

enter image description here To obtain the values in this average altitude difference raster above, the 'conceptual' way is that for a given cell x, you subtract values of neighbouring cells one by one (y1, y2..., sigma and take a mean (n = number neighbouring cells != x). As just stated in the brackets, if the value of a neighbour is the same, do not use in the sigma or the mean. E.g. for cell 5 above: ((4-2)+(4-2)+(4-3)+(4-2))/4 = 1.75. (There are 5 neighbours, not 8, and 1 is discarded)

# 3) This is then geoCorrected to slope, i.e. altitude diff / distance travelled
slope <- geoCorrection(hd,scl=FALSE)
plot(raster(slope), sub="slope (height diff / distance))", cex.sub=0.8)
text(raster(slope),digits=3)

enter image description here How did the cells change to these new values? The non-corrected values now have distance taken into account as a function, i.e. the diagonals, so this is a mean slope to this cell, from the surrounding 8, taking into account distance.

Whereas before, in the cell 5 example, we did: ((4-2)+(4-2)+(4-3)+(4-2))/4 = 1.75, we now adjust for distance (by simple pythagoras, the diagonal distance is sqrt(2) );

( ((4-2)*1/1) + ((4-2)*1/sqrt(2)) + ((4-3)*1/1) + ((4-2)*1/1) )/4 = 1.604

# 4) define Tobler function, get adjacent cell pairs and copy slope to pace
tobler.hm <- function(x){ 1 / ((6 * exp(-3.5 * abs(x + 0.05))) * 1000) }
adj <- adjacent(x=r, cells=1:ncell(r), direction=8)
pace <- slope

adj is the 84 direction-dependent cell combinations (1 to 2, 10 to 14 etc.). adj will always be ncells x 8 - (4 corners cells x 5) - (remaining border cells x 3). That is 16 cells x 8 directions minus 44 border cell pairings that don't exist.

# 5) The following line will apply Tobler's function *to the matrix of slope*, not
# the raster values you see when plotting slope, to calculate pace
pace[adj] <- tobler.hm(slope[adj]) 
plot(raster(pace), sub="m per hour", cex.sub=0.8)
text(raster(slope),digits=3)

enter image description here All the values used to calculate slope in step 3 have had the Tobler function applied, and a new surface created. This time you DO include the cells of same value as you are computing average pace to that cell, not slope. Conceptually:

(tobler.hm(((4-2)*1/1)) + tobler.hm(((4-2)*1/sqrt(2))) + 
tobler.hm(((4-3)*1/1)) + tobler.hm(((4-4)*1/1)) + tobler.hm(((4-2)*1/1)) ) / 5 = 0.09404

Inspecting the transition matrix values is a big help in these stages. Values in the raster represent an average (which average is dependent on the stage) of the values 'coming into' that cell

# inspect a transition matrix for slope, for example
round(transitionMatrix(slope),3)
16 x 16 sparse Matrix of class "dgCMatrix"

 [1,]  .      .     .     .     2.000  0.707  .      .      .     .      .      .      .      .      .      .
 [2,]  .      .     1.000 .     1.414  1.000  .      .      .     .      .      .      .      .      .      .
 [3,]  .     -1.000 .     .     .      .     -1.000 -0.707  .     .      .      .      .      .      .      .
 [4,]  .      .     .     .     .      .     -0.707 -1.000  .     .      .      .      .      .      .      .
 [5,] -2.000 -1.414 .     .     .     -1.000  .      .     -2.000 .      .      .      .      .      .      .
 [6,] -0.707 -1.000 .     .     1.000  .     -1.000  .     -0.707 1.000 -1.414  .      .      .      .      .
 [7,]  .      .     1.000 0.707 .      1.000  .      .      .     1.414 -1.000  .      .      .      .      .
 [8,]  .      .     0.707 1.000 .      .      .      .      .     .     -0.707  .      .      .      .      .
 [9,]  .      .     .     .     2.000  0.707  .      .      .     2.000  .      .     -1.000  0.707  .      .
[10,]  .      .     .     .     .     -1.000 -1.414  .     -2.000 .     -3.000  .     -2.121 -1.000 -0.707  .
[11,]  .      .     .     .     .      1.414  1.000  0.707  .     3.000  .      1.000  .      1.414  2.000  .
[12,]  .      .     .     .     .      .      .      .      .     .     -1.000  .      .      .      0.707 -1
[13,]  .      .     .     .     .      .      .      .      1.000 2.121  .      .      .      2.000  .      .
[14,]  .      .     .     .     .      .      .      .     -0.707 1.000 -1.414  .     -2.000  .      .      .
[15,]  .      .     .     .     .      .      .      .      .     0.707 -2.000 -0.707  .      .      .     -2
[16,]  .      .     .     .     .      .      .      .      .     .      .      1.000  .      .      2.000  .

It is perhaps easier to visualise for cell 5, i.e. column 5, the values to mean and that have the Tobler function applied to.

Finally, the pace matrix will then need geoCorrecting again to account for distance to create a time matrix. See the doc linked below.

Hope this all helps. I'm guessing section 9 of this document is your friend.

  • Wow! Thanks for taking the time to provide me with such an extensive reply. It will take some time to go through it. I will comment back in case I need further elucidations. – NewAtGis May 26 '18 at 9:36
  • no problem, let me know how you get on – Sam Jun 6 '18 at 8:27
  • I went through your well-explained reply, and I could follow most of it; some aspects are above my head. My bad. I read again the gdistance guide, but I am still puzzled by the conductance values eventually returned by the described routine (page 14 of the PDF). Those values are 1/travel time. Since my ultimate goal is to calculate travel time and plot isochrones, I am still wondering how can I adapt to code to get what I am after (travel time). I have also another question (i.e. accounting for true topographic distance between cells) which I will possibly ask in a different thread. – NewAtGis Jun 13 '18 at 7:31
  • I re-read your reply, and it is almost all clear now. Thanks again. There is something I do not fully grasp in the very last bit of your reply. It is all clear till the step in which we apply the reciprocal of the Tobler's function to the matrix of slope (i.e. slope[adj]) to get the matrix of pace. So far so good. Then you say that we need to apply geo-corretion to the pace matrix. This is what is haunting me. As I read in the package vignette (page 14), the geocorrection will divide the pace values by the inter-cell distance. Are we sure that what we will obtain will be still the pace? – NewAtGis Jun 26 '18 at 16:51
  • we know that speed is space/time, and pace (being its reciprocal) is time/space. If we divide (as geocorrection does if we apply geocorrection to the pace matrix) time/space (i.e., pace) by space (i.e., inter-cell distance) we will end up having time/space^2 (if my math is not wrong). This is what is haunting me. Maybe I am wrong. Any elucidation will be appreciated. – NewAtGis Jun 26 '18 at 16:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.