R/movecomp.R

#' 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
#'
#' Like '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 'choice' parameter (see the examples below).
#' If a DTM is not provided, 'movecomp()' downloads elevation data from online sources for the area enclosed by the polygon fed via
#' the 'studyplot' parameter (see \code{\link{movecost}} for more details). Under the hood, movecomp()' relies on 'movecost()' and implements the same cost functions:
#' see the help documentation of 'movecost()' for further information.\cr
#'
#' 'movecomp()' produces a plot representing the input DTM overlaid by an hillshade 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 'return.base' to TRUE.\cr
#'
#' The function returns the LCPs and (if requested by the user) the LCPs back to the origin. If the DTM has been aquired online, it will be returned as well.
#' The LCPs (and the LCPs back to the origin) will store two attributes: the length of each path and an abbreviation corresponding to the cost function used to
#' generate each LCPs. This might prove handy if the user wants to compare the distribution of the length of groups of LCPs generate by different cost functions
#' (see the example in the 'Examples' section further down).
#' The mentioned data can be exported by setting the 'export' parameter to TRUE.\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 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(s) around which the boundary(ies) 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 aquired (see \code{\link{movecost}}); NULL is default.
#' @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{-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{-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 (for all the mentioned cost functions, see Details);\cr
#' @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 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 cex.leg 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 hillshade raster that is plotted over the DTM (0.5 by default).
#' @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 aquired online}
##'  \item{LCPs: }{estimated least-cost paths ('SpatialLinesDataFrame' class)}
##'  \item{LCPs.back: }{estimated least-cost paths back to the origin ('SpatialLinesDataFrame' class)}
##' }
##'
#' @keywords movecomp
#' @export
#' @importFrom raster ncell raster hillShade terrain
#' @importFrom elevatr get_elev_raster
#' @importFrom grDevices terrain.colors topo.colors grey
#'
#' @examples
#' # load a sample Digital Terrain Model
#' volc <- raster::raster(system.file("external/maungawhau.grd", package="gdistance"))
#'
#'
#' # load the sample destination locations on the above DTM
#' data(destin.loc)
#'
#'
#' # compare the LCPs generated using the Tobler hiking function (for on-path walk)
#' # and the wheeled vehicle cost function
#'
#' result <- movecomp(dtm=volc, origin=volc.loc, destin=destin.loc, choice=c("t", "wcs"), move=16)
#'
#' #the distribution of the length of the LCPs can be easily compared using boxplots for instance:
#' boxplot(result$LCPs$length ~ result$LCPs$funct, xlab="cost function", ylab="length (meter)")
#'
#'
#' @seealso \code{\link{movecost}}
#'
#'
movecomp <- function (dtm=NULL, origin, destin, studyplot=NULL, choice, move=16, cogn.slp=FALSE, sl.crit=10, W=70, L=0, N=1, V=1.2, z=9, return.base=FALSE, leg.pos="topright", cex.leg=0.75, transp=0.5, 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()

  #set the time to "h"; this is not essential, but need to be set so that the 'movecost()' used
  #inside loops to come can work
  time <- "h"

  #run the 'movecost' function one time for each number of cost-function to be compared
  #and store the LPCs in the list previously created;
  #also, at each loop, store the corresponding function's abbreviation in the newly created 'funct' field

  message("Wait...calculating the LCPs...")

  for (i in 1:n) {
    lcp.paths[[i]] <- movecost::movecost(dtm=dtm, origin=origin, destin=destin, funct=choice[i], time=time, move=move, cogn.slp=cogn.slp, sl.crit=sl.crit, W=W, L=L, N=N, V=V, z=z, return.base=FALSE, graph.out=FALSE)$LCPs
    lcp.paths[[i]]$funct <- choice[i]
  }

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

  #if the user wants the LCPs back to the origin,
  #do the same as before, setting the movecost's parameter 'return.base' to TRUE
  if(return.base==TRUE) {
    message("Wait...calculating the LCPs back to the origin...")
    lcp.paths.back <- list()

    for (i in 1:n) {
      lcp.paths.back[[i]] <- movecost::movecost(dtm=dtm, origin=origin, destin=destin, funct=choice[i], time=time, move=move, cogn.slp=cogn.slp, sl.crit=sl.crit, W=W, L=L, N=N, V=V, z=z, return.base=TRUE, graph.out=FALSE)$LCPs.back
      lcp.paths.back[[i]]$funct <- choice[i]
    }

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

  #produce the ingredients for the hillshade raster
  #to be used in both the rendered plots
  slope <- raster::terrain(dtm, opt = "slope")
  aspect <- raster::terrain(dtm, opt = "aspect")
  hill <- raster::hillShade(slope, aspect, angle = 45, direction = 0)


  #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 hillshade
  raster::plot(hill,
               col = grey(0:100/100),
               legend = FALSE,
               alpha=transp,
               add=TRUE)

  #add the (merged) LCPs
  raster::plot(merged.paths,
               add=T,
               lty=as.numeric(as.factor(merged.paths$funct)))

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

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

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

  #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 hillshade
    raster::plot(hill,
                 col = grey(0:100/100),
                 legend = FALSE,
                 alpha=transp,
                 add=TRUE)

    #add the (merged) LCPs back to the origin
    raster::plot(merged.paths.back,
                 add=T,
                 lty=as.numeric(as.factor(merged.paths.back$funct)))

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

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

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

  #if export is TRUE, export the LCPs as a shapefile
  if(export==TRUE){
    rgdal::writeOGR(merged.paths, ".", "LCPs", driver="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){
    rgdal::writeOGR(merged.paths.back, ".", "LCPs_back", driver="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)
}

Try the movecost package in your browser

Any scripts or data that you put into this service are public.

movecost documentation built on Sept. 13, 2021, 1:06 a.m.