Nothing
#' 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)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.