R/hf_shortPath.R

In kbvernon/rHike: Functions for Computing Hiking Efficiency and Velocity

#' Least Cost Path Between Points
#'
#' A generalized method for calculating the "shortest" or
#' least-cost paths between multiple start and end points.
#'
#' @usage hf_shortPath(x, origin, goal, fun = function(x) 1/mean(x), add_cost = FALSE)
#'
#' @param x A \code{RasterLayer} with velocity or efficiency values.
#' @param origin A \code{POINT} object from package \code{sf} specifying start locations
#' @param goal A \code{POINT} object from package \code{sf} specifying end locations
#' @param fun A function to calculate transition weights from the raster.
#' @param add_cost Logical. Whether to return as a column the cost of travel between points.
#'
#' @return An \code{sf} object with least-cost line features. If \code{add_cost = TRUE},
#' the sf object also includes a cost column.
#' @export
#'
#' @examples
#' \dontrun{
#'
#' library(raster)
#' library(sf)
#'
#' dem <- raster(system.file("extdata/slc.tif", package = "rHike"))
#'
#' slope <- hf_slope(dem)
#'
#' velocity <- hf_velocity(slope, hf = "campbell", decile = 30)
#'
#' start_points <- st_sf(id = 1:3,
#'                       geometry = st_sfc(st_point(c(424350, 4514200)),
#'                                         st_point(c(426000, 4515000)),
#'                                         st_point(c(429000, 4516500)),
#'                                         crs = 26912))
#'
#' end_points <- st_sf(id = 1:2,
#'                     geometry = st_sfc(st_point(c(426200, 4516400)),
#'                                       st_point(c(427600, 4515600)),
#'                                       crs = 26912))
#'
#' short_paths <- hf_shortPath(velocity, start_points, end_points)
#'
#' plot(dem)
#' plot(st_geometry(start_points), color = "red", add = TRUE)
#' plot(st_geometry(end_points), color = "blue", add = TRUE)
#' plot(st_geometry(short_paths), lty = 2, add = TRUE)
#'
#' }
#'
hf_shortPath <- function(x, origin, goal, fun = function(x) 1/mean(x), add_cost = FALSE){

  origin <- sf::st_coordinates(origin)
  origin <- as.list(data.frame(t(origin))) # https://stackoverflow.com/questions/2471188/converting-a-matrix-to-a-list

  goal <- sf::st_coordinates(goal)

  # generate transition object from friction surface
  transition <- gdistance::transition(x, fun, 8)

  # account for diagonal travel in Moore's neighborhood
  transition <- gdistance::geoCorrection(transition)

  short_paths <- lapply(origin, function(z) gdistance::shortestPath(transition, z, goal, output = "SpatialLines"))

  # convert to sf object
  short_paths <- lapply(short_paths, sf::st_as_sf)

  # some 'paths' are single points, so convert to line with same start and end point
  # unfortunately, it's hard to iterate over sf geometries with apply functions
  short_paths <- lapply(short_paths, function(x){

    vertices <- lapply(sf::st_geometry(x), function(x) unlist(sf::st_geometry(x), use.names = FALSE))

    npts <- sapply(vertices, function(x) length(x) / 2)

    if(any(npts == 1)){

      for(i in which(npts == 1)){

        new_sfg <- matrix(rep(vertices[[i]], 2),
                          ncol = 2,
                          byrow = TRUE)

        new_sfg <- sf::st_linestring(new_sfg)

        x[i,] <- sf::st_set_geometry(x[i,], value = sf::st_sfc(new_sfg))

      }

    }

    return(x)

  })

  # combine list into one sf object (https://github.com/r-spatial/sf/issues/798)
  short_paths <- sf::st_as_sf(data.table::rbindlist(short_paths))

  # add indices
  short_paths <- cbind(short_paths,
                       origin = rep(seq_along(origin), each = nrow(goal)),
                       goal = rep(1:nrow(goal), times = length(origin)))

  if(add_cost){

    cost <- lapply(origin, function(z) gdistance::costDistance(transition, z, goal))

    short_paths <- cbind(short_paths,
                         cost = unlist(cost, use.names = FALSE))

  }

  return(short_paths)

}