Transition Function Inputs

Question

I would like to define my own gdistance::transition function. However, there are two issues that I aren't clear to me yet.

First, I do not yet understand what the input vector to this transition() function looks like. Transition layers based on elevation often use x[2] – x[1]. Others use mean(x), max(x) etc. I assume that the elements of vector x are cell values from the raster layer conditional on the directions option. But which cell is at which position? And what would that function look like if the cost of travelling through a cell only depends on the cell’s value itself, independent of neighbouring cells. Would the function just be x?

Second, how can I feed multiple raster layers into the transition function? Assume that I would like to compute the transition costs between cells based on elevation and temperature. Is there a way to account for the two variables in a single transition function? The documentation indicates that I can use multiple raster bricks but then I don't know how to access the specific elements in my transition function.

The package documentation is illustrative in many regards but unfortunately remains scarce concerning the definition of the actual transition functions. I'm really stuck here.

Answer 1

The answer to your question is indeed not obvious from the documentation. In figuring out what the transition function actually refers to, it helps to look into the source code of that function. In the following I walk you through that code and clarify the steps relevant to your question using an example.

Question 1

setMethod("transition", 
      signature(x = "RasterLayer"), 
      def = function(x, transitionFunction, directions, 
                     symm=TRUE, intervalBreaks=NULL)
    {
        if(class(transitionFunction)=="character") 
        {
            if(transitionFunction != "barriers" & transitionFunction != "areas")
            {
                stop("argument transitionFunction invalid")
            }
            if(transitionFunction=="barriers")
            {
                return(.barriers(x, directions, symm, intervalBreaks))
            }
            if(transitionFunction=="areas")
            {
                return(.areas(x, directions))
            }
        } else {
            return(.TfromR(x, transitionFunction, directions, symm))
        }
    }
)

This chunk of code is in line what the documentation explains. If you feed a single raster layer to transition(), the transitionFunction argument allows for three options: "areas", "barriers" or your own customized transition function. In your case, the customized transition function, the .TfromR() function is called.

.TfromR <- function(x, transitionFunction, directions, symm)
{
    tr <- new("TransitionLayer",
        nrows=as.integer(nrow(x)),
        ncols=as.integer(ncol(x)),
        extent=extent(x),
        crs=projection(x, asText=FALSE),
        transitionMatrix = Matrix(0,ncell(x),ncell(x)),
        transitionCells = 1:ncell(x))
    transitionMatr <- transitionMatrix(tr)
    Cells <- which(!is.na(getValues(x)))
    adj <- adjacent(x, cells=Cells, pairs=TRUE, 
                    target=Cells, 
                    directions=directions)
    if(symm){adj <- adj[adj[,1] < adj[,2],]}
    dataVals <- cbind(getValues(x)[adj[,1]],
                      getValues(x)[adj[,2]])
    transition.values <- apply(dataVals,1,transitionFunction)

    if(!all(transition.values>=0)){
      warning("transition function gives negative values")
    }

    transitionMatr[adj] <- as.vector(transition.values)
    if(symm)
    {
        transitionMatr <- forceSymmetric(transitionMatr)
    }
    transitionMatrix(tr) <- transitionMatr
    matrixValues(tr) <- "conductance"
    return(tr)
}

I am going to illustrate what it does using an unprojected example raster layer of 4 x 4 pixels.

x <- raster(matrix(data = c(10, 20, 7, 8, 15, 18, 16, 5, 1, 21, 3, 15, 22, 17, 12, 14), nrow = 4, ncol = 4, byrow = T))

Then given this input layer, the code first generates an empty transition layer with an empty 16 x 16 sparse transition matrix.

tr <- new("TransitionLayer",
        nrows=as.integer(nrow(x)),
        ncols=as.integer(ncol(x)),
        extent=extent(x),
        crs=projection(x, asText=FALSE),
        transitionMatrix = Matrix(0,ncell(x),ncell(x)),
        transitionCells = 1:ncell(x))
transitionMatr <- transitionMatrix(tr)

transitionMatr
16 x 16 sparse Matrix of class "dsCMatrix"

 [1,] . . . . . . . . . . . . . . . .
 [2,] . . . . . . . . . . . . . . . .
 [3,] . . . . . . . . . . . . . . . .
 [4,] . . . . . . . . . . . . . . . .
 [5,] . . . . . . . . . . . . . . . .
 [6,] . . . . . . . . . . . . . . . .
 [7,] . . . . . . . . . . . . . . . .
 [8,] . . . . . . . . . . . . . . . .
 [9,] . . . . . . . . . . . . . . . .
[10,] . . . . . . . . . . . . . . . .
[11,] . . . . . . . . . . . . . . . .
[12,] . . . . . . . . . . . . . . . .
[13,] . . . . . . . . . . . . . . . .
[14,] . . . . . . . . . . . . . . . .
[15,] . . . . . . . . . . . . . . . .
[16,] . . . . . . . . . . . . . . . .

If you could move from any pixel to any other pixel, all transition matrix cells would receive a value - except for the diagonal along which departure and destination pixel are the same. However, movement through the grid is spatially constrained. You can only move between adjacent cells and therefore only need to calculate transition costs for a small fraction of these pairwise connections. What counts as adjacent is defined by the directions argument in transition(). Here I am going to use directions = 8 (queen's case contiguity).

Transition costs are based on cell values. The function, therefore, checks which cells are not NA (in our case all 16) and then generates a 2-column matrix listing all feasible pairwise connections.

Cells <- which(!is.na(getValues(x)))
Cells
[1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16

adj <- adjacent(x, cells=Cells, pairs=TRUE, 
                target=Cells, 
                directions=directions)
adj
      from to
 [1,]    6  1
 [2,]    7  2
 [3,]    8  3
 [4,]   10  5
 [5,]   11  6
 [6,]   12  7
 [7,]   14  9
 [8,]   15 10
 [9,]   16 11
[10,]    2  1
      ...

The input grid's cells are counted from left to right and from top to bottom. Thus, the first row contains cells 1 to 4, the second row cells 5 to 8 etc. adjacent() lists the available pairwise connections using these cell numbers. The resulting adj matrix has 84 rows in total. The four corner pixels have three neighbors each, the eight edge pixels have five neighbors each and the four center pixels have eight neighbors each (4 * 3 + 8 * 5 + 4 * 8 = 84). Expressed in indices the first two rows of adj states that you can move from cell [2,2] to cell [1,1] and from cell [2,3] to cell [1,2]. The adj matrix list connections in both directions, i.e. from cell 1 to cell 2 and from cell 2 to cell 1. If you assume transition costs to be symmetric (and set symm to true), the length of adj drops to 42 as each bilateral connection is only listed once. Symmetry in your elevation grid would mean that going uphill incurs the same cost as going downhill.

if(symm){adj <- adj[adj[,1] < adj[,2],]}

In this example I assume assymetry (symm == F).

In the next stage, the function generates the dataVals object replacing the cell numbers in adj with the respective cell values.

dataVals <- cbind(getValues(x)[adj[,1]],
                  getValues(x)[adj[,2]])
dataVals
      [,1] [,2]
 [1,]   18   10
 [2,]   16   20
 [3,]    5    7
 [4,]   21   15
 [5,]    3   18
 [6,]   15   16
 [7,]   17    1
 [8,]   12   21
 [9,]   14    3
[10,]   20   10
      ...

And this is how the input to the transition function is generated. The transition function is applied to each row of the dataVals matrix.

transition.values <- apply(dataVals,1,transitionFunction)

The x in your transition function refers to these two elements in each row. x[1] is the value of your departure pixel and x[2] the value of your destination pixel. If we apply your elevation transition function of x[2] - x[1] to this data, the first ten elements of the transition.values vector are: -8, 4, 2, -6, 15, 1, -16, 9, -11, -10. The transition cost refers to the altitude difference, with downhill movements resulting in negative and uphill movements in positive numbers. The corresponding outcomes of mean(x) and max(x) are "14.0, 18.0, 6.0, 18.0, 10.5, 15.5, 9.0, 16.5, 8.5, 15.0" and "18, 20, 7, 21, 18, 16, 17, 21, 14, 20" respectively. Once you know how the input raster enters into your transition function, mofifications become intuitive.

Negative values in a transition layer, as generated by your elevation transition function, can be an issue and lead to errors in other gdistance functions. transition() accordingly prints a warning. And it is advisable to use transition functions that produce positive values only.

if(!all(transition.values>=0)){
  warning("transition function gives negative values")
}

The rest of the code then just plugs the transition values into the transition matrix and the transition matrix into the transition layer. The transition layer is like a raster layer - a matrix with some additional information on projection, extent etc.

transitionMatr[adj] <- as.vector(transition.values)
transitionMatr
16 x 16 sparse Matrix of class "dgCMatrix"

 [1,]   . 10   .  .  5  8  .   .   .  .   .  .  .  .  .  .
 [2,] -10  . -13  . -5 -2 -4   .   .  .   .  .  .  .  .  .
 [3,]   . 13   .  1  . 11  9  -2   .  .   .  .  .  .  .  .
 [4,]   .  .  -1  .  .  .  8  -3   .  .   .  .  .  .  .  .
 [5,]  -5  5   .  .  .  3  .   . -14  6   .  .  .  .  .  .
 [6,]  -8  2 -11  . -3  . -2   . -17  3 -15  .  .  .  .  .
 [7,]   .  4  -9 -8  .  2  . -11   .  5 -13 -1  .  .  .  .
 [8,]   .  .   2  3  .  . 11   .   .  .  -2 10  .  .  .  .
 [9,]   .  .   .  . 14 17  .   .   . 20   .  . 21 16  .  .
[10,]   .  .   .  . -6 -3 -5   . -20  . -18  .  1 -4 -9  .
[11,]   .  .   .  .  . 15 13   2   . 18   . 12  . 14  9 11
[12,]   .  .   .  .  .  .  1 -10   .  . -12  .  .  . -3 -1
[13,]   .  .   .  .  .  .  .   . -21 -1   .  .  . -5  .  .
[14,]   .  .   .  .  .  .  .   . -16  4 -14  .  5  . -5  .
[15,]   .  .   .  .  .  .  .   .   .  9  -9  3  .  5  .  2
[16,]   .  .   .  .  .  .  .   .   .  . -11  1  .  . -2  .

if(symm)
{
    transitionMatr <- forceSymmetric(transitionMatr)
}
transitionMatrix(tr) <- transitionMatr
matrixValues(tr) <- "conductance"
return(tr)

If you set symm to true, the function computes transition costs only for one direction and then pastes the result into both directions in the transition matrix. With a transition function of mean(x) and assumed symmetry the outcome is the following.

transitionMatr
16 x 16 sparse Matrix of class "dsCMatrix"

 [1,]  .   15.0  .    .   12.5 14.0  .    .    .    .    .    .    .    .    .    .  
 [2,] 15.0  .   13.5  .   17.5 19.0 18.0  .    .    .    .    .    .    .    .    .  
 [3,]  .   13.5  .    7.5  .   12.5 11.5  6.0  .    .    .    .    .    .    .    .  
 [4,]  .    .    7.5  .    .    .   12.0  6.5  .    .    .    .    .    .    .    .  
 [5,] 12.5 17.5  .    .    .   16.5  .    .    8.0 18.0  .    .    .    .    .    .  
 [6,] 14.0 19.0 12.5  .   16.5  .   17.0  .    9.5 19.5 10.5  .    .    .    .    .  
 [7,]  .   18.0 11.5 12.0  .   17.0  .   10.5  .   18.5  9.5 15.5  .    .    .    .  
 [8,]  .    .    6.0  6.5  .    .   10.5  .    .    .    4.0 10.0  .    .    .    .  
 [9,]  .    .    .    .    8.0  9.5  .    .    .   11.0  .    .   11.5  9.0  .    .  
[10,]  .    .    .    .   18.0 19.5 18.5  .   11.0  .   12.0  .   21.5 19.0 16.5  .  
[11,]  .    .    .    .    .   10.5  9.5  4.0  .   12.0  .    9.0  .   10.0  7.5  8.5
[12,]  .    .    .    .    .    .   15.5 10.0  .    .    9.0  .    .    .   13.5 14.5
[13,]  .    .    .    .    .    .    .    .   11.5 21.5  .    .    .   19.5  .    .  
[14,]  .    .    .    .    .    .    .    .    9.0 19.0 10.0  .   19.5  .   14.5  .  
[15,]  .    .    .    .    .    .    .    .    .   16.5  7.5 13.5  .   14.5  .   13.0
[16,]  .    .    .    .    .    .    .    .    .    .    8.5 14.5  .    .   13.0  .  

Question 2

The multivariate part of the functions requires the input to be a raster brick. Thus, I generate another 4 x 4 raster layer and merged it to the first one.

y <- raster(matrix(data = c(690, 530, 673, 442, 750, 620, 680, 491, 467, 512, 624, 590, 554, 675, 727, 462), nrow = 4, ncol = 4, byrow = T))
x <- brick(x, y)

The description in the function documentation: "This method serves to summarize several layers of data in a single distance measure. The distance between adjacent cells is the normalized reciprocal of the Mahalanobis distance (mean distance / (mean distance + distance ij)". In the source code it translates to the following.

setMethod("transition", signature(x = "RasterBrick"), 
      def = function(x, transitionFunction="mahal", directions)
    {
        if(transitionFunction != "mahal")
        {
            stop("only Mahalanobis distance method",
                 " implemented for RasterBrick")
        }
        xy <- cbind(1:ncell(x),getValues(x))
        xy <- na.omit(xy)
        dataCells <- xy[,1]
        adj <- adjacent(x, cells=dataCells, pairs=TRUE,
                        target=dataCells, directions=directions)
        x.minus.y <- xy[adj[,1],-1]-xy[adj[,2],-1]
        cov.inv <- solve(cov(xy[,-1]))
        mahaldistance <- apply(x.minus.y,1,function(x){sqrt((x%*%cov.inv)%*%x)})
        mahaldistance <- mean(mahaldistance)/(mahaldistance+mean(mahaldistance))
        transitiondsC <- new("dsCMatrix", 
                p = as.integer(rep(0,ncell(x)+1)),
                Dim = as.integer(c(ncell(x),ncell(x))),
                Dimnames = list(as.character(1:ncell(x)),as.character(1:ncell(x)))
        )
        transitiondsC[adj] <- mahaldistance
        tr <- new("TransitionLayer",
                  nrows=as.integer(nrow(x)),
                  ncols=as.integer(ncol(x)),
                  extent = extent(x),
                  crs=projection(x, asText=FALSE),
                  matrixValues="conductance", 
                  transitionMatrix = transitiondsC)
        return(tr)
    }
)

In the beginning the function follows the same procedure as in the univariate case above. It extracts the cell numbers of those pixels without missing values in any of the input layers. And it computes the adj matrix.

xy <- cbind(1:ncell(x),getValues(x))
xy <- na.omit(xy)
dataCells <- xy[,1]
adj <- adjacent(x, cells=dataCells, pairs=TRUE,
                        target=dataCells, directions=directions)

As there are no missing values in any of the two layers dataCells is equal to Cells in Question 1 and adj is also the same as above.

The function then subtracts the value of the destination cell from the value of the departure cell, separately in each raster layer. This corresponds to a transition function of x[1] - x[2] in the univariate application.

x.minus.y <- xy[adj[,1],-1]-xy[adj[,2],-1]
x.minus.y
      layer.1 layer.2
 [1,]       8     -70
 [2,]      -4     150
 [3,]      -2    -182
 [4,]       6    -238
 [5,]     -15       4
 [6,]      -1     -90
 [7,]      16     208
 [8,]      -9     215
 [9,]      11    -162
[10,]      10    -160
          ...

The first row is derived via 18 - 10 in layer 1 and 620 - 690 in layer 2.

The next lines are a bit technical and refer to what the package authors mention in the documentation: the "normalized reciprocal of the Mahalanobis distance".

cov.inv <- solve(cov(xy[,-1]))
mahaldistance <- apply(x.minus.y,1,function(x){sqrt((x%*%cov.inv)%*%x)})
mahaldistance <- mean(mahaldistance)/(mahaldistance+mean(mahaldistance))

The reast works like in the univariate case. The function generates an empty sparse matrix, plugs the derived transition values into it and plugs that matrix into a transition layer.

transitiondsC <- new("dsCMatrix", 
                p = as.integer(rep(0,ncell(x)+1)),
                Dim = as.integer(c(ncell(x),ncell(x))),
                Dimnames = list(as.character(1:ncell(x)),as.character(1:ncell(x)))
)
transitiondsC[adj] <- mahaldistance

transitiondsC
16 x 16 sparse Matrix of class "dsCMatrix"
   [[ suppressing 16 column names ‘1’, ‘2’, ‘3’ ... ]]

1  .         0.4505593 .         .         0.6804280 0.5638611 .         .         .         .         .         .         .         .         .         .        
2  0.4505593 .         0.4257341 .         0.4451164 0.6629630 0.5342325 .         .         .         .         .         .         .         .         .        
3  .         0.4257341 .         0.4549351 .         0.5119735 0.5834103 0.5193701 .         .         .         .         .         .         .         .        
4  .         .         0.4549351 .         .         .         0.4350499 0.7300553 .         .         .         .         .         .         .         .        
5  0.6804280 0.4451164 .         .         .         0.5752941 .         .         0.3693584 0.4219209 .         .         .         .         .         .        
6  0.5638611 0.6629630 0.5119735 .         0.5752941 .         0.7324411 .         0.4045133 0.6125833 0.4546135 .         .         .         .         .        
7  .         0.5342325 0.5834103 0.4350499 .         0.7324411 .         0.4511992 .         0.5006844 0.4906662 0.6860347 .         .         .         .        
8  .         .         0.5193701 0.7300553 .         .         0.4511992 .         .         .         0.5809755 0.5302805 .         .         .         .        
9  .         .         .         .         0.3693584 0.4045133 .         .         .         0.3874451 .         .         0.3739122 0.3923429 .         .        
10 .         .         .         .         0.4219209 0.6125833 0.5006844 .         0.3874451 .         0.3820783 .         0.8201156 0.5169856 0.4175535 .        
11 .         .         .         .         .         0.4546135 0.4906662 0.5809755 .         0.3820783 .         0.5009822 .         0.4734515 0.5459049 0.4371435
12 .         .         .         .         .         .         0.6860347 0.5302805 .         .         0.5009822 .         .         .         0.5642306 0.6058838
13 .         .         .         .         .         .         .         .         0.3739122 0.8201156 .         .         .         0.5612243 .         .        
14 .         .         .         .         .         .         .         .         0.3923429 0.5169856 0.4734515 .         0.5612243 .         0.6627038 .        
15 .         .         .         .         .         .         .         .         .         0.4175535 0.5459049 0.5642306 .         0.6627038 .         0.4184627
16 .         .         .         .         .         .         .         .         .         .         0.4371435 0.6058838 .         .         0.4184627 .        

tr <- new("TransitionLayer",
          nrows=as.integer(nrow(x)),
          ncols=as.integer(ncol(x)),
          extent = extent(x),
          crs=projection(x, asText=FALSE),
          matrixValues="conductance", 
          transitionMatrix = transitiondsC)
return(tr)

If you would like to use another mathematical approach, you can simply replace the few lines on the normalized reciprocal of the Mahalanobis distance - including x.minus.y - with something else.