/
shortestPath.R
146 lines (132 loc) · 4.26 KB
/
shortestPath.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
#' Shortest path
#'
#' Calculates the shortest path from an origin to a goal
#'
#' @name shortestPath
#' @aliases shortestPath
#' @aliases shortestPath,TransitionLayer,Coords,Coords-method
#'
#' @param x \code{TransitionLayer} object
#' @param origin SpatialPoints, vector or matrix with coordinates,
#' at the moment only the first cell is taken into account
#' @param goal SpatialPoints, vector or matrix with coordinates
#' @param ... Additional argument: output
#' @return Transition object
#' @details
#' As an additional argument output, the desired output object can be specified:
#' either \dQuote{TransitionLayer} (default),
#' \dQuote{TransitionStack} or \dQuote{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.
#' @author Jacob van Etten
#' @seealso \code{\link{costDistance}}, \code{\link{accCost}}
#' @examples
#' #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))
#'
#' @exportMethod shortestPath
setGeneric("shortestPath", function(x, origin, goal, ...) {
standardGeneric("shortestPath")
}
)
#check if Transition and RasterLayers coincide, etc.
#' @export
setMethod("shortestPath",
signature(x = "TransitionLayer",
origin = "Coords",
goal = "Coords"),
def = function(x, origin, goal,
output="TransitionLayer")
{
origin <- .coordsToMatrix(origin)
goal <- .coordsToMatrix(goal)
return(.shortestPath(x, origin, goal, output))
}
)
.shortestPath <- function(x, origin, goal, output)
{
originCells <- raster::cellFromXY(x, origin)
goalCells <- raster::cellFromXY(x, goal)
indexOrigin <- originCells
indexGoal <- goalCells
y <- transitionMatrix(x)
if(isSymmetric(y)) {
mode <- "undirected"
}else{
mode <- "directed"
}
adjacencyGraph <- igraph::graph.adjacency(y, mode=mode, weighted=TRUE)
E(adjacencyGraph)$weight <- 1/E(adjacencyGraph)$weight
shortestPaths <- get.shortest.paths(adjacencyGraph,
indexOrigin, indexGoal)$vpath
if(output=="TransitionLayer")
{
result <- x
transitionMatrix(result) <- Matrix::Matrix(0, ncol=ncell(x), nrow=ncell(x))
for(i in 1:length(shortestPaths))
{
sPVector <- shortestPaths[[i]]
adj <- cbind(sPVector[-(length(sPVector))], sPVector[-1])
adj <- rbind(adj,cbind(adj[,2], adj[,1]))
transitionMatrix(result)[adj] <- 1/length(shortestPaths) + transitionMatrix(result)[adj]
}
}
if(output=="TransitionStack")
{
result <- x
transitionMatrix(result) <- Matrix::Matrix(0, ncol=ncell(x), nrow=ncell(x))
for(i in 1:length(shortestPaths))
{
resultNew <- result
sPVector <- shortestPaths[[i]]
adj <- cbind(sPVector[-(length(sPVector))], sPVector[-1])
adj <- rbind(adj,cbind(adj[,2], adj[,1]))
transitionMatrix(resultNew)[adj] <- 1/length(shortestPaths)
result <- raster::stack(result, resultNew)
}
result <- result[[2:nlayers(result)]]
}
if(output=="SpatialLines")
{
linesList <- vector(mode="list", length=length(shortestPaths))
for(i in 1:length(shortestPaths))
{
sPVector <- shortestPaths[[i]]
coords <- raster::xyFromCell(x, sPVector)
linesList[[i]] <- Line(coords)
}
# Suggested by Sergei Petrov
LinesObject <- mapply(Lines,
slinelist = linesList,
ID = as.character(1:length(shortestPaths)),
SIMPLIFY = FALSE)
result <- sp::SpatialLines(LinesObject, proj4string = sp::CRS(projection(x)))
}
return(result)
}