,
Kauê de Sousa [ctb] ,
Andrew Marx [cre, ctb] | Title: | Distances and Routes on Geographical Grids |
|---|---|
| Description: | Provides classes and functions to calculate various distance measures and routes in heterogeneous geographic spaces represented as grids. The package implements measures to model dispersal histories first presented by van Etten and Hijmans (2010) <doi:10.1371/journal.pone.0012060>. Least-cost distances as well as more complex distances based on (constrained) random walks can be calculated. The distances implemented in the package are used in geographical genetics, accessibility indicators, and may also have applications in other fields of geospatial analysis. |
| Authors: | Jacob van Etten [aut] ,
Kauê de Sousa [ctb] ,
Andrew Marx [cre, ctb] |
| Maintainer: | Andrew Marx <[email protected]> |
| License: | GPL (>=2) |
| Version: | 1.6.4 |
| Built: | 2024-08-31 03:20:19 UTC |
| Source: | https://github.com/agrdatasci/gdistance |
Calculates the accumulated cost surface from one or more origins.
accCost(x, fromCoords)accCost(x, fromCoords)
x |
object of class |
fromCoords |
origin point locations (SpatialPoints, matrix or numeric class) |
If more than one coordinate is supplied in fromCoords, the function calculates the minimum least-cost distance from any origin point.
The function uses Dijkstra's algorithm (as implemented in the igraph package).
a RasterLayer object
Jacob van Etten
E.W. Dijkstra. 1959. A note on two problems in connexion with graphs. Numerische Mathematik 1, 269 - 271.
library("raster") # example equivalent to that in the documentation on r.cost in GRASS r <- raster(nrows=6, ncols=7, xmn=0, xmx=7, ymn=0, ymx=6, crs="+proj=utm +units=m") r[] <- c(2, 2, 1, 1, 5, 5, 5, 2, 2, 8, 8, 5, 2, 1, 7, 1, 1, 8, 2, 2, 2, 8, 7, 8, 8, 8, 8, 5, 8, 8, 1, 1, 5, 3, 9, 8, 1, 1, 2, 5, 3, 9) # 1/mean: reciprocal to get permeability tr <- transition(r, function(x) 1/mean(x), 8) tr <- geoCorrection(tr) c1 <- c(5.5,1.5) c2 <- c(1.5,5.5) A <- accCost(tr, c1) plot(A) text(A)library("raster") # example equivalent to that in the documentation on r.cost in GRASS r <- raster(nrows=6, ncols=7, xmn=0, xmx=7, ymn=0, ymx=6, crs="+proj=utm +units=m") r[] <- c(2, 2, 1, 1, 5, 5, 5, 2, 2, 8, 8, 5, 2, 1, 7, 1, 1, 8, 2, 2, 2, 8, 7, 8, 8, 8, 8, 5, 8, 8, 1, 1, 5, 3, 9, 8, 1, 1, 2, 5, 3, 9) # 1/mean: reciprocal to get permeability tr <- transition(r, function(x) 1/mean(x), 8) tr <- geoCorrection(tr) c1 <- c(5.5,1.5) c2 <- c(1.5,5.5) A <- accCost(tr, c1) plot(A) text(A)
Identify pairs of cells that are adjacent.
adjacencyFromTransition(x)adjacencyFromTransition(x)
x |
TransitionLayer |
Extracts the indices of those cells that are connected (e.g. cells i,j that have a non-zero value in the transition matrix).
Cell numbers are unique indices of cells in the original grid. By convention, cell numbers start with 1 in the upper-left corner of the grid and increase from left to right and from top to bottom.
A two column matrix with each row containing a pair of adjacent cells.
Jacob van Etten
library("raster") r <- raster(nrows = 6, ncols = 7, xmn = 0, xmx = 7, ymn = 0, ymx = 6, crs = "+proj=utm +units=m") r[] <- runif(6 * 7) tr <- transition(r, function(x) 1 / mean(x), 8) aft <- adjacencyFromTransition(tr) head(aft)library("raster") r <- raster(nrows = 6, ncols = 7, xmn = 0, xmx = 7, ymn = 0, ymx = 6, crs = "+proj=utm +units=m") r[] <- runif(6 * 7) tr <- transition(r, function(x) 1 / mean(x), 8) aft <- adjacencyFromTransition(tr) head(aft)
Standard arithmetic operators for computations with Transition*
objects and numeric values. Transition objects must have the same extent
and resolution. All arithmetic and mathematical operations that work on
the sparse matrices are available for Transition* objects.
## S4 method for signature 'TransitionLayer,TransitionLayer' Arith(e1, e2)## S4 method for signature 'TransitionLayer,TransitionLayer' Arith(e1, e2)
e1 |
objects |
e2 |
objects |
a Transition* object or numeric.
Jacob van Etten
#create a new raster and set all its values to unity. raster <- raster(nrows=18, ncols=36) raster <- setValues(raster,rep(1,ncell(raster))) #create TransitionLayer objects tr1 <- transition(raster,mean,4) tr2 <- tr1 #arithmetic operations tr3 <- tr1 * tr2 tr4 <- tr3 * 4 #mathematical operations tr5 <- sqrt(tr4)#create a new raster and set all its values to unity. raster <- raster(nrows=18, ncols=36) raster <- setValues(raster,rep(1,ncell(raster))) #create TransitionLayer objects tr1 <- transition(raster,mean,4) tr2 <- tr1 #arithmetic operations tr3 <- tr1 * tr2 tr4 <- tr3 * 4 #mathematical operations tr5 <- sqrt(tr4)
Calculates the resistance distance between points.
commuteDistance(x, coords)commuteDistance(x, coords)
x |
object of class |
coords |
point locations coordinates (of SpatialPoints, matrix or numeric class) |
This function calculates the expected random-walk commute time between nodes in a graph. It is defined as the effective distance (resistance distance) between the selected nodes multiplied by the volume of the graph, which is the sum of the conductance weights of all the edges in the graph (Chandra et al. 1997). The result represents the average number of steps that is needed to commute between the nodes during a random walk.
The function implements the algorithm given by Fouss et al. (2007).
Before calculating commute-time distances from a TransitionLayer
object, see if you need to apply the function
geoCorrection
distance matrix (S3 class dist or matrix)
Jacob van Etten
Chandra, A.K., Raghavan, P., Ruzzo, W.L., Smolensy, R. & Tiwari, P. 1996. The electrical resistance of a graph captures its commute and cover times. Computational Complexity, 6(4), 312-340.
Fouss, F., Pirotte, A., Renders, J.-M. & Saerens, M. 2007. Random-walk computation of similarities between nodes of a graph with application to collaborative recommendation. IEEE Transactions on Knowledge and Data Engineering, 19(3), 355-369.
McRae, B.H. 2006. Isolation by resistance. Evolution 60(8), 1551-1561. https://circuitscape.org/
library("raster") # Create a new raster and set all its values to unity. r <- raster(nrows=18, ncols=36) r <- setValues(r,rep(1,ncell(raster))) #Create a Transition object from the raster tr <- transition(r, function(x) 1/mean(x),4) # Create two sets of coordinates library("sp") sP1 <- SpatialPoints(cbind(c(65,5,-65),c(55,35,-35))) sP2 <- SpatialPoints(cbind(c(50,15,-40),c(80,20,-5))) #Calculate the resistance distance between the points commuteDistance(tr, sP1)library("raster") # Create a new raster and set all its values to unity. r <- raster(nrows=18, ncols=36) r <- setValues(r,rep(1,ncell(raster))) #Create a Transition object from the raster tr <- transition(r, function(x) 1/mean(x),4) # Create two sets of coordinates library("sp") sP1 <- SpatialPoints(cbind(c(65,5,-65),c(55,35,-35))) sP2 <- SpatialPoints(cbind(c(50,15,-40),c(80,20,-5))) #Calculate the resistance distance between the points commuteDistance(tr, sP1)
This is a class union, providing a convenient class for coordinates in several formats.
The class accepts coordinates in any of the following formats: 1. SpatialPoints 2. two-columned matrix 3. vector of length 2
showClass("Coords")showClass("Coords")
Calculate the least-cost distance between points.
costDistance(x, fromCoords, toCoords)costDistance(x, fromCoords, toCoords)
x |
object of class |
fromCoords |
first set of point locations
(of |
toCoords |
optional, second set of point
locations (of |
Cost units between cells are defined as the reciprocal of the values in the transition matrix.
The function uses Dijkstra's algorithm, as implemented in the igraph package.
A projection correction is needed for accuracy in the case of grid data for a longlat raster (see function geoCorrection).
distance matrix (S3 class dist or matrix)
Jacob van Etten
E.W. Dijkstra. 1959. A note on two problems in connexion with graphs. Numerische Mathematik 1, 269 - 271.
library("raster") # create a new raster and set all its values to unity. r <- raster(nrows=18, ncols=36) r <- setValues(r,runif(ncell(r),0,1)) # create a Transition object from the raster tr <- transition(r,function(x) 1/mean(x),8) # asymmetric ncf <- function(x) max(x) - x[1] + x[2] tr2 <- transition(r,ncf,8, symm=FALSE) # create two sets of coordinates sP1 <- cbind(c(65,5,-65),c(55,35,-35)) sP2 <- cbind(c(50,15,-40),c(80,20,-5)) # from and to identical costDistance(tr,sP1) costDistance(tr2,sP1) # from and to different costDistance(tr,sP1,sP2) costDistance(tr2,sP1,sP2)library("raster") # create a new raster and set all its values to unity. r <- raster(nrows=18, ncols=36) r <- setValues(r,runif(ncell(r),0,1)) # create a Transition object from the raster tr <- transition(r,function(x) 1/mean(x),8) # asymmetric ncf <- function(x) max(x) - x[1] + x[2] tr2 <- transition(r,ncf,8, symm=FALSE) # create two sets of coordinates sP1 <- cbind(c(65,5,-65),c(55,35,-35)) sP2 <- cbind(c(50,15,-40),c(80,20,-5)) # from and to identical costDistance(tr,sP1) costDistance(tr2,sP1) # from and to different costDistance(tr,sP1,sP2) costDistance(tr2,sP1,sP2)
Provides classes and functions to calculate various distance measures and routes in heterogeneous geographic spaces represented as grids. The package implements measures to model dispersal histories first presented by van Etten and Hijmans (2010) doi:10.1371/journal.pone.0012060. Least-cost distances as well as more complex distances based on (constrained) random walks can be calculated. The distances implemented in the package are used in geographical genetics, accessibility indicators, and may also have applications in other fields of geospatial analysis.
Jacob van Etten and Kauê de Sousa
Useful links:
Package's paper: doi:10.18637/jss.v076.i13
Development repository: https://github.com/AgrDataSci/gdistance
Static documentation: https://github.com/AgrDataSci/gdistance
Report bugs: https://github.com/AgrDataSci/gdistance/issues
genDist: Genetic distances between 23 populations of Y-chromosome haplogroup R1b1b2 in Europe, calculated from haplotype microsatellite data (see source).
dist object (genDist) and data.frame (popCoord)
The distance measure used is the negative logarithm of the shared proportion of alleles.
popCoord: Geographic coordinates of the same populations.
Distances calculated from Table S1 in Balaresque et al. (2010). Coordinates from Table 1.
Balaresque P., et al. 2010. A predominantly Neolithic origin for European paternal lineages. PLoS Biology 8(1): e1000285.
Correct 'TransitionLayer' objects taking into account local distances
geoCorrection(x, type, ...)geoCorrection(x, type, ...)
x |
object of class |
type |
type of correction: "c", "r", or missing (only required for lonlat, see Details) |
... |
additional arguments passed to methods. |
Geographic correction is necessary for all objects of the class Transition that are either: (1) based on a grid in a geographic (lonlat) projection and covering a large area; (2) made with directions > 4.
The function will correct for map distortion, as well as for diagonal connections between grid cells (which cover a longer distance than vertical or horizontal connections).
When working with lonlat grids, users should also anticipate whether they
will use methods based on either least-cost or random walks, and set the
type argument accordingly. In the case of least-cost distances, the
correction is only done in East-West direction. In the case of random walks
there should be an additional correction which reduces the conductance in
North-South direction (type="r").
The correction is done by dividing conductance values by the inter-cell
distance. Distances are calculated as great-circle distances for lonlat
grids (see function [raster]{isLonLat}) and Euclidean distances for all other grids.
In the case of randomised shortest paths, the need for correction is somewhat in between these two correction methods. We have not developed an analytical solution for this problem yet. With very low values for theta, you may choose the correction for random walks, and for high values the one for least-cost paths. Users who want to work with intermediate values of theta are encouraged to experiment with different solutions.
The values are scaled to get values near 1 if the argument scl is set to
TRUE. This is desirable for subsequent calculations involving random walk
calculations. Values are divided by the W-E inter-cell distance (at the
centre of the grid).
a 'Transition*' object
Jacob van Etten
library("raster") r <- raster(ncol=36,nrow=18) r <- setValues(r,rep(1,times=ncell(r))) tr <- transition(r, mean, directions=8) # directly tr1 <- geoCorrection(tr, type="c", multpl=FALSE) tr1 # the same, but with a separate correction matrix trCorr <- geoCorrection(tr, type="c", multpl=TRUE) tr2 <- tr * trCorr tr2library("raster") r <- raster(ncol=36,nrow=18) r <- setValues(r,rep(1,times=ncell(r))) tr <- transition(r, mean, directions=8) # directly tr1 <- geoCorrection(tr, type="c", multpl=FALSE) tr1 # the same, but with a separate correction matrix trCorr <- geoCorrection(tr, type="c", multpl=TRUE) tr2 <- tr * trCorr tr2
Normalize the transition matrix.
normalize(x, ...)normalize(x, ...)
x |
object of class |
... |
optional argument |
argument method passed through ...
a character for the normalization method
the default is 'row' users can set the optional method argument
to either "col" or "symm"
Normalization of the weighted adjacency matrix in the Transition* object. Matrix values are divided by their respective row-sums, column-sums, or the product of the square-roots of both (symmetric normalization).
The default method is row-normalization. To use the other normalization
methods, users can set the optional method argument
to either "col" or "symm".
For random walk calculations a symmetric matrix is needed (method = "symm").
an object of class TransitionLayer
Jacob van Etten
von Luxburg, U. 2007. A tutorial on spectral clustering. Statistics and Computing 17(4), 395-416. doi:10.1007/s11222-007-9033-z
Chung, F. 1997. Spectral Graph Theory. Conference Board of the Mathematical Sciences, Washington.
library("raster") r <- raster(ncol=36,nrow=18) r <- setValues(r,rep(1,times=ncell(r))) tr <- transition(r, mean, directions=8) normalize(tr) normalize(tr, method="symm")library("raster") r <- raster(ncol=36,nrow=18) r <- setValues(r,rep(1,times=ncell(r))) tr <- transition(r, mean, directions=8) normalize(tr) normalize(tr, method="symm")
Special functions to calculate the degree of overlap and nonoverlap between trajectories
overlap(a, b)overlap(a, b)
a |
|
b |
|
These functions are used by the pathInc() as defaults.
A matrix object
The functions are provided here to assist the user in defining alternative measures of overlap / nonoverlap.
Jacob van Etten
Calculates for each cell the number of passages of a random-walk or randomised shortest paths with given origin(s) and destination(s). Either the total or the net number of passages can be calculated. In the case of multiple origins or destinations, each receives equal weight.
passage(x, origin, goal, theta, ...)passage(x, origin, goal, theta, ...)
x |
Object of class |
origin |
|
goal |
|
theta |
If zero or missing, a random walk results. If a numeric value 0 < theta < 20 is given, randomised shortest paths are calculated; theta is the degree from which the path randomly deviates from the shortest path |
... |
Additional arguments: totalNet ("total" or "net"), and output ("RasterLayer" or "Transition") |
The net number of passages between i and j is defined as: abs( passages from i to j - passages from j to i ).
Defaults for additional argument totalNet is "net"
and for output it is "RasterLayer".
Random walk requires a symmetric transition matrix.
RasterLayer or Transition object, depending on the output argument
Jacob van Etten. Implementation of randomised shortest paths based on Matlab code by Marco Saerens
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.
#create a new raster and set all its values to unity. raster <- raster(nrows=18, ncols=36) raster <- setValues(raster,rep(1,ncell(raster))) #create a Transition object from the raster tr <- transition(raster,mean,4) trC <- geoCorrection(tr, type="c", scl=TRUE) trR <- geoCorrection(tr, type="r", scl=TRUE) #create two coordinates sP1 <- SpatialPoints(cbind(-105,55)) sP2 <- SpatialPoints(cbind(105,-55)) #randomised shortest paths with theta = 2 rSPraster <- passage(trC, sP1, sP2, 2) plot(rSPraster) points(sP1) points(sP2) #randomised shortest paths with theta = 0.05 rSPraster <- passage(trC, sP1, sP2, 0.05) plot(rSPraster) points(sP1) points(sP2) #randomised shortest paths with theta = 0.05 #and total rSPraster <- passage(trC, sP1, sP2, 0.05, totalNet = "total") plot(rSPraster) points(sP1) points(sP2) #random walk rwraster <- passage(trR, sP1, sP2) plot(rwraster) points(sP1) points(sP2)#create a new raster and set all its values to unity. raster <- raster(nrows=18, ncols=36) raster <- setValues(raster,rep(1,ncell(raster))) #create a Transition object from the raster tr <- transition(raster,mean,4) trC <- geoCorrection(tr, type="c", scl=TRUE) trR <- geoCorrection(tr, type="r", scl=TRUE) #create two coordinates sP1 <- SpatialPoints(cbind(-105,55)) sP2 <- SpatialPoints(cbind(105,-55)) #randomised shortest paths with theta = 2 rSPraster <- passage(trC, sP1, sP2, 2) plot(rSPraster) points(sP1) points(sP2) #randomised shortest paths with theta = 0.05 rSPraster <- passage(trC, sP1, sP2, 0.05) plot(rSPraster) points(sP1) points(sP2) #randomised shortest paths with theta = 0.05 #and total rSPraster <- passage(trC, sP1, sP2, 0.05, totalNet = "total") plot(rSPraster) points(sP1) points(sP2) #random walk rwraster <- passage(trR, sP1, sP2) plot(rwraster) points(sP1) points(sP2)
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, ...)pathInc(x, origin, from, to, theta, weight, ...)
x |
transition matrix of class |
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. |
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.
list of dist objects or list of matrices
Jacob van Etten. Implementation of randomised shortest paths based on Matlab code by Marco Saerens.
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.
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) # 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) # Random walk pathInc(trR, origin=sP0, from=sP1) # TransitionStack as weights pathInc(trR, origin=sP0, from=sP1, weight=ts) # 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))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) # 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) # Random walk pathInc(trR, origin=sP0, from=sP1) # TransitionStack as weights pathInc(trR, origin=sP0, from=sP1, weight=ts) # 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))
Create a RasterLayer from a TransitionLayer with a call to
the generic function raster. The n x n transition matrix of
the TransitionLayer is transformed to form the values n cells of a raster.
## S4 method for signature 'TransitionLayer' raster(x, reduceMethod = "NZcolMeans")## S4 method for signature 'TransitionLayer' raster(x, reduceMethod = "NZcolMeans")
x |
an object of class |
reduceMethod |
character for the method to reduce the transition matrix. See details |
The following methods to ‘reduce’ the transition matrix are
available with the optional argument reduceMethod):
colSums
rowSums
colMeans
rowMeans
NZcolMeans
NZrowMeans
The latter two methods only take into account the non-zero entries in the transition matrix. The default is NZcolMeans.
a RasterLayer
Jacob van Etten
#create a new raster and set all its values to unity. r <- raster(nrows=18, ncols=36) r <- setValues(r,runif(ncell(r),0,1)) #create a Transition object from the raster tr1 <- transition(r,mean,8) #asymmetric asf <- function(x) max(x) - x[1] + x[2] tr2 <- transition(r,asf,8, symm=FALSE) #create RasterLayer objects r1 <- raster(tr1) r2 <- raster(tr2) r3 <- raster(tr1, "colMeans")#create a new raster and set all its values to unity. r <- raster(nrows=18, ncols=36) r <- setValues(r,runif(ncell(r),0,1)) #create a Transition object from the raster tr1 <- transition(r,mean,8) #asymmetric asf <- function(x) max(x) - x[1] + x[2] tr2 <- transition(r,asf,8, symm=FALSE) #create RasterLayer objects r1 <- raster(tr1) r2 <- raster(tr2) r3 <- raster(tr1, "colMeans")
Calculates the randomized shortest path distance between points.
rSPDistance(x, from, to, theta, totalNet = "net", method = 1)rSPDistance(x, from, to, theta, totalNet = "net", method = 1)
x |
|
from |
point locations coordinates (of SpatialPoints, matrix or numeric class) |
to |
point locations coordinates (of SpatialPoints, matrix or numeric class) |
theta |
theta is the degree from which the path randomly deviates from the shortest path, 0 < theta < 20 |
totalNet |
total or net movements between cells |
method |
method 1 (as defined in Saerens et al.) or method 2 (a modified version, see below in Details) |
The function implements the algorithm given by Saerens et al. (2009).
Method 1 implements the method as it is. Method 2 uses W = exp(-theta * ln(P)).
distance matrix (S3 class dist or matrix)
Jacob van Etten
Saerens M., L. Yen, F. Fouss, and Y. Achbany. 2009. Randomized shortest-path problems: two related models. Neural Computation, 21(8):2363-2404.
#Create a new raster and set all its values to unity. r <- raster(nrows=18, ncols=36) r <- setValues(r,rep(1,ncell(raster))) #Create a Transition object from the raster tr <- transition(r,mean,4) #Create two sets of coordinates sP1 <- SpatialPoints(cbind(c(65,5,-65),c(55,35,-35))) sP2 <- SpatialPoints(cbind(c(50,15,-40),c(80,20,-5))) #Calculate the RSP distance between the points rSPDistance(tr, sP1, sP2, 1)#Create a new raster and set all its values to unity. r <- raster(nrows=18, ncols=36) r <- setValues(r,rep(1,ncell(raster))) #Create a Transition object from the raster tr <- transition(r,mean,4) #Create two sets of coordinates sP1 <- SpatialPoints(cbind(c(65,5,-65),c(55,35,-35))) sP2 <- SpatialPoints(cbind(c(50,15,-40),c(80,20,-5))) #Calculate the RSP distance between the points rSPDistance(tr, sP1, sP2, 1)
Calculates the shortest path from an origin to a goal
shortestPath(x, origin, goal, ...)shortestPath(x, origin, goal, ...)
x |
|
origin |
SpatialPoints, vector or matrix with coordinates, at the moment only the first cell is taken into account |
goal |
SpatialPoints, vector or matrix with coordinates |
... |
Additional argument: output |
As an additional argument output, the desired output object can be specified: either “TransitionLayer” (default), “TransitionStack” or “SpatialLines”.
If there is more than one path either (1) transition values in the TransitionLayer get values of 1 / number of paths or (2) the SpatialLines object will contain more than one line.
Transition object
Jacob van Etten
#example equivalent to that in the documentation on r.cost/r.drain in GRASS r <- raster(nrows=6, ncols=7, xmn=0, xmx=7, ymn=0, ymx=6, crs="+proj=utm +units=m") r[] <- c(2, 2, 1, 1, 5, 5, 5, 2, 2, 8, 8, 5, 2, 1, 7, 1, 1, 8, 2, 2, 2, 8, 7, 8, 8, 8, 8, 5, 8, 8, 1, 1, 5, 3, 9, 8, 1, 1, 2, 5, 3, 9) tr <- transition(r, function(x) 1/mean(x), 8) # 1/mean: reciprocal to get permeability tr <- geoCorrection(tr) c1 <- c(5.5,1.5) c2 <- c(1.5,5.5) #make a SpatialLines object for visualization sPath1 <- shortestPath(tr, c1, c2, output="SpatialLines") plot(r) lines(sPath1) #make a TransitionLayer for further calculations sPath2 <- shortestPath(tr, c1, c2) plot(raster(sPath2))#example equivalent to that in the documentation on r.cost/r.drain in GRASS r <- raster(nrows=6, ncols=7, xmn=0, xmx=7, ymn=0, ymx=6, crs="+proj=utm +units=m") r[] <- c(2, 2, 1, 1, 5, 5, 5, 2, 2, 8, 8, 5, 2, 1, 7, 1, 1, 8, 2, 2, 2, 8, 7, 8, 8, 8, 8, 5, 8, 8, 1, 1, 5, 3, 9, 8, 1, 1, 2, 5, 3, 9) tr <- transition(r, function(x) 1/mean(x), 8) # 1/mean: reciprocal to get permeability tr <- geoCorrection(tr) c1 <- c(5.5,1.5) c2 <- c(1.5,5.5) #make a SpatialLines object for visualization sPath1 <- shortestPath(tr, c1, c2, output="SpatialLines") plot(r) lines(sPath1) #make a TransitionLayer for further calculations sPath2 <- shortestPath(tr, c1, c2) plot(raster(sPath2))
The following summary methods are available: mean, Median, max, min, range, prod, sum, any, all
## S4 method for signature 'TransitionStack' Summary(x, ..., na.rm = FALSE)## S4 method for signature 'TransitionStack' Summary(x, ..., na.rm = FALSE)
x |
objects |
... |
further arguments passes to or from methods |
na.rm |
logical, should missing values be removed? |
a TransitionLayer
These methods compute a summary statistic based on cell values of layers in a TransitionStack. The result of these methods is always a single TransitionLayer.
Jacob van Etten
#Create a new raster and set all its values to unity. raster <- raster(nrows=18, ncols=36) raster <- setValues(raster,rep(1,ncell(raster))) #Create a Transition object from the raster tr <- transition(raster,mean,4) trS <- stack(tr, tr*2) #Apply a Summary method trSum <- sum(trS) #plot(raster(trMean))#Create a new raster and set all its values to unity. raster <- raster(nrows=18, ncols=36) raster <- setValues(raster,rep(1,ncell(raster))) #Create a Transition object from the raster tr <- transition(raster,mean,4) trS <- stack(tr, tr*2) #Apply a Summary method trSum <- sum(trS) #plot(raster(trMean))
Reciprocal of the sum of the reciprocals of conductance Transition* objects
sumReciprocal(x1, x2)sumReciprocal(x1, x2)
x1 |
|
x2 |
|
To calculate the total resistance of two resistors that are serially
connected, we should add their resistance values. However, if we work
with conductance values, we need to take the reciprocal of the summed
reciprocals of the conductance values. This function does that when
adding two TransitionLayers with conductance values
(matrixValues(tr) == "conductance").
For a TransitionLayer with resistance values
(matrixValues(tr) == "resistance"), the function will not take
reciprocals for that object, but will still take a reciprocal for the
final product (which will consequently have conductance values).
TransitionLayer object containing conductance values.
#Create a new raster and set all its values to unity. raster <- raster(nrows=18, ncols=36) raster <- setValues(raster,rep(1,ncell(raster))) #Create TransitionLayer objects tr1 <- transition(raster,mean,4) tr2 <- tr1 matrixValues(tr1) #Set one to resistance matrixValues(tr2) <- "resistance" #Sum the two objects sumReciprocal(tr1,tr2)#Create a new raster and set all its values to unity. raster <- raster(nrows=18, ncols=36) raster <- setValues(raster,rep(1,ncell(raster))) #Create TransitionLayer objects tr1 <- transition(raster,mean,4) tr2 <- tr1 matrixValues(tr1) #Set one to resistance matrixValues(tr2) <- "resistance" #Sum the two objects sumReciprocal(tr1,tr2)
Create a Transition object from a RasterLayer or RasterBrick object. Transition values are calculated with a user-defined function from the grid values.
transition(x, transitionFunction, directions, ...)transition(x, transitionFunction, directions, ...)
x |
|
transitionFunction |
Function to calculate transition values from grid values |
directions |
Directions in which cells are connected
(4, 8, 16, or other), see |
... |
additional arguments, passed to methods |
Users may use one of three methods to construct a Transition* object with this function.
1) TransitionLayer from RasterLayer
transition(x, transisitonFunction, directions, symm)
When a symmetric transition matrix is required, the user should supply a transitionFunction f that obeys f(i,j) = f(j,i) (a commutative function).
The function transition does no commutativity check.
To obtain an asymmetric transition matrix, a non-commutative function should be supplied and an additional argument ‘symm’ should be set to FALSE.
2) TransitionLayer from RasterBrick
transition(x, transitionFunction = "mahal", directions)
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).
3) TransitionStack from RasterLayer
In contrast with the above methods, this method produces resistance matrices by default.
a) Continuous variables - barriers
transition(x, transitionFunction = "barriers",
directions, symm, intervalBreaks)
This method creates a TransitionStack with each layer
containing a discrete boundary between areas in x.
Areas are defined by intervals in x.
The argument intervalBreaks is a vector of interval
breaks corresponding to the values in x.
If between a pair of cells i and j, min(i,j) < break AND max(i,j) > break,
then the value ij in the transition matrix becomes 1.
All other values in the transition matrix remain 0. The package classInt offers several methods to define intervals. If symm is changed from the default (TRUE) to "up" or "down", it will give either only the upslope (symm = "up") or downslope (symm = "down") barriers.
b) Categorical variables - barriers
transition(x, transitionFunction = "barriers", directions)
In this case, areas are defined as categories in the input raster.
A raster with a categorical variable can be created with
asFactor().
The layers of the resulting TransitionStack contain
all possible combinations of categories.
Which layer contains the combination of categories i and j
out of n categories, can be determined with these formulae:
if symm is TRUE: layer(i,j) = n*(j-1) - j*(j-1)/2 + i-j.
if symm is FALSE and i>j: layer(i,j) = ((n*(j-1) - j*(j-1)/2 + i-j) * 2) - 1.
if symm is FALSE and i<j: layer(i,j) = (n*(j-1) - j*(j-1)/2 + i-j) * 2.
c) Categorical variables - areas
transition(x, transitionFunction = "areas", directions)
Here, areas are also a categorical variable (see under 3b). The layers in the resulting TransitionStack represent each one area. Connections between two cells which are each inside the area are set to 1. Connections between a cell inside and a cell outside the area are set to 0.5. Connections between two cells outside the area are set to 0.
TransitionLayer and TransitionStack (or Transition*) are the core
classes of the package gdistance. They are the main input into the
functions to calculate distances and routes. An object of the class
TransitionLayer contains two main elements:
a. a transition matrix with transition values between connected cells in a raster - an object of class sparseMatrix (package Matrix); b. information on the extent, resolution and projection of the underlying raster - an object of class Raster (package raster).
All slots belong to these two elements from other package, except two additional slots: 1. slot transitionCells, which is only used internally in the package; 2. slot matrixValues indicates if the nonzero values of the transition matrix contains conductance or resistance values.
Class TransitionStack contains one or more transition matrices.
Class Transition is the union of TransitionLayer and TransitionStack.
transitionMatrixObject of class "sparseMatrix"
transitionCellsObject of class "integer"
matrixValuesObject of class "character"
ncolsObject of class "integer"
nrowsObject of class "integer"
crsObject of class "CRS" from sp package
extentObject of class "Extent"
layernameObject of class "vector"
Objects can be created by calls of the form new("Transition", nrows, ncols, xmin, xmax, ymin, ymax, projection).
Class Raster
showClass("TransitionLayer") tr <- new("TransitionLayer", nrows=as.integer(36), ncols=as.integer(18), extent=extent(c(xmin=-180,xmax=180, ymin=-90,ymax=90)), crs=CRS("+proj=longlat +datum=WGS84")) tr <- new("TransitionLayer",nrows=as.integer(36),ncols=as.integer(18), extent=extent(c(xmin=-180, xmax=180, ymin=-90,ymax=90)), crs=CRS(""))showClass("TransitionLayer") tr <- new("TransitionLayer", nrows=as.integer(36), ncols=as.integer(18), extent=extent(c(xmin=-180,xmax=180, ymin=-90,ymax=90)), crs=CRS("+proj=longlat +datum=WGS84")) tr <- new("TransitionLayer",nrows=as.integer(36),ncols=as.integer(18), extent=extent(c(xmin=-180, xmax=180, ymin=-90,ymax=90)), crs=CRS(""))
Methods for functions [ and [<- for object of
the class TransitionLayer.
transitionMatrix(x) <- valuetransitionMatrix(x) <- value
x |
an object of class |
value |
the value to assign |
Methods for functions [[ and [[<- for object of
the class TransitionStack.
Also see adjacencyFromTransition.
#Create a new raster and set all its values to unity. r <- raster(nrows=18, ncols=36) r <- setValues(r,rep(1,ncell(r))) #Create TransitionLayer objects tr1 <- transition(r,mean,4) tr2 <- tr1 #Extracting and replacing tr1[cbind(1:9,1:9)] <- tr2[cbind(1:9,1:9)] tr1[1:9,1:9] <- tr2[1:9,1:9] tr1[1:5,1:5]#Create a new raster and set all its values to unity. r <- raster(nrows=18, ncols=36) r <- setValues(r,rep(1,ncell(r))) #Create TransitionLayer objects tr1 <- transition(r,mean,4) tr2 <- tr1 #Extracting and replacing tr1[cbind(1:9,1:9)] <- tr2[cbind(1:9,1:9)] tr1[1:9,1:9] <- tr2[1:9,1:9] tr1[1:5,1:5]
These functions are to be used to access slots of Transition* objects.
transitionMatrix(x, inflate)transitionMatrix(x, inflate)
x |
object of class |
inflate |
logical (default is |