R/movecomp.r

In movecost: Calculation of Slope-Dependant Accumulated Cost Surface, Least-Cost Paths, Least-Cost Corridors, Least-Cost Networks Related to Human Movement Across the Landscape

#' R function for comparing least-cost paths generated using different cost functions
#'
#' The function provides the facility to calculate LCPs using different cost functions and to plot them in the same visual output to allow
#' comparability. See, for instance, fig. 14.2 in Parcero-Oubina C. et al, Footprints and Cartwheels on a Pixel Road: On the Applicability of GIS for
#' the Modelling of Ancient (Roman) Routes (2019). In Verhagen P., Joyce J., Groenhuijzen M.R. (eds), Finding the
#' Limits of the Limes. Modelling Demography, Economy and Transport on the Edge of the Roman Empire, Springer,  291-311.\cr
#' Visit this \href{https://drive.google.com/file/d/1gLDrkZFh1b_glzCEqKdkPrer72JJ9Ffa/view?usp=sharing}{LINK} to access the package's vignette.\cr
#'
#' Like \code{movecost()}, the function just requires an input DTM ('RasterLayer' class), and an origin and destination dataset  ('SpatialPointsDataFrame' class).
#' The cost functions to be used have to be entered into 'movecomp()' via a character vector fed via the \code{choice} parameter (see the examples below).
#' If a DTM is not provided, \code{movecomp()} downloads elevation data from online sources for the area enclosed by the polygon fed via
#' the \code{studyplot} parameter (see \code{\link{movecost}} for more details). Under the hood, \code{movecomp()} relies on \code{movecost()} and implements the same cost functions:
#' see the help documentation of \code{movecost()} for further information.\cr
#'
#' \code{movecomp()} produces a plot representing the input DTM overlaid by a slopeshade raster, whose transparency can be adjusted using
#' the 'transp' parameter. On the rendered plot, the LPCs ('SpatialLinesDataFrame' class) generated by the different input cost functions are given a different line type; a legend
#' indicates which line type corresponds to which cost function. LCPs back to the origin can be calculated (and plotted) setting the
#' parameter \code{return.base} to \code{TRUE}.\cr
#'
#' The function returns the LCPs and (if requested by the user) the LCPs back to the origin. If the DTM has been acquired online, it will be returned as well.
#' The LCPs (and the LCPs back to the origin) will store three variables: the length of each path, the cost of each path, and an abbreviation corresponding to the cost function used to
#' generate the LCPs. The mentioned data can be exported by setting the \code{export} parameter to \code{TRUE}.\cr
#'
#' If the users want to compare the distribution of the length of the LCPs generate by different cost functions,
#' it suffices to set the \code{add.chart} parameter to \code{TRUE}. Two charts featuring boxplots will be rendered:
#' one plotting the distribution of the LCPs length by cost function; one portaying the distribution of the cost by cost function.\cr
#'
#' The following example uses in-built datasets to compare the LCPs generated using two cost functions: the Tobler hiking function , the
#' wheeled vehicle cost function, and the Pantolf. et al's cost function with correction factor.
#' LCPs back to the origin location will be calculated as well. The origin and destination locations are close to Mt Etna (Sicily, Italy).
#' Note that elevation data are acquired online for the area enclosed by the polygon fed via the \code{studyplot} parameter:\cr
#'
#' result <- movecomp(origin=Etna_start_location, destin=Etna_end_location, choice=c("t", "wcs", "pcf"), studyplot = Etna_boundary, return.base=TRUE)\cr
#'
#'
#' @param dtm Digital Terrain Model (RasterLayer class); if not provided, elevation data will be acquired online for the area enclosed by the 'studyplot' parameter (see \code{\link{movecost}}).
#' @param origin location from which least-cost path(s) is calculated (SpatialPointsDataFrame class).
#' @param destin location(s) to which least-cost path(s) is calculated (SpatialPointsDataFrame class).
#' @param studyplot polygon (SpatialPolygonDataFrame class) representing the study area for which online elevation data are acquired (see \code{\link{movecost}}); NULL is default.
#' @param barrier area where the movement is inhibited (SpatialLineDataFrame or SpatialPolygonDataFrame class) (see \code{\link{movecost}}).
#' @param plot.barrier TRUE or FALSE (default) if the user wants or does not want the barrier to be plotted (see \code{\link{movecost}}).
#' @param irregular.dtm TRUE or FALSE (default) if the input DTM features irregular margins (see \code{\link{movecost}}).
#' @param choice character vector indicating the cost functions to be compared (for details on each of the following, see \code{\link{movecost}}):\cr
#'
#' \strong{-functions expressing cost as walking time-}\cr
#' \strong{t} (default) uses the on-path Tobler's hiking function;\cr
#' \strong{tofp} uses the off-path Tobler's hiking function;\cr
#' \strong{mp} uses the Marquez-Perez et al.'s modified Tobler's function;\cr
#' \strong{icmonp} uses the Irmischer-Clarke's hiking function (male, on-path);\cr
#' \strong{icmoffp} uses the Irmischer-Clarke's hiking function (male, off-path);\cr
#' \strong{icfonp} uses the Irmischer-Clarke's hiking function (female, on-path);\cr
#' \strong{icfoffp} uses the Irmischer-Clarke's hiking function (female, off-path);\cr
#' \strong{ug} uses the Uriarte Gonzalez's walking-time cost function;\cr
#' \strong{ma} uses the Marin Arroyo's walking-time cost function;\cr
#' \strong{alb} uses the Alberti's Tobler hiking function modified for pastoral foraging excursions;\cr
#' \strong{gkrs} uses the Garmy, Kaddouri, Rozenblat, and Schneider's hiking function;\cr
#' \strong{r} uses the Rees' hiking function;\cr
#' \strong{ks} uses the Kondo-Seino's hiking function;\cr
#' \strong{trp} uses the Tripcevich's hiking function;\cr
#'
#' \strong{-functions for wheeled-vehicles-}\cr
#' \strong{wcs} uses the wheeled-vehicle critical slope cost function;\cr
#'
#' \strong{-functions expressing abstract cost-}\cr
#' \strong{ree} uses the relative energetic expenditure cost function;\cr
#' \strong{b} uses the Bellavia's cost function;\cr
#' \strong{e} uses the Eastman's cost function;\cr
#'
#' \strong{-functions expressing cost as metabolic energy expenditure-}\cr
#' \strong{p} uses the Pandolf et al.'s metabolic energy expenditure cost function;\cr
#' \strong{pcf} uses the Pandolf et al.'s cost function with correction factor for downhill movements;\cr
#' \strong{m} uses the Minetti et al.'s metabolic energy expenditure cost function;\cr
#' \strong{hrz} uses the Herzog's metabolic energy expenditure cost function;\cr
#' \strong{vl} uses the Van Leusen's metabolic energy expenditure cost function;\cr
#' \strong{ls} uses the Llobera-Sluckin's metabolic energy expenditure cost function;\cr
#' \strong{a} uses the Ardigo et al.'s metabolic energy expenditure cost function;\cr
#' \strong{h} uses the Hare's metabolic energy expenditure cost function (for all the mentioned cost functions, see \code{\link{movecost}}).\cr
#' @param time time-unit to be used if Tobler's and other time-related cost functions are used; h' for hour, 'm' for minutes;
#' @param move number of directions in which cells are connected: 4 (rook's case), 8 (queen's case), 16 (knight and one-cell queen moves; default).
#' @param field value assigned to the cells coinciding with the barrier (0 by default) (see \code{\link{movecost}}.
#' @param cogn.slp  TRUE or FALSE (default) if the user wants or does not want the 'cognitive slope' to be used in place of the real slope (see \code{\link{movecost}}).
#' @param sl.crit critical slope (in percent), typically in the range 8-16 (10 by default) (used by the wheeled-vehicle cost function; see \code{\link{movecost}}).
#' @param W walker's body weight (in Kg; 70 by default; used by the Pandolf's and Van Leusen's cost function; see \code{\link{movecost}}).
#' @param L carried load weight (in Kg; 0 by default; used by the Pandolf's and Van Leusen's cost function; see \code{\link{movecost}}).
#' @param N coefficient representing ease of movement (1 by default) (see \code{\link{movecost}}).
#' @param V speed in m/s (1.2 by default) (used by the Pandolf et al.'s, Pandolf et al.s with correction factor, Van Leusen's, and Ardigo et al.'s cost function; if set to 0, it is internally worked out on the basis of Tobler on-path hiking function (see \code{\link{movecost}}).
#' @param z zoom level for the elevation data downloaded from online sources (from 0 to 15; 9 by default) (see \code{\link{movecost}} and \code{\link[elevatr]{get_elev_raster}}).
#' @param return.base TRUE or FALSE (default) if the user wants or does not want the least-cost paths back to the origin to be calculated and plotted (as dashed lines).
#' @param leg.pos set the position of the legend in the plotted cost allocation raster; 'topright' by default (other options: "bottomright", "bottom", "bottomleft", "left", "topleft", "top", "topright", "right", "center").
#' @param leg.cex set the size of the labels used in the legend displayed in the rendered plot (0.75 by default).
#' @param transp set the transparency of the slopeshade raster that is plotted over the DTM (0.5 by default).
#' @param add.chart TRUE or FALSE (default) is the user wants or does not want boxplots visualising LCPs length/cost vs cost function to be rendered.
#' @param oneplot TRUE (default) or FALSE if the user wants or does not want the plots displayed in a single window.
#' @param export TRUE or FALSE (default) if the user wants or does not want the LCPs to be exported as a shapefile; the DTM is exported only if it was not provided by the user
#' and downloaded by the function from online sources.
#'
#' @return The function returns a list storing the following components \itemize{
##'  \item{dtm: }{Digital Terrain Model ('RasterLayer' class); returned only if acquired online}
##'  \item{LCPs: }{estimated least-cost paths ('SpatialLinesDataFrame' class); three variables are stored: 'length', 'cost', and
##'  'funct'.}
##'  \item{LCPs.back: }{estimated least-cost paths back to the origin ('SpatialLinesDataFrame' class); three variables are stored; see above.}
##' }
##'
#' @keywords movecomp
#' @export
#' @importFrom raster ncell raster terrain
#' @importFrom terra writeVector vect
#' @importFrom elevatr get_elev_raster
#' @importFrom grDevices terrain.colors topo.colors grey
#' @importFrom graphics boxplot
#'
#' @examples
#' # load a sample Digital Terrain Model
#' data(volc)
#'
#'
#' # load the sample destination locations on the above DTM
#' data(destin.loc)
#'
#'
#' # compare the LCPs generated using different walking-time cost functions (time in minutes)
#'
#' result <- movecomp(volc, volc.loc, destin.loc, choice=c("t", "ug", "gkrs"), time="m", move=8)
#'
#' # the distribution of the length and cost of the LCPs by cost function can be easily compared
#' # using the 'add.chart' parameter:
#' #result <- movecomp(volc, volc.loc, destin.loc, choice=c("t", "ug", "gkrs"), time="m",
#' #move=8, add.chart=T)
#'
#'
#' @seealso \code{\link{movecost}}
#'
#'
movecomp <- function (dtm=NULL, origin, destin, studyplot=NULL, barrier=NULL, plot.barrier=FALSE, irregular.dtm=FALSE, choice, time="h", move=16, field=0, cogn.slp=FALSE, sl.crit=10, W=70, L=0, N=1, V=1.2, z=9, return.base=FALSE, leg.pos="topright", leg.cex=0.75, transp=0.5, add.chart=FALSE, oneplot=TRUE, export=FALSE){

  #if no dtm is provided
  if (is.null(dtm)==TRUE) {
    #get the elvation data using the elevatr's get_elev_raster() function, using the studyplot dataset (SpatialPolygonDataFrame)
    #to select the area whose elevation data are to be downloaded;
    #z sets the resolution of the elevation datataset
    elev.data <- elevatr::get_elev_raster(studyplot, z = z, verbose=FALSE, override_size_check = TRUE)
    #crop the elevation dataset to the exact boundary of the studyplot dataset
    dtm <- raster::crop(elev.data, studyplot)
  }

  #get the number of elements in the vector of cost-functions entered by the user
  n <- length(choice)

  #create an empty list where the results will be stored
  lcp.paths <- list()

  #run the 'movecost' function one time for each number of cost-functions to be compared
  #and store the LPCs in the list previously created;
  #also, at each loop, store the corresponding function's cost and abbreviation in the newly created 'cost' and 'funct' field
  for (i in 1:n) {
    message(paste0("Wait...calculating the LCP(s) on the basis of '", choice[i],"' function"))
    output.data <- movecost::movecost(dtm=dtm, origin=origin, destin=destin, barrier=barrier, irregular.dtm=irregular.dtm, funct=choice[i], time=time, move=move, field=field, cogn.slp=cogn.slp, sl.crit=sl.crit, W=W, L=L, N=N, V=V, z=z, return.base=FALSE, graph.out=FALSE)
    lcp.paths[[i]] <- output.data$LCPs
    lcp.paths[[i]]$cost <- output.data$dest.loc.w.cost$cost
    lcp.paths[[i]]$funct <- choice[i]
  }

  #merge all the LCPs; each one will be identified by the function's abbreviation stored in the 'funct' field,
  #as per previous step
  merged.paths <- do.call(rbind, lcp.paths)

  #if the user wants the LCPs back to the origin...
  if(return.base==TRUE) {

    #create two empty listes to be used later on
    lcp.paths.back <- list()
    merged.path.back.temp <- list()

    #use two nested loops; the inner one loops through the destinations using them as origin
    #this proves admittedly slower than the previous loop
    for (i in 1:n) {
      message(paste0("Wait...calculating the LCP(s) back to the origin on the basis of '", choice[i],"' function"))
      for (j in 1:length(destin)) {
        output.data <- movecost::movecost(dtm=dtm, origin=destin[j,], destin=origin, barrier=barrier, irregular.dtm=irregular.dtm, funct=choice[i], time=time, move=move, field=field, cogn.slp=cogn.slp, sl.crit=sl.crit, W=W, L=L, N=N, V=V, z=z, return.base=TRUE, graph.out=FALSE)
        lcp.paths.back[[j]] <- output.data$LCPs
        lcp.paths.back[[j]]$cost <- output.data$dest.loc.w.cost$cost
        lcp.paths.back[[j]]$funct <- choice[i]
        merged.path.back.temp[[i]] <- do.call(rbind, lcp.paths.back)
      }
    }

    #merge all the LCPs back to the origin
    merged.paths.back <- do.call(rbind, merged.path.back.temp)
  }

  #produce the ingredient for the slopeshade raster
  #to be used in both the rendered plots
  slope <- raster::terrain(dtm, opt = "slope")

  #if the user wants the LCPs back to the origin AND 'oneplot' is TRUE,
  #set the layout in just one visualization
  if(return.base==TRUE & oneplot==TRUE){
    m <- rbind(c(1,2))
    layout(m)
  }

  #plot the DTM raster
  raster::plot(dtm,
               main="Comparison of Least-Cost Paths",
               sub="black dot=start location\n red dot(s)=destination location(s)",
               cex.main=0.95,
               cex.sub=0.75,
               legend.lab="Elevation (masl)")

  #add the slopeshade
  raster::plot(slope,
               col = rev(grey(0:100/100)),
               legend = FALSE,
               alpha=transp,
               add=TRUE)

  #add the (merged) LCPs
  #create a loop to plot the LCPs by function (i.e., iteratively subsetting the 'merged.paths' dataset
  #and plotting the subsets additively)
  for (i in 1:n) {
    raster::plot(merged.paths[merged.paths$funct==choice[i],], add=T, lty=i)
  }

  #add the legend
  legend(leg.pos,
         legend=choice,
         lty=1:n,
         cex=leg.cex)

  #add the origin
  raster::plot(origin,
               pch=20,
               add=TRUE)

  #add the destinations
  raster::plot(destin,
               pch=20,
               col="red",
               add=TRUE)

  #if the barrier is provided AND if plot.barrier is TRUE, add the barrier
  if(is.null(barrier)==FALSE & plot.barrier==TRUE) {
    raster::plot(barrier, col="blue", add=TRUE)
  }

  #if the user wants the LCPs back to the origin
  if(return.base==TRUE) {
    #plot the DTM raster
    raster::plot(dtm,
                 main="Comparison of Least-Cost Paths back to the origin",
                 sub="black dot=start location\n red dot(s)=destination location(s)",
                 cex.main=0.95,
                 cex.sub=0.75,
                 legend.lab="Elevation (masl)")

    #add the slopeshade
    raster::plot(slope,
                 col = rev(grey(0:100/100)),
                 legend = FALSE,
                 alpha=transp,
                 add=TRUE)

    #add the (merged) LCPs abck to the origin
    #create a loop to plot the LCPs by function (i.e., iteratively subsetting the 'merged.paths.back' dataset
    #and plotting the subsets additively)
    for (i in 1:n) {
      raster::plot(merged.paths.back[merged.paths.back$funct==choice[i],], add=T, lty=i)
    }

    #add the legend
    legend(leg.pos,
           legend=choice,
           lty=1:n,
           cex=leg.cex)

    #add the origin
    raster::plot(origin,
                 pch=20,
                 add=TRUE)

    #add the destinations
    raster::plot(destin,
                 pch=20,
                 col="red",
                 add=TRUE)

    #if the barrier is provided AND if plot.barrier is TRUE, add the barrier
    if(is.null(barrier)==FALSE & plot.barrier==TRUE) {
      raster::plot(barrier, col="blue", add=TRUE)
    }
  }

  #if 'add.chart' is TRUE, render the boxplots
  if(add.chart==TRUE) {
    graphics::boxplot(merged.paths$length ~ merged.paths$funct,
                      xlab="cost function",
                      ylab="length",
                      main="Boxplot of LCPs length by cost function",
                      cex.main=0.85,
                      cex.lab=0.95)

    graphics::boxplot(merged.paths$cost ~ merged.paths$funct,
                      xlab="cost function",
                      ylab="cost",
                      main="Boxplot of LCPs cost by cost function",
                      cex.main=0.85,
                      cex.lab=0.95)

    #if 'return.base' is TRUE, render the boxplots for the LCPs back to the origin
    if(return.base==TRUE) {
      graphics::boxplot(merged.paths.back$length ~ merged.paths.back$funct,
                        xlab="cost function",
                        ylab="length",
                        main="Boxplot of length by cost function (LCPs back to the origin)",
                        cex.main=0.85,
                        cex.lab=0.95)

      graphics::boxplot(merged.paths.back$cost ~ merged.paths.back$funct,
                        xlab="cost function",
                        ylab="cost",
                        main="Boxplot of cost by cost function (LCPs back to the origin)",
                        cex.main=0.85,
                        cex.lab=0.95)
    }
  }

  #if export is TRUE, export the LCPs as a shapefile
  if(export==TRUE){
    terra::writeVector(vect(merged.paths), filename="LCPs", filetype="ESRI Shapefile")
  }

  #if export is TRUE and return.base is also TRUE, export the LCPs back to the origin as a shapefile
  if(export==TRUE & return.base==TRUE){
    terra::writeVector(vect(merged.paths.back), filename="LCPs_back", filetype="ESRI Shapefile")
  }

  #if no DTM was provided (i.e., if 'studyplot' is not NULL), export the downloaded DTM as a raster file
  if(export==TRUE & is.null(studyplot)==FALSE){
    raster::writeRaster(dtm, "dtm", format="GTiff")
  }

  #if studyplot is NULL (i.e., if the DTM has been provided)..
  if(is.null(studyplot)==TRUE){
    #set dtm to NULL, i.e. there is no DTM to return
    dtm <- NULL
  } else{
    #otherwise (i.e., if studyplot was not NULL), set the dtm to the downloaded dtm (i.e., there is
    #something to return)
    dtm <- dtm
  }

  #same as the step above with regards to the LCPs back
  if(return.base==TRUE){
    merged.paths.back <- merged.paths.back
  } else{
    merged.paths.back <- NULL
  }

  #restore the original graphical device's settings if previously modified
  if(return.base==TRUE & oneplot==TRUE){
    par(mfrow = c(1,1))
  }

  results <- list("dtm"=dtm,
                  "LCPs"=merged.paths,
                  "LPCs.back"=merged.paths.back)

  return(results)
}