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.

Internal function.

directionFromEachPoint(from, to = NULL, landscape)

.pointDirection(from, to)

Arguments

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.

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.

landscape

(optional) RasterLayer or SpatRaster. This is only used if to = NULL, in which case all cells are considered to.

Value

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.

Details

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.

See also

distanceFromEachPoint(), which will also return directions if angles = TRUE.

Author

Eliot McIntire

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)