Calculate the overlap and non-overlap of paths departing from a common origin. Two algorithms are available: random walk and randomised shortest paths.

pathInc(x, origin, from, to, theta, weight, ...)

Arguments

x

transition matrix of class Transition*

origin

coordinates of the origin (one point location, SpatialPoints, matrix or numeric class)

from

coordinates of the destinations (SpatialPoints, matrix or numeric class)

to

second set of coordinates of the destinations (can be missing)

theta

value > 0 and < 20 (randomised shortest paths) or missing (random walk)

weight

matrix – Reciprocals of the non-zero values are used as weights. If missing, reciprocals of the transition matrix are used.

...

additional arguments passed to methods. See Details.

Value

list of dist objects or list of matrices

Details

This is a complex wrapper function that calculates to what extent dispersal routes overlap or do not overlap.

First, the function calculates the trajectories for all "from" and "to" locations, starting from a single "origin" location. These trajectories can either be based on random walks or randomised shortest paths (giving a value to theta).

Then, for all unique pairs of trajectories, it calculates the extent to which these trajectories overlap or diverge. These values are given back to the user as a list of (distance) matrices.

If only "from" coordinates are given, the function calculates symmetric distance matrices for all combinations of points in "from". If both "from" and "to" coordinates are given, the function calculates a matrix of values with rows for all locations in "from" and columns for all locations in "to".

Overlap is currently calculated as the minimum values of each pair of trajectories compared. Non-overlap uses the following formula: Nonoverlap = max(0,max(a,b)*(1-min(a,b))-min(a,b)) (see van Etten and Hijmans 2010). See the last example to learn how to use an alternative function.

References

McRae B.H., B.G. Dickson, and T. Keitt. 2008. Using circuit theory to model connectivity in ecology, evolution, and conservation. Ecology 89:2712-2724.

Saerens M., L. Yen, F. Fouss, and Y. Achbany. 2009. Randomized shortest-path problems: two related models. Neural Computation, 21(8):2363-2404.

van Etten, J., and R.J. Hijmans. 2010. A geospatial modelling approach integrating archaeobotany and genetics to trace the origin and dispersal of domesticated plants. PLoS ONE 5(8): e12060.

Author

Jacob van Etten. Implementation of randomised shortest paths based on Matlab code by Marco Saerens.

Examples

library("raster")
library("sp")

# Create TransitionLayer
r <- raster(ncol=36,nrow=18)
r <- setValues(r,rep(1,times=ncell(r)))
tr <- transition(r,mean,directions=4)
#> The extent and CRS indicate this raster is a global lat/lon raster. This means that transitions going off of the East or West edges will 'wrap' to the opposite edge.
#> Global lat/lon rasters are not supported under new optimizations for 4 and 8 directions with custom transition functions. Falling back to old method.

# Two different types of correction are required
trR <- geoCorrection(tr, type="r", multpl=FALSE)
trC <- geoCorrection(tr, type="c", multpl=FALSE)

# Create TransitionStack
ts <- stack(trR, trR)

# Points for origin and coordinates between which to calculate path (non)overlaps
sP0 <- SpatialPoints(cbind(0,0))
sP1 <- SpatialPoints(cbind(c(65,5,-65),c(-55,35,-35)))

# Randomised shortest paths
# rescaling is needed: exp(-theta * trC) should give reasonable values
# divide by median of the non-zero values
trC <- trC / median(transitionMatrix(trC)@x)
pathInc(trC, origin=sP0, from=sP1, theta=2)
#> $function1layer1
#>             1           2
#> 2 0.009152689            
#> 3 0.690511868 0.008950098
#> 
#> $function2layer1
#>          1        2
#> 2 13.37958         
#> 3 16.94227 13.33754
#> 

# Random walk
pathInc(trR, origin=sP0, from=sP1)
#> $function1layer1
#>          1        2
#> 2 19864845         
#> 3 27837551 19721595
#> 
#> $function2layer1
#>          1        2
#> 2 36590892         
#> 3 34131717 34685800
#> 

# TransitionStack as weights
pathInc(trR, origin=sP0, from=sP1, weight=ts)
#> $function1layer1
#>          1        2
#> 2 19864845         
#> 3 27837551 19721595
#> 
#> $function1layer2
#>          1        2
#> 2 19864845         
#> 3 27837551 19721595
#> 
#> $function2layer1
#>          1        2
#> 2 36590892         
#> 3 34131717 34685800
#> 
#> $function2layer2
#>          1        2
#> 2 36590892         
#> 3 34131717 34685800
#> 

# Demonstrate use of an alternative function
#The current default is to take the minimum of each pair of layers

altoverlap <- function(a, b)
{
  aV <- as.vector(a[,rep(1:ncol(a), each=ncol(b))])
  bV <- as.vector(b[,rep(1:ncol(b), times=ncol(a))])
  result <- matrix(aV * bV, nrow = nrow(a), ncol=ncol(a)*ncol(b))
  return(result)
}

pathInc(trR, origin=sP0, from=sP1, weight=ts, functions=list(altoverlap))
#> $function1layer1
#>         1       2
#> 2 1281948        
#> 3 1547342 1317062
#> 
#> $function1layer2
#>         1       2
#> 2 1281948        
#> 3 1547342 1317062
#>