# Understanding the values from transition layers produced by the R package `gdistance`

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): 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))
`````` (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 - x}
hd <- transition(subR, altDiff, 8, symm=FALSE)
slope <- geoCorrection(hd,scl=FALSE)
plot(raster(slope), sub="slope (rise over run)", cex.sub=0.8)
`````` (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) }
pace <- slope
plot(raster(pace), sub="m per hour", cex.sub=0.8)
plotvalues(raster(pace))
`````` 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.

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)
`````` ``````# 2) define altitude difference function and create a transition matrix from the original
altDiff <- function(x){x - x}
hd <- transition(r, altDiff, 8, symm=FALSE)
plot(raster(hd), sub="non geo-corr", cex=0.8)
text(raster(hd),digits=3)
`````` 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)
`````` 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) }
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
plot(raster(pace), sub="m per hour", cex.sub=0.8)
text(raster(slope),digits=3)
`````` 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. 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. 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? 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. Jun 26 '18 at 16:59