lcp_over_surface.Rd

% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/lcps.R
\name{lcp_over_surface}
\alias{lcp_over_surface}
\title{Calculate shortest path(s) and/or distance(s) over a surface between origin and destination coordinates}
\usage{
lcp_over_surface(
  origin,
  destination,
  surface,
  crop = NULL,
  buffer = NULL,
  aggregate = NULL,
  mask = NULL,
  mask_inside = TRUE,
  plot = TRUE,
  goal = 1,
  combination = "pair",
  method = "cppRouting",
  cppRouting_algorithm = "bi",
  cl = NULL,
  varlist = NULL,
  use_all_cores = FALSE,
  check = TRUE,
  verbose = TRUE
)
}
\arguments{
\item{origin}{A matrix which defines the coordinates (x, y) of the starting location(s). Coordinates should lie on a plane (i.e., Universal Transverse Mercator projection).}

\item{destination}{A matrix which defines the coordinates (x, y) of the finishing location(s). Coordinates should lie on a plane (i.e., Universal Transverse Mercator projection).}

\item{surface}{A \code{\link[raster]{raster}} over which the object (e.g., individual) must move from \code{origin} to \code{destination}. The \code{surface} must be planar (i.e., Universal Transverse Mercator projection) with units of metres in x, y and z directions (m). The \code{surface}'s \code{\link[raster]{resolution}} is taken to define the distance between horizontally and vertically connected cells and must be the same in both x and y directions (for \code{surface}'s with unequal horizontal resolution, \code{\link[raster]{resample}} can be used to equalise resolution: see Examples). Any cells with NA values (e.g., due to missing data) are treated as `impossible' to move though by the algorithm. In this case, the \code{surface} might need to be pre-processed so that NAs are replaced/removed before implementing the function, depending on their source.}

\item{crop}{(optional) An \code{\link[raster]{extent}} object that is used to \code{\link[raster]{crop}} the extent of the \code{surface}, before the least-cost algorithms are implemented. This may be useful for large rasters to reduce memory requirements/computation time.}

\item{buffer}{(optional) A named list of arguments, passed to \code{\link[rgeos]{gBuffer}} (e.g. \code{buffer = list(width = 1000)}) (m) that is used to define a buffer around a Euclidean transect connecting the \code{origin} and \code{destination}. (This option can only be implemented for a single \code{origin} and \code{destination} pair.) The \code{surface} is then cropped to the extent of this buffer before the least-cost algorithms are implemented. This may be useful for large rasters to reduce memory requirements and/or computation time.}

\item{aggregate}{(optional) A named list of arguments, passed to \code{\link[raster]{aggregate}}, to aggregate raster cells before the least-cost algorithms are implemented. This may be useful for large rasters to reduce memory requirements and/or computation time.}

\item{mask}{(optional) A Raster or Spatial \code{\link[raster]{mask}} that is used to prevent movement over `impossible' areas on the \code{surface}. This must also lie on a planar surface (i.e., Universal Transverse Mercator projection). For example, for marine animals, \code{mask} might be a \code{\link[sp]{SpatialPolygonsDataFrame}} which defines the coastline. The effect of the \code{mask} depends on \code{mask_inside} (see below).}

\item{mask_inside}{A logical input that defines whether or not to mask the \code{surface} inside (\code{TRUE}) or outside (\code{FALSE}) of the \code{mask} (see \code{\link[flapper]{mask_io}}).}

\item{plot}{A logical input that defines whether or not to plot the inputted and processed surfaces. If \code{TRUE}, the inputted and processed plots are produced side-by-side. For the inputted surface, the \code{mask} and the region selected (via \code{crop} and/or \code{buffer}) are shown along with the \code{origin} and \code{destination}. For the processed surface, the surface and the \code{origin} and \code{destination} are shown, along with the shortest path(s) (if and once computed: see \code{goal}). This is useful for checking that any \code{surface} processing steps have been applied correctly and the \code{origin} and \code{destination} are positioned correctly on the \code{surface}.}

\item{goal}{An integer that defines the output of the function: \code{goal = 1} computes shortest distances, \code{goal = 2} computes shortest paths and \code{goal = 3} computes both shortest paths and the corresponding distances. Note that \code{goal = 3} results in least-cost algorithms being implemented twice, which will be inefficient for large problems; in this case, use \code{goal = 2} to compute shortest paths and then calculate their distance using outputs returned by the function (see Value).}

\item{combination}{A character string (\code{"pair"} or \code{"matrix"}) that defines whether or not to compute shortest distances/paths for (a) each sequential \code{origin} and \code{destination} pair of coordinates (\code{combination = "pair"}) or (b) all combinations of \code{origin} and \code{destination} coordinates (\code{combination = "matrix"}). This argument is only applicable if there is more than one \code{origin} and \code{destination}. For \code{combination = "pair"}, the number of \code{origin} and \code{destination} coordinates needs to be the same, since each \code{origin} is matched with each \code{destination}.}

\item{method}{A character string (\code{"cppRouting"} or \code{"gdistance"}) that defines the method used to compute the shortest distances between the \code{origin} and the \code{destination}. \code{"cppRouting"} is the default \code{method}. Under this option, functions in the \code{cppRouting} package are used to compute the shortest paths (\code{\link[cppRouting]{get_path_pair}} or \code{\link[cppRouting]{get_multi_paths}} for each pair of coordinates or for all combinations of coordinates, respectively) and/or distances (\code{\link[cppRouting]{get_distance_pair}} or \code{\link[cppRouting]{get_distance_matrix}}). This package implements functions written in C++ and massively outperforms the other \code{method = "gdistance"} for large problems. Otherwise, if \code{method = "gdistance"}, functions in the \code{\link[gdistance]{gdistance}} are called iteratively to compute shortest paths (via \code{\link[gdistance]{shortestPath}}) or distances (via \code{\link[gdistance]{costDistance}}).}

\item{cppRouting_algorithm}{A character string that defines the algorithm used to compute shortest paths or distances. This is only applicable if \code{method = "cppRouting"}: \code{method = "gdistance"} implements Dijkstra's algorithm only. For shortest paths or their distances between pairs of coordinates, the options are \code{"Dijkstra"}, \code{"bi"}, \code{"A*"} or \code{"NBA"} for the uni-directional Dijkstra, bi-directional Dijkstra, A star unidirectional search or new bi-directional A star algorithms respectively (see \code{\link[cppRouting]{get_path_pair}} or \code{\link[cppRouting]{get_distance_pair}}). For shortest paths between all combinations of coordinates, \code{cppRouting_algorithm} is ignored and the Dijkstra algorithm is implemented recursively. For shortest distances between all combinations of coordinates, the options are \code{"phast"} or \code{"mch"} (see \code{\link[cppRouting]{get_distance_matrix}}).}

\item{use_all_cores, cl, varlist}{Parallelisation options for \code{method = "cppRouting"} (\code{use_all_cores}) or \code{method = "gdistance"} (\code{cl} and \code{varlist}) respectively. If \code{method = "cppRouting"}, parallelisation is implemented via \code{use_all_cores} for computing shortest distances only (not computing shortest paths). \code{use_all_cores} is a logical input that defines whether or not to use all cores for computing shortest distance(s). If \code{method = "gdistance"}, parallelisation is implemented via \code{cl} and \code{varlist} for both shortest paths and distances function calls. \code{cl} is (a) a cluster object from \code{\link[parallel]{makeCluster}} or (b) an integer that defines the number of child processes. \code{varlist} is a character vector of variables for export (see \code{\link[flapper]{cl_export}}). Exported variables must be located in the global environment. If a cluster is supplied, the connection to the cluster is closed within the function (see \code{\link[flapper]{cl_stop}}). For further information, see \code{\link[flapper]{cl_lapply}} and \code{\link[flapper]{flapper-tips-parallel}}.}

\item{check}{A logical input that defines whether or not to check function inputs. If \code{TRUE}, internal checks are implemented to check user-inputs and whether or not inputted coordinates are in appropriate places on the processed \code{surface} (for instance, to ensure inputted coordinates do not lie over masked areas). This helps to prevent intractable error messages. If \code{FALSE}, these checks are not implemented, so function progress may be faster initially (especially for large \code{origin}/\code{destination} coordinate matrices).}

\item{verbose}{A logical input that defines whether or not to print messages to the console to monitor function progress. This is especially useful with a large \code{surface} since the algorithms are computationally intensive.}
}
\value{
\subsection{A named list}{The function returns a named list. The most important element(s) of this list are `path_lcp' and/or `dist_lcp', the shortest path(s) and/or distance(s) (m) between \code{origin} and \code{destination} coordinate pairs/combinations. `path_lcp' is returned if \code{goal = 2} or \code{goal = 3} and `dist_lcp' is returned if \code{goal = 1} or \code{goal = 3}. `path_lcp' contains (a) a dataframe with the cells comprising each path (`cells'), (b) a named list containing a \code{\link[sp]{SpatialLines}} object for each path (`SpatialLines') and (c) a named list of matrices of the coordinates of each path (`coordinates'). `dist_lcp' is a (a) numeric vector or (b) matrix with the distances (m) between each pair or combination of coordinates respectively. If `dist_lcp' is computed, `dist_euclid', the Euclidean distances (m) between the \code{origin} and \code{destination}, is also returned for comparison.
}

\subsection{Common elements}{Other elements of the list record important outputs at sequential stages of the algorithm's progression. These include the following elements: `args', a named list of user inputs; `time', a dataframe that defines the times of sequential stages in the algorithm's progression; `surface', the surface over which shortest distances are computed (this may differ from the inputted surface if any of the processing options, such as \code{crop}, have been implemented); `surface_param', a named list that defines the cell IDs, the number of rows, the number of columns, the coordinates of the implemented surface and the cell IDs of the \code{origin} and \code{destination} nodes; `cost', a named list of arguments that defines the distances (m) between connected cells under a rook's or bishop's movement (`dist_rook' and `dist_bishop'), the planar and vertical distances between connected cells (`dist_planar' and `dist_vertical') and the total distance between connected cells (`dist_total'); and `cppRouting_param' or `gdistance_param', a named list of arguments used to compute shortest paths/distances via \code{cppRouting} or \code{\link[gdistance]{gdistance}} (see below).
}

\subsection{Method-specific elements}{If \code{method = "cppRouting"}, the `cppRouting_param' list contains a named list of arguments passed to \code{\link[cppRouting]{makegraph}} (`makegraph_param') as well as \code{\link[cppRouting]{get_path_pair}} (`get_path_pair_param') or \code{\link[cppRouting]{get_multi_paths}} (`get_multi_paths_param') and/or \code{\link[cppRouting]{get_distance_pair}} (`get_distance_pair_param') or \code{\link[cppRouting]{get_distance_matrix}} (`get_distance_matrix_param'), depending on whether or not shortest paths and/or distances have been computed (see \code{goal}) and whether or not shortest paths/distances have been computed for each pair of coordinates or all combinations of coordinates. If \code{method = "gdistance"}, this list contains a named list of arguments passed iteratively, for each pair/combination of coordinates, to \code{\link[gdistance]{shortestPath}} (`shortestPath_param') or \code{\link[gdistance]{costDistance}} (`costDistance_param'). This includes an object of class TransitionLayer (see \code{\link[gdistance]{Transition-classes}}), in which the \code{transitionMatrix} slot contains a (sparse) matrix that defines the ease of moving between connected cells (the reciprocal of the `dist_total' matrix).
}

\subsection{Plot}{If \code{plot = TRUE}, a plot is also produced of the inputted and processed surfaces that are used in the calculations, along with the shortest path(s) (if and once computed).
}
}
\description{
This function computes the shortest path(s) and/or distance(s) over a \code{surface} between \code{origin} and \code{destination} coordinates. To implement this function, \code{origin} and \code{destination} coordinates need to be specified as matrices and the surface over which movement occurs should be supplied as a \code{\link[raster]{raster}}. Since determining shortest paths can be computationally and memory-intensive, the \code{surface} can be reduced in size and/or resolution before these are computed, by (a) cropping the surface within user-defined extents; (b) focusing on a buffer zone along a Euclidean transect connecting \code{origin} and \code{destination} coordinates; (c) aggregating the surface to reduce the resolution; and/or (d) masking out areas over which movement is impossible (e.g., land for marine animals). Then, the function computes distances between connected cells, given (a) the planar distances between connected cells and (b) their difference in elevation. These distances are taken as a measure of `cost'. For each pair of \code{origin} and \code{destination} coordinates, or for all combinations of coordinates, these distances are used to compute the least-cost (i.e., shortest) path and/or the distance of this path, using functions in the \code{cppRouting} or \code{\link[gdistance]{gdistance}} package. The function returns the shortest path(s) and/or their distance(s) (m) along with a plot and a list of objects involved in the calculations.
}
\details{
\subsection{Methods}{
This function was motivated by the need to determine the shortest paths and their distances between points for benthic animals, which must move over the seabed to navigate from A to B. For these animals, especially in areas with heterogeneous bathymetric landscapes and/or coastline, the shortest path that an individual must travel to move from A and B may differ substantially from the Euclidean path that is often used as a proxy for distance in biological studies. However, this function can still be used in situations where the surface over which an individual must move is irrelevant (e.g., for a pelagic animal), by supplying a flat surface; then shortest paths/distances simply depend on the planar distances between locations and any barriers (e.g., the coastline). (However, this process will be somewhat inefficient.)

The function conceptualises an object moving across a landscape as a queen on a chessboard which can move, in eight directions around its current position, across this surface. Given the potentially large number of possible paths between an \code{origin} and \code{destination}, the surface may be reduced in extent or size before the game begins. To determine shortest path/distance over the surface between each \code{origin} and \code{destination} pair/combination, the function first considers the distance that an object must travel between pairs of connected cells. This depends on the planar distances between cells and their differences in elevation. Planar distances (\eqn{d_p}, m) depend on the movement type: under a rook's movement (i.e., horizontally or vertically), the distance (\eqn{d_{p,r}}) between connected cells is extracted from the raster's resolution (which is assumed to be identical in the x and y directions); under a bishop's movement (i.e., diagonally), the distance between connected cells \eqn{d_{p,b}} is given by Pythagoras' Theorem: \eqn{d_{p,b} = \sqrt{(d_{p, r}^2 + d_{p, r}^2)}}. Vertical distances (\eqn{d_v}, m) are simply the differences in height between cells. The total distance (\eqn{d_t}) between any two connected cells is a combination of these distances given by Pythagoras' Theorem: \eqn{d_t = \sqrt{(d_p^2 + d_v^2)}}. These distances are taken to define the `cost' of movement between connected cells. Thus, `costs' are symmetric (i.e., the cost of moving from A to B equals the cost of moving from B to A).

This cost surface is then used to compute the shortest path and/or distance of the shortest path between each \code{origin} and \code{destination} pair/combination using functions in the \code{cppRouting} or \code{\link[gdistance]{gdistance}} package. The functions implemented depend on the \code{goal} (i.e., whether the aim is to compute shortest paths, shortest distances or both) and, if there is more than one \code{origin}/\code{destination}, the \code{combination} type (i.e., whether to compute shortest paths/distances for each sequential pair of coordinates or all possible combinations of coordinates).
}

\subsection{Warnings}{
The function returns a warning produced by \code{\link[gdistance]{transition}} which is implemented to facilitate the definition of the cost surface, before shortest paths/distances are computed by either method: `In .TfromR(x, transitionFunction, directions, symm) : transition function gives negative values'. This warning arises because the height differences between connecting cells can be negative. It can be safely ignored.
}
}
\examples{
#### Example types
# Shortest distances between a single origin and a single destination
# Shortest paths between a single origin and a single destination
# Shortest distances/paths between origin/destination pairs
# Shortest distances/paths between all origin/destination combinations

#### Simulate a hypothetical landscape
# Define a miniature, blank landscape with appropriate dimensions
proj_utm <- sp::CRS(SRS_string = "EPSG:32629")
r <- raster::raster(
  nrows = 3, ncols = 3,
  crs = proj_utm,
  resolution = c(5, 5),
  ext = raster::extent(0, 15, 0, 15)
)
# Define a matrix of hypothetical values for the landscape
mat <- matrix(c(
  5, 10, 3,
  2, 1, 4,
  5, 6, 6
), ncol = 3, nrow = 3, byrow = TRUE)
r[] <- mat
# Visualise simulated landscape
raster::plot(r)
raster::text(r)
# Extract coordinates of cells
rxy <- raster::coordinates(r)


############################################################################
#### Shortest distances between a single origin and a single destination

#### Example (1): Find the distance between a single origin and destination
# ... using the "cppRouting" method:
out1 <- lcp_over_surface(
  origin = rxy[1, , drop = FALSE],
  destination = rxy[2, , drop = FALSE],
  surface = r
)
# Extract shortest distance
out1$dist_lcp

#### Example (1) continued: An explanation of function outputs
# The function returns a list:
# The 'args' element simply contains all inputted arguments
out1$args
# The 'surface' element contains the surface used to compute shortest distances
# ... This may differ from $args$surface if cropped, buffered etc.
out1$surface
# The 'surface_param' element contains the cell IDs, number of rows, cells and coordinates
# ... of this surface
out1$surface_param
# The 'cost' element is a list of objects that define the cost matrix:
# ... 'dist_rook' and 'dist_bishop' are matrices which define the distance of planar
# ... ... movement from one cell to any other cell under a rook's or bishop's movement.
# ... ... For example, the planar distance of moving from cell 1 to cell 2 to is 5 m:
out1$cost$dist_rook
out1$cost$dist_bishop
# ... 'dist_planar' gives the planar distance between connected cell combinations
# ... ... under a queen's movements:
out1$cost$dist_planar
# ... 'dist_vertical' gives the vertical distance between connected cells
out1$cost$dist_vertical
# ... and 'dist_total' gives the total distance between connected cells
out1$cost$dist_total
# 'dist_euclid' is the Euclidean distance between the origin and destination
out1$dist_euclid
# 'cppRouting_param' contains lists of parameters passed to (a) makegraph() and (b)
# ... get_distance_matrix() to compute shortest distances using cppRouting
utils::str(out1$cppRouting_param)
# 'time' records the time of each stage
out1$time

#### Example (2): Find the distance between a single origin and destination
# ... using the "gdistance" method
out2 <- lcp_over_surface(
  origin = rxy[1, , drop = FALSE],
  destination = rxy[2, , drop = FALSE],
  surface = r,
  method = "gdistance"
)
# Extract distance
out2$dist_lcp
# Elements of the returned list are the same apart from 'gdistance_param'
# ... which contains a list of arguments passed to costDistance() to compute
# ... shortest distances
utils::str(out2$gdistance_param)

#### Example (3): Find the distances between other origins and destinations
## Implement function to determine shortest distances:
# shortest distance between cell 1 and 6 via gdistance
lcp_over_surface(
  origin = rxy[1, , drop = FALSE],
  destination = rxy[6, , drop = FALSE],
  surface = r,
  method = "gdistance",
  verbose = FALSE
)$dist_lcp
# shortest distance between cell 1 and 6 via cppRouting
lcp_over_surface(
  origin = rxy[1, , drop = FALSE],
  destination = rxy[6, , drop = FALSE],
  surface = r,
  method = "cppRouting",
  verbose = FALSE
)$dist_lcp
## Compare to manually computed distances
# The shortest distance from cell 1 to cell 6 is to move diagonally
# ... from cell 1 to 5 and then 6.
# Define planar distances of moving like a rook (d_pr) or bishop (d_pb)
d_pr <- 5
d_pb <- sqrt(5^2 + 5^2)
# Define total distance travelled along shortest path, as computed by the
# ... algorithm to demonstrate we obtain the same value:
sqrt((5 - 1)^2 + d_pb^2) + sqrt((4 - 1)^2 + d_pr^2)

#### Example (4): Find the shortest distances around NAs
## Force the 5th cell to be NA
rtmp <- r
rtmp[5] <- NA
raster::plot(rtmp)
raster::text(rtmp)
## Compute shortest distances via algorithm:
# Now compute shortest distances, which we can see have increased the distance
# ... since movement through an NA cell is not allowed. Therefore, if this NA
# ... reflects missing data, it may be appropriate to interpolate NAs using
# ... surrounding cells (e.g., see raster::approxNA()) so that movement
# ... is possible through these cells
out1 <- lcp_over_surface(
  origin = rxy[1, , drop = FALSE],
  destination = rxy[6, , drop = FALSE],
  surface = rtmp,
  method = "gdistance",
  verbose = FALSE
)
out2 <- lcp_over_surface(
  origin = rxy[1, , drop = FALSE],
  destination = rxy[6, , drop = FALSE],
  surface = rtmp,
  method = "cppRouting",
  verbose = FALSE
)
out1$dist_lcp
out2$dist_lcp
## Compare to manual calculations of shortest route:
# Route option one: cell 1 to 2, 2 to 6
sqrt((5 - 10)^2 + d_pr^2) + sqrt((10 - 4)^2 + d_pb^2)
# Or, using the numbers computed in the dist_total object:
out1$cost$dist_total[1, 2] + out1$cost$dist_total[2, 6]
## Compare to effect of making a value in the landscape extremely large
# In the same way, we can force the shortest path away from particular areas
# ... by making the height of the landscape in those areas very large or Inf:
rtmp[5] <- 1e20
raster::plot(rtmp)
raster::text(rtmp)
lcp_over_surface(
  origin = rxy[1, , drop = FALSE],
  destination = rxy[6, , drop = FALSE],
  surface = rtmp,
  method = "cppRouting",
  verbose = FALSE
)$dist_lcp
lcp_over_surface(
  origin = rxy[1, , drop = FALSE],
  destination = rxy[6, , drop = FALSE],
  surface = rtmp,
  method = "gdistance",
  verbose = FALSE
)$dist_lcp

#### Example (4): Find the distances between points on a real landscape
## We will use some example bathymetry data:
dat_gebco_oban <- prettyGraphics::dat_gebco
raster::plot(dat_gebco_oban)
## Process bathymetry data before function implementation
# (a) Define utm coordinates:
dat_gebco_utm <- raster::projectRaster(dat_gebco_oban, crs = proj_utm)
raster::res(dat_gebco_utm)
# (b) Resample so that the resolution in the x and y directions is identical
dat_gebco_utm_planar <- raster::raster(
  crs = proj_utm,
  ext = raster::extent(dat_gebco_utm),
  resolution = 250
)
dat_gebco_utm_planar <- raster::resample(dat_gebco_utm, dat_gebco_utm_planar, method = "bilinear")
# Examine processed raster
pp <- par(mfrow = c(1, 2))
raster::plot(dat_gebco_utm, main = "UTM raster")
raster::plot(dat_gebco_utm_planar, main = "UTM raster with equal res")
par(pp)
## Define example origin and destination
set.seed(1)
dat_gebco_utm_planar_xy <- raster::coordinates(dat_gebco_utm_planar)
index <- sample(1:nrow(dat_gebco_utm_planar_xy), 2)
origin <- dat_gebco_utm_planar_xy[index[1], , drop = FALSE]
destination <- dat_gebco_utm_planar_xy[index[2], , drop = FALSE]
## Implement function to compute shortest distances
out_gebco1 <- lcp_over_surface(
  origin = origin,
  destination = destination,
  surface = dat_gebco_utm_planar,
  method = "gdistance"
)
out_gebco2 <- lcp_over_surface(
  origin = origin,
  destination = destination,
  surface = dat_gebco_utm_planar,
  method = "cppRouting"
)
# Compare Euclidean and shortest distances
out_gebco1$dist_euclid
out_gebco1$dist_lcp
out_gebco2$dist_euclid
out_gebco2$dist_lcp

#### Example (5A): Reduce the complexity of the landscape by cropping
ext <- raster::extent(dat_gebco_utm_planar)
ext[3] <- 6255000
ext[4] <- 6265000
out_gebco1 <- lcp_over_surface(
  origin = origin,
  destination = destination,
  surface = dat_gebco_utm_planar,
  method = "gdistance",
  crop = ext
)
out_gebco1$dist_lcp

#### Example (5B): Reduce the complexity of the landscape around a buffer
out_gebco1 <- lcp_over_surface(
  origin = origin,
  destination = destination,
  surface = dat_gebco_utm_planar,
  method = "gdistance",
  buffer = list(width = 1000)
)
out_gebco1$dist_lcp

#### Example (5C): Reduce the complexity of the landscape via aggregation
out_gebco1 <- lcp_over_surface(
  origin = origin,
  destination = destination,
  surface = dat_gebco_utm_planar,
  method = "gdistance",
  buffer = list(width = 1000),
  aggregate = list(fact = 5, fun = mean, na.rm = TRUE)
)
out_gebco1$dist_lcp

#### Example (6C): Implement a mask
# Define coastline
dat_coast_around_oban <- prettyGraphics::dat_coast_around_oban
coastline <- sp::spTransform(dat_coast_around_oban, proj_utm)
# Visualise bathymetry and coastline
raster::plot(dat_gebco_utm_planar)
raster::lines(coastline)
# Define example origin and destination within the sea
origin <- matrix(c(714000, 6260000), ncol = 2)
destination <- matrix(c(721000, 6265000), ncol = 2)
# Implement algorithm
out_gebco1 <- lcp_over_surface(
  origin = origin,
  destination = destination,
  surface = dat_gebco_utm_planar,
  method = "gdistance",
  crop = raster::extent(coastline),
  mask = coastline,
  mask_inside = TRUE
)
out_gebco2 <- lcp_over_surface(
  origin = origin,
  destination = destination,
  surface = dat_gebco_utm_planar,
  method = "cppRouting",
  crop = raster::extent(coastline),
  mask = coastline,
  mask_inside = TRUE
)
# Compare Euclidean and least-cost distances
out_gebco1$dist_euclid
out_gebco1$dist_lcp
out_gebco2$dist_euclid
out_gebco2$dist_lcp

#### Example (7) Implement shortest distance algorithms in parallel:
# With the default method ("cppRouting"), use use_all_cores = TRUE
out_gebco1 <- lcp_over_surface(
  origin = origin,
  destination = destination,
  surface = dat_gebco_utm_planar,
  method = "cppRouting",
  crop = raster::extent(coastline),
  mask = coastline,
  mask_inside = TRUE,
  use_all_cores = TRUE
)
# With method = "gdistance" use cl argument
out_gebco2 <- lcp_over_surface(
  origin = origin,
  destination = destination,
  surface = dat_gebco_utm_planar,
  method = "cppRouting",
  crop = raster::extent(coastline),
  mask = coastline,
  mask_inside = TRUE,
  cl = parallel::makeCluster(2L)
)
out_gebco1$dist_lcp
out_gebco2$dist_lcp


############################################################################
#### Shortest paths between a single origin and a single destination

#### Example (8) Shortest paths (goal = 2) only using default method
# Implement function
out1 <- lcp_over_surface(
  origin = rxy[1, , drop = FALSE],
  destination = rxy[6, , drop = FALSE],
  surface = r,
  goal = 2
)
# The path is stored in path_lcp, This includes:
# ... (a) a dataframe with the cells comprising each path:
# ... (b) a SpatialLines object of the path
# ... (c) a matrix of the coordinates of the path
out1$path_lcp
# For method = "cppRouting", paths between pairs of coordinates are computed
# ... by cppRouting::get_path_pair(), the arguments of which are retained in
# ... this list:
out1$cppRouting_param$get_path_pair_param
# Note the path is also added to the plot produced, if plot = TRUE.

#### Example (9) Shortest distances and paths (goal = 3) using default method
out1 <- lcp_over_surface(
  origin = rxy[1, , drop = FALSE],
  destination = rxy[6, , drop = FALSE],
  surface = r,
  goal = 3
)
out1$dist_lcp
out1$path_lcp

#### Example (10) Shortest distances and paths (goal = 3) via gdistance
out1 <- lcp_over_surface(
  origin = rxy[1, , drop = FALSE],
  destination = rxy[6, , drop = FALSE],
  surface = r,
  goal = 3,
  method = "gdistance"
)
# As above, paths are stored in path_lcp and distances in dist_lcp
out1$path_lcp
# For method = "gdistance", paths between pairs of coordinates are computed
# ... by repeated calls to gdistance::shortestPath(), for each pair of coordinates
# ... in the following list of arguments:
out1$gdistance_param$shortestPath_param

#### Example (11): Parallelisation proceeds as described above via
# ... use_all_cores or cl arguments. For cppRouting, parallelisation
# ... is only implemented for distance calculations (so not if goal = 2),
# ... while parallelisation is implemented for both distance and shortest
# ... paths for method = "gdistance"
out1 <- lcp_over_surface(
  origin = rxy[1, , drop = FALSE],
  destination = rxy[6, , drop = FALSE],
  surface = r,
  goal = 3,
  use_all_cores = TRUE
)
out1 <- lcp_over_surface(
  origin = rxy[1, , drop = FALSE],
  destination = rxy[6, , drop = FALSE],
  surface = r,
  goal = 2,
  method = "gdistance",
  cl = parallel::makeCluster(2L)
)
out1 <- lcp_over_surface(
  origin = rxy[1, , drop = FALSE],
  destination = rxy[6, , drop = FALSE],
  surface = r,
  goal = 3,
  method = "gdistance",
  cl = parallel::makeCluster(2L)
)


############################################################################
#### Shortest distances/paths between origin/destination pairs

#### Example (12): Shortest distances/paths computed in sequence:
# cppRouting method
out1 <- lcp_over_surface(
  origin = rxy[1:2, , drop = FALSE],
  destination = rxy[5:6, , drop = FALSE],
  surface = r,
  goal = 3
)
# gdistance method
out2 <- lcp_over_surface(
  origin = rxy[1:2, , drop = FALSE],
  destination = rxy[5:6, , drop = FALSE],
  surface = r,
  goal = 3,
  method = "gdistance"
)
out1$dist_lcp
out1$path_lcp
out2$dist_lcp
out2$path_lcp

#### Example (13): Shortest distances/paths computed in parallel:
# cppRouting method for goal 3
out1 <- lcp_over_surface(
  origin = rxy[1:2, , drop = FALSE],
  destination = rxy[5:6, , drop = FALSE],
  surface = r,
  goal = 3,
  use_all_cores = TRUE
)
# gdistance method for goal 2
out1 <- lcp_over_surface(
  origin = rxy[1:2, , drop = FALSE],
  destination = rxy[5:6, , drop = FALSE],
  surface = r,
  goal = 2,
  method = "gdistance",
  cl = parallel::makeCluster(2L)
)
# gdistance method for goal 3
out1 <- lcp_over_surface(
  origin = rxy[1:2, , drop = FALSE],
  destination = rxy[5:6, , drop = FALSE],
  surface = r,
  goal = 3,
  method = "gdistance",
  cl = parallel::makeCluster(2L)
)


############################################################################
#### Shortest distances/paths between all origin/destination combinations

#### Example (14) Compute all combinations via combination = "matrix"
# cppRouting goal 3
out1 <- lcp_over_surface(
  origin = rxy[1:2, , drop = FALSE],
  destination = rxy[5:6, , drop = FALSE],
  surface = r,
  goal = 3,
  cppRouting_algorithm = "phast",
  combination = "matrix"
)
out1$dist_euclid
out1$dist_lcp
out1$path_lcp
# cppRouting goal 3 parallelised
out1 <- lcp_over_surface(
  origin = rxy[1:2, , drop = FALSE],
  destination = rxy[5:6, , drop = FALSE],
  surface = r,
  goal = 3,
  combination = "matrix",
  cppRouting_algorithm = "phast",
  use_all_cores = TRUE
)
# gdistance goal 3
out1 <- lcp_over_surface(
  origin = rxy[1:2, , drop = FALSE],
  destination = rxy[5:6, , drop = FALSE],
  surface = r,
  goal = 3,
  method = "gdistance",
  combination = "matrix",
  cl = parallel::makeCluster(2L)
)
out1$dist_euclid
out1$dist_lcp
out1$path_lcp

#### Example (15): Real world example with multiple origins/destinations
## Zoom in on an area of interest
ext <- raster::extent(715000, 720000, 6250000, 6260000)
dat_gebco_utm_planar_zoom <- raster::crop(dat_gebco_utm_planar, ext)
## Define example origins/destinations
# Define available coordinates
dat_gebco_utm_planar_zoom_xy <- raster::coordinates(dat_gebco_utm_planar_zoom)
# Sample random origins
set.seed(2019)
index <- sample(1:nrow(dat_gebco_utm_planar_zoom_xy), 2)
origin <- dat_gebco_utm_planar_zoom_xy[index, ]
# Sample random destinations
set.seed(2020)
index <- sample(1:nrow(dat_gebco_utm_planar_zoom_xy), 3)
destination <- dat_gebco_utm_planar_zoom_xy[index, ]
# Implement algorithm
out_gebco1 <- lcp_over_surface(
  origin = origin,
  destination = destination,
  surface = dat_gebco_utm_planar,
  method = "gdistance",
  goal = 3,
  combination = "matrix"
)
out_gebco2 <- lcp_over_surface(
  origin = origin,
  destination = destination,
  surface = dat_gebco_utm_planar,
  method = "cppRouting",
  goal = 3,
  combination = "matrix",
  cppRouting_algorithm = "phast"
)
# The function returns distances between all combinations
out_gebco1$dist_lcp
out_gebco2$dist_lcp
# The two outputs are the same (though the ordering of factor levels differs)
str(out_gebco1$path_lcp$cells)
str(out_gebco2$path_lcp$cells)

}
\author{
Edward Lavender
}