#' Calculate distances and directions between many points and many grid cells
#'
#' This is a modification of [terra::distance()] for the case of many points.
#' This version can often be faster for a single point because it does not return a `RasterLayer`.
#' This is different than [terra::distance()] because it does not take the
#' minimum distance from the set of points to all cells.
#' Rather this returns the every pair-wise point distance.
#' As a result, this can be used for doing inverse distance weightings, seed rain,
#' cumulative effects of distance-based processes etc.
#' If memory limitation is an issue, `maxDistance` will keep memory use down,
#' but with the consequences that there will be a maximum distance returned.
#' This function has the potential to use a lot of memory if there are a lot of
#' `from` and `to` points.
#'
#' This function is cluster aware if the `raster` package is available.
#' If there is a cluster running, it will use it.
#' To start a cluster use `raster::beginCluster()`, with `N` being
#' the number of cores to use. See examples in `SpaDES.core::experiment`.
#'
#' @param from Numeric matrix with 2 or 3 or more columns. They must include x and y,
#' representing x and y coordinates of "from" cell. If there is a column
#' named "id", it will be "id" from `to`, i.e,. specific pair distances.
#' All other columns will be included in the return value of the function.
#'
#' @param to Numeric matrix with 2 or 3 columns (or optionally more, all of which
#' will be returned),
#' x and y, representing x and y coordinates of "to" cells, and
#' optional "id" which will be matched with "id" from `from`. Default is all cells.
#'
#' @param landscape `RasterLayer`. optional. This is only used if `to` is NULL, in which case
#' all cells are considered `to`.
#'
#' @param angles Logical. If `TRUE`, then the function will return angles in radians,
#' as well as distances.
#'
#' @param maxDistance Numeric in units of number of cells. The algorithm will build
#' the whole surface (from `from` to `to`), but will
#' remove all distances that are above this distance.
#' Using this will keep memory use down.
#'
#' @param cumulativeFn A function that can be used to incrementally accumulate
#' values in each `to` location, as the function iterates
#' through each `from`. See Details.
#'
#' @param distFn A function. This can be a function of `landscape`,
#' `fromCell` (single integer value of a from pixel),
#' `toCells` (integer vector value of all the to pixel indices),
#' and `dist`.
#' If `cumulativeFn` is supplied, this will be used to convert
#' the distances to some other set of units that will be accumulated
#' by the `cumulativeFn`. See Details and examples.
#'
#' @param ... Any additional objects needed for `distFn`.
#'
#' @inheritParams splitRaster
#'
#' @return A sorted matrix on `id` with same number of rows as `to`,
#' but with one extra column, `"dists"`, indicating the distance
#' between `from` and `to`.
#'
#' @seealso [rings()], [cir()], [terra::distance()],
#' which can all be made to do the same thing, under specific combinations of arguments.
#' But each has different primary use cases. Each is also faster under different conditions.
#' For instance, if `maxDistance` is relatively small compared to the number of cells
#' in the `landscape`, then [cir()] will likely be faster. If a minimum
#' distance from all cells in the `landscape` to any cell in `from`, then
#' `distanceFromPoints` will be fastest. This function scales best when there are
#' many `to` points or all cells are used `to = NULL` (which is default).
#'
#' @details
#'
#' If the user requires an id (indicating the `from` cell for each `to` cell)
#' to be returned with the function, the user must add an identifier to the
#' `from` matrix, such as `"id"`.
#' Otherwise, the function will only return the coordinates and distances.
#'
#' `distanceFromEachPoint` calls `.pointDistance`, which is not intended to be called
#' directly by the user.
#'
#' This function has the potential to return a very large object, as it is doing pairwise
#' distances (and optionally directions) between from and to.
#' If there are memory limitations because there are many
#' `from` and many `to` points, then `cumulativeFn` and `distFn` can be used.
#' These two functions together will be used iteratively through the `from` points. The
#' `distFn` should be a transformation of distances to be used by the
#' `cumulativeFn` function. For example, if `distFn` is `1 / (1+x)`,
#' the default, and `cumulativeFn` is ``+``, then it will do a sum of
#' inverse distance weights.
#' See examples.
#'
#' @export
#' @importFrom terra ncell xyFromCell
#' @importFrom parallel clusterApply mclapply
#'
#' @example inst/examples/example_distanceFromEachPoint.R
#'
distanceFromEachPoint <- function(from, to = NULL, landscape, angles = NA_real_,
maxDistance = NA_real_, cumulativeFn = NULL,
distFn = function(dist) 1 / (1 + dist), cl, ...) {
matched <- FALSE
fromColNames <- colnames(from)
otherFromCols <- is.na(match(fromColNames, c("x", "y", "id")))
if ("id" %in% fromColNames) {
ids <- unique(from[, "id"])
}
if ("id" %in% colnames(to)) {
matched <- TRUE
}
if (is.null(to)) {
to <- xyFromCell(landscape, 1:ncell(landscape))
}
if (!is.null(cumulativeFn)) {
forms <- names(formals(distFn))
# browser()
fromC <- "fromCell" %in% forms
if (fromC) fromCell <- cellFromXY(landscape, from[, xycolNames])
toC <- "toCells" %in% forms
if (toC) toCells <- cellFromXY(landscape, to[, xycolNames])
land <- "landscape" %in% forms
distFnArgs <- if (land) list(landscape = landscape[]) else list()
if (length(list(...)) > 0) distFnArgs <- append(distFnArgs, list(...))
xDist <- "dist" %in% forms
xDir <- "angle" %in% forms
if (is.character(cumulativeFn)) {
cumulativeFn <- get(cumulativeFn)
}
}
if (!matched) {
nrowFrom <- NROW(from)
if (nrowFrom >= 1) {
if (is.null(cumulativeFn)) {
# if (any(otherFromCols) | isTRUE(angles)) {
out <- lapply(seq_len(nrowFrom), function(k) {
out <- .pointDistance(from = from[k, , drop = FALSE], to = to,
angles = angles, maxDistance = maxDistance,
otherFromCols = otherFromCols)
})
out <- do.call(rbind, out)
# } else {
# maxDistance2 <- if (is.na(maxDistance)) Inf else maxDistance
# browser()
# out <- pointDistance3(fromX = from[, "x"], toX = to[, "x"],
# fromY = from[, "y"], toY = to[, "y"],
# maxDistance = maxDistance2)
# }
} else {
# if there is a cluster, then there are two levels of cumulative function,
# inside each cluster and outside, or "within and between clusters".
# This is the outer one.
# The inner one is the one defined by the user argument.
# outerCumFun <- function(x, from, fromCell, landscape, to, angles, maxDistance, xDir,
# distFnArgs, fromC, toC, xDist, cumulativeFn, distFn, evalEnv) {
#
# fromCell <- eval(fromCell, envir = evalEnv)
# to <- eval(to, envir = evalEnv)
# distFn <- eval(distFn, envir = evalEnv)
#
# fromCell <- rlang::eval_tidy(fromCell)
#
# cumVal <- rep_len(0, NROW(to))
# needAngles <- isTRUE(angles) && isTRUE(xDir)
#
# for (k in seq_len(nrowFrom)) {
# out <- .pointDistance(from = from[k, , drop = FALSE], to = to,
# angles = angles, maxDistance = maxDistance,
# otherFromCols = otherFromCols)
# if (toC)
# toCells <- cellFromXY(landscape, out[, xycolNames])
# if (k == 1) {
# if (fromC) distFnArgs <- append(distFnArgs, list(fromCell = fromCell[k]))
# if (toC) distFnArgs <- append(distFnArgs, list(toCells = toCells))
# if (xDist) distFnArgs <- append(distFnArgs, list(dist = out[, "dists", drop = FALSE]))
# if (needAngles) distFnArgs <- append(distFnArgs, list(angle = out[, "angles", drop = FALSE]))
# } else {
# if (fromC) distFnArgs[["fromCell"]] <- fromCell[k]
# if (toC) distFnArgs[["toCells"]] <- toCells
# if (xDist) distFnArgs[["dist"]] <- out[, "dists"]
# if (needAngles) distFnArgs[["angle"]] <- out[, "angles", drop = FALSE]
# }
# if (!is.null(distFnArgs$landscape))
# distFnArgs$landscape
#
# # call inner cumulative function
# if (isTRUE(!anyNA(maxDistance))) {
# distFnOut <- docall(distFn, args = distFnArgs)
# cumVal[out[, "keptIndex"]] <- docall(
# cumulativeFn, args = list(cumVal[out[, "keptIndex"]], distFnOut)
# )
# } else {
# cumVal <- docall(cumulativeFn, args = list(cumVal, docall(distFn, args = distFnArgs)))
# }
# # cumVal <- docall(cumulativeFn, args = list(cumVal, docall(distFn, args = distFnArgs)))
#
# }
# return(cumVal)
# }
if (missing(cl)) {
cl <- tryCatch(raster::getCluster(), error = function(e) NULL)
on.exit(if (!is.null(cl)) raster::returnCluster(), add = TRUE)
}
outerCumFunArgs <- list(landscape = landscape, to = to, angles = angles,
maxDistance = maxDistance, distFnArgs = distFnArgs,
fromC = fromC, toC = toC, xDist = xDist, xDir = xDir,
cumulativeFn = cumulativeFn, distFn = distFn,
nrowFrom = nrowFrom, otherFromCols = otherFromCols)#, evalEnv = evalEnv)
parFunFun <- function(x) {
# this is a slightly tweaked version of outerCumFun, doing all calls
docall(outerCumFun, append(list(x = x, from = fromList[[x]],
if (fromC) fromCell = fromCellList[[x]]), # nolint
outerCumFunArgs))
}
if (!is.null(cl)) {
if (is.numeric(cl)) {
parFun <- "mclapply"
cl <- seq(cl)
} else {
parFun <- "clusterApply"
}
seqLen <- seq_len(min(nrowFrom, length(cl)))
inds <- rep(seq_along(cl), length.out = nrowFrom)
fromList <- lapply(seqLen, function(ind) {
from[inds == ind, , drop = FALSE]
})
if (fromC) fromCellList <- lapply(seqLen, function(ind) {
fromCell[inds == ind]
})
parFunArgs <- if (is.numeric(cl)) {
list(mc.cores = max(cl), X = seqLen, FUN = parFunFun)
} else {
list(cl = cl, x = seqLen, fun = parFunFun)
}
} else {
parFun <- "lapply"
fromList <- list(from)
if (fromC) fromCellList <- list(fromCell)
parFunArgs <- list(X = 1, FUN = parFunFun)
}
# The actual call
# cumVal <- do.call(get(parFun), args = parFunArgs, quote = TRUE)
cumVal <- docall(get(parFun), args = parFunArgs)
# must cumulativeFn the separate cluster results
if (length(cumVal) >= 1) {
cumVal <- list(Reduce(cumulativeFn, cumVal))
# cumVal[[2]] <- docall(cumulativeFn, cumVal[1:2])
# cumVal[[1]] <- NULL
}
cumVal <- cumVal[[1]]
out <- cbind(to, val = cumVal)
}
} else {
out <- .pointDistance(from = from, to = to, angles = angles,
maxDistance = maxDistance, otherFromCols = otherFromCols)
}
} else {
out <- lapply(ids, function(k) {
.pointDistance(from = from[from[, "id"] == k, , drop = FALSE],
to = to[to[, "id"] == k, , drop = FALSE],
angles = angles, maxDistance = maxDistance,
otherFromCols = otherFromCols)
})
out <- do.call(rbind, out)
}
return(out)
}
#' Alternative point distance (and direction) calculations
#'
#' These have been written with speed in mind.
#'
#' @param from TODO: description needed
#' @param to TODO: description needed
#' @param angles TODO: description needed
#' @param maxDistance TODO: description needed
#' @param otherFromCols TODO: description needed
#'
#' @aliases pointDistance
#' @export
#' @name .pointDistance
#' @rdname distances
.pointDistance <- function(from, to, angles = NA, maxDistance = NA_real_, otherFromCols = FALSE) {
if (!is.na(maxDistance)) {
keep3 <- which(abs(to[, "x"] - from[, "x"]) <= maxDistance)
keep4 <- which(abs(to[keep3, "y"] - from[, "y"]) <= maxDistance)
keep <- keep3[keep4]
# keepOrig <- which((abs(to[, "x"] - from[, "x"]) <= maxDistance) &
# (abs(to[, "y"] - from[, "y"]) <= maxDistance))
# if (!identical(keepOrig, keep)) browser()
to <- to[keep, , drop = FALSE]
}
dists <- cbind(to, dists = sqrt((from[, "x"] - to[, "x"])^2 + (from[, "y"] - to[, "y"])^2))
if (isTRUE(angles)) {
dists <- cbind(dists, angles = .pointDirectionInner(from = from, to = to))
}
if (!is.na(maxDistance)) {
keep2 <- which(dists[, "dists"] <= maxDistance)
dists <- dists[keep2, , drop = FALSE]
}
if (any(otherFromCols)) {
colNums <- seq_len(ncol(dists))
dists <- cbind(dists = dists, from[, otherFromCols])
colnames(dists)[-colNums] <- colnames(from)[otherFromCols]
}
if (!is.na(maxDistance)) {
dists <- cbind(dists, keptIndex = keep[keep2])
}
return(dists)
}
#' Calculate matched point directions
#'
#' Internal function
#'
#' @aliases matchedPointDirection
#' @keywords internal
#' @name .matchedPointDirection
#' @rdname matchedPointDirection
.matchedPointDirection <- function(to, from) {
ids <- unique(from[, "id"])
orig <- order(to[, "id", drop = FALSE], to[, "to", drop = FALSE])
to <- to[orig, , drop = FALSE]
angls <- lapply(ids, function(i) {
m1 <- to[to[, "id"] == i, xycolNames, drop = FALSE]
m2 <- from[from[, "id"] == i, xycolNames, drop = FALSE]
.pointDirection(m2, m1)
})
do.call(rbind, angls)
}
#' Calculate distances and directions between many points and many grid cells
#'
#' This is a modification of [terra::distance()] for the case
#' of many points.
#' This version can often be faster for a single point because it does not return
#' a `RasterLayer`.
#' This is different than [terra::distance()] because it does
#' not take the minimum distance from the set of points to all cells.
#' Rather this returns the every pair-wise point distance.
#' As a result, this can be used for doing inverse distance weightings, seed rain,
#' cumulative effects of distance-based processes etc.
#' If memory limitation is an issue, `maxDistance` will keep memory use down,
#' but with the consequences that there will be a maximum distance returned.
#' This function has the potential to use a lot of memory if there are a lot of
#' `from` and `to` points.
#'
#' `directionFromEachPoint` calls `.pointDirection`, which is
#' not intended to be called directly by the user.
#'
#' If knowing the which from cell matches with which to cell is important,
#' put a column "id" (e.g., starting cell) in the `from` matrix.
#'
#' @param from matrix with 2 or 3 columns, x and y, representing x and y coordinates
#' of `from` cell, and optional `id`, which will be returned,
#' and if `id` column is in `to`, it will be matched with that.
#' @param to matrix with 2 or 3 columns (or optionally more, all of which will be returned),
#' x and y, representing x and y coordinates of `to` cells, and
#' optional `id` which will be matched with `id` from `from`.
#' It makes no sense to have `id` column here with no `id` column
#' in `from`.
#' @param landscape (optional) `RasterLayer` or `SpatRaster`.
#' This is only used if `to = NULL`, in which case all cells are considered `to`.
#'
#' @return A sorted matrix on `id` with same number of rows as `to`,
#' but with one extra column, `angles` indicating the angle in radians between from and to.
#' For speed, this angle will be between `-pi/2` and `3*pi/2`.
#' If the user wants this between say, `0` and `2*pi`,
#' then `angles \%\% (2*pi)` will do the trick. See example.
#'
#' @seealso [distanceFromEachPoint()], which will also return directions if `angles = TRUE`.
#'
#' @examples
#' library(terra)
#'
#' origDTThreads <- data.table::setDTthreads(2L)
#' origNcpus <- options(Ncpus = 2L)
#'
#' N <- 2
#' dirRas <- terra::rast(terra::ext(0,40,0,40), res = 1)
#' coords <- cbind(x = round(runif(N, xmin(dirRas), xmax(dirRas))) + 0.5,
#' y = round(runif(N, xmin(dirRas), xmax(dirRas))) + 0.5,
#' id = 1:N)
#'
#' dirs1 <- directionFromEachPoint(from = coords, landscape = dirRas)
#'
#' # helper function for converting radians to degrees
#' deg2 <- function(radian) (radian * 180)/pi
#' dirs1[, "angles"] <- deg2(dirs1[,"angles"] %% (2*pi))
#' indices <- cellFromXY(dirRas,dirs1[, c("x", "y")])
#' minDir <- tapply(dirs1[, "angles"], indices, function(x) min(x)) # minimum angle
#' dirRas[] <- as.vector(minDir)
#'
#' if (interactive()) {
#' terra::plot(dirRas)
#' start <- terra::vect(coords[, c("x", "y"), drop = FALSE])
#' terra::plot(start, add = TRUE)
#' }
#'
#' # clean up
#' data.table::setDTthreads(origDTThreads)
#' options(Ncpus = origNcpus)
#'
#' @export
#' @importFrom terra ncell xyFromCell
#' @rdname directions
directionFromEachPoint <- function(from, to = NULL, landscape) {
matched <- FALSE
nrowFrom <- NROW(from)
if ("id" %in% colnames(from)) {
ids <- unique(from[, "id"])
} else if (nrowFrom >= 1) {
ids <- seq_len(nrowFrom)
}
if ("id" %in% colnames(to)) {
matched <- TRUE
}
if (is.null(to))
to <- xyFromCell(landscape, 1:ncell(landscape))
if (!matched) {
if (nrowFrom >= 1) {
out <- lapply(seq_len(nrowFrom), function(k) {
out <- .pointDirection(from = from[k, , drop = FALSE], to = to)
cbind(out, id = ids[k])
})
out <- do.call(rbind, out)
} else {
out <- .pointDirection(from = from, to = to)
}
} else {
out <- lapply(ids, function(k) {
.pointDirection(from = from[from[, "id"] == k, , drop = FALSE],
to = to[to[, "id"] == k, , drop = FALSE])
})
out <- do.call(rbind, out)
}
}
#' Calculate the direction from a point to a set of points
#'
#' Internal function.
#'
#' @author Eliot McIntire
#' @keywords internal
#' @rdname directions
.pointDirection <- function(from, to) {
angls <- .pointDirectionInner(from, to)
cbind(to, angles = angls)
}
.pointDirectionInner <- function(from, to) {
rise <- to[, "y"] - from[, "y"]
run <- to[, "x"] - from[, "x"]
pi / 2 - atan2(rise, run) # Convert to geographic 0 = North
}
docall <- function(what, args, quote = FALSE, envir = parent.frame()) {
if (quote) {
args <- lapply(args, enquote)
}
namsArgs <- names(args)
if (is.null(namsArgs) || is.data.frame(args)) {
argn <- args
args <- list()
} else {
hasName <- namsArgs != ""
argn <- lapply(namsArgs[hasName], as.name)
names(argn) <- namsArgs[hasName]
argn <- c(argn, args[namsArgs == ""])
args <- args[hasName]
}
if ("character" %in% class(what)) {
if (is.character(what)) {
fn <- strsplit(what, "[:]{2,3}")[[1]]
what <- if (length(fn) == 1) {
get(fn[[1]], envir = envir, mode = "function")
} else {
get(fn[[2]], envir = asNamespace(fn[[1]]), mode = "function")
}
}
call <- as.call(c(list(what), argn))
} else if ("function" %in% class(what)) {
f_name <- deparse(substitute(what))
call <- as.call(c(list(as.name(f_name)), argn))
args[[f_name]] <- what
} else if ("name" %in% class(what)) {
call <- as.call(c(list(what, argn)))
}
eval(call, envir = args, enclos = envir)
}
outerCumFun <- function(x, from, fromCell, landscape, to, angles, maxDistance, xDir, distFnArgs,
fromC, toC, xDist, cumulativeFn, distFn, nrowFrom, otherFromCols) {
cumVal <- rep_len(0, NROW(to))
needAngles <- isTRUE(angles) && isTRUE(xDir)
nrowFrom <- NROW(from)
for (k in seq_len(nrowFrom)) {
out <- .pointDistance(from = from[k, , drop = FALSE], to = to,
angles = angles, maxDistance = maxDistance,
otherFromCols = otherFromCols)
if (NROW(out) > 0) {
if (toC)
toCells <- cellFromXY(landscape, out[, xycolNames])
if (k == 1) {
if (fromC) distFnArgs <- append(distFnArgs, list(fromCell = fromCell[k]))
if (toC) distFnArgs <- append(distFnArgs, list(toCells = toCells))
if (xDist) distFnArgs <- append(distFnArgs, list(dist = out[, "dists", drop = FALSE]))
if (needAngles) distFnArgs <- append(distFnArgs, list(angle = out[, "angles", drop = FALSE]))
} else {
if (fromC) distFnArgs[["fromCell"]] <- fromCell[k]
if (toC) distFnArgs[["toCells"]] <- toCells
if (xDist) distFnArgs[["dist"]] <- out[, "dists"]
if (needAngles) distFnArgs[["angle"]] <- out[, "angles", drop = FALSE]
}
if (!is.null(distFnArgs$landscape))
distFnArgs$landscape
# call inner cumulative function
if (isTRUE(!anyNA(maxDistance))) {
distFnOut <- docall(distFn, args = distFnArgs)
if (anyNA(distFnOut)) browser()
cumVal[out[, "keptIndex"]] <- cumulativeFn(cumVal[out[, "keptIndex"]], distFnOut)
# cumVal[out[, "keptIndex"]] <- docall(
# cumulativeFn, args = list(cumVal[out[, "keptIndex"]], distFnOut)
# )
} else {
cumVal <- docall(cumulativeFn, args = list(cumVal, docall(distFn, args = distFnArgs)))
}
# cumVal <- docall(cumulativeFn, args = list(cumVal, docall(distFn, args = distFnArgs)))
}
}
return(cumVal)
}
spiralDistances <- function(pixelGroupMap, maxDis, cellSize) {
.requireNamespace("raster")
spiral <- which(raster::focalWeight(pixelGroupMap, maxDis, type = "circle") > 0, arr.ind = TRUE) -
ceiling(maxDis/cellSize) - 1
spiral <- cbind(spiral, dists = sqrt( (0 - spiral[, 1]) ^ 2 + (0 - spiral[, 2]) ^ 2))
spiral <- spiral[order(spiral[, "dists"], apply(abs(spiral), 1, sum),
abs(spiral[, 1]), abs(spiral[, 2])), NULL, drop = FALSE]
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.