The function provides the facility to calculate the anisotropic accumulated cost of movement around a starting location and to optionally calculate least-cost path(s) toward
one or multiple destinations. It implements different cost estimations related to human movement across the landscape.
The function takes as input a Digital Terrain Model ('RasterLayer' class) and a point feature ('SpatialPointsDataFrame' class), the latter representing
the starting location, i.e. the location from which the accumulated cost is calculated.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28
movecost( dtm = NULL, origin, destin = NULL, studyplot = NULL, funct = "t", time = "h", outp = "r", move = 16, cogn.slp = FALSE, sl.crit = 10, W = 70, L = 0, N = 1, V = 1.2, z = 9, return.base = FALSE, rb.lty = 2, breaks = NULL, cont.lab = TRUE, destin.lab = TRUE, cex.breaks = 0.6, cex.lcp.lab = 0.6, graph.out = TRUE, transp = 0.5, oneplot = TRUE, export = FALSE )
Digital Terrain Model (RasterLayer class); if not provided, elevation data will be acquired online for the area enclosed by the 'studyplot' parameter (see Details).
location from which the cost surface is calculated (SpatialPointsDataFrame class).
location(s) to which least-cost path(s) is calculated (SpatialPointsDataFrame class).
polygon (SpatialPolygonDataFrame class) representing the study area for which online elevation data are aquired (see Details); NULL is default.
cost function to be used:
time-unit expressed by the accumulated raster and by the isolines if Tobler's and other time-related cost functions are used; 'h' for hour, 'm' for minutes.
type of output: 'raster' or 'contours' (see Details).
number of directions in which cells are connected: 4 (rook's case), 8 (queen's case), 16 (knight and one-cell queen moves; default).
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 Details).
critical slope (in percent), typically in the range 8-16 (10 by default) (used by the wheeled-vehicle cost function; see Details).
walker's body weight (in Kg; 70 by default; used by the Pandolf's and Van Leusen's cost function; see Details).
carried load weight (in Kg; 0 by default; used by the Pandolf's and Van Leusen's cost function; see Details).
coefficient representing ease of movement (1 by default) (see Details).
speed in m/s (1.2 by default) (used by the Pandolf's and Van Leusen's cost function; see Details).
zoom level for the elevation data downloaded from online sources (0 to 15; 9 by default) (see Details and
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).
Line type used to represent the least-cost paths back to the origin in the returned plot (2 by default; dashed line; see 'lty' parameter in
isolines interval; if no value is supplied, the interval is set by default to 1/10 of the range of values of the accumulated cost surface.
if set to TRUE (default) display the labels of the contours over the accumulated cost surface.
if set to TRUE (default) display the label(s) indicating the cost at the destination location(s).
set the size of the time labels used in the isochrones plot (0.6 by default).
set the size of the labels used in least-cost path(s) plot (0.6 by default).
TRUE (default) or FALSE if the user wants or does not want a graphical output to be generated.
set the transparency of the hillshade raster that is plotted over the rendered plots (0.5 by default).
TRUE (default) or FALSE if the user wants or does not want the plots displayed in a single window.
TRUE or FALSE (default) if the user wants or does not want the outputs to be exported; if TRUE, the DTM, the cost-surface, and the accumulated cost surface are exported as a GeoTiff file, while the isolines, the least-cost path(s), and a copy of the input destination locations (storing the cost measured at each location) are exported as shapefile; all the exported files (excluding the DTM) will bear a suffix corresponding to the cost function selected by the user. Note that the DTM is exported only if it was not provided by the user and downloaded by the function from online sources.
If the parameter 'destin' is fed with a dataset representing destination location(s) ('SpatialPointsDataFrame' class), the function also calculates
least-cost path(s) plotted on the input DTM; the length of each path is saved under the variable 'length' stored in the 'LCPs' dataset ('SpatialLines' class) returned by the function.
In the produced plot, the red dot(s) representing the destination location(s) are labelled with numeric values representing
the cost value at the location(s).
The cost value is also appended to the updated destination locations dataset returned by the function, which
stores a new variable labelled 'cost'. If the cost is expressed in terms of walking time, the labels accompaining each destinaton location will
express time in sexagesimal numbers (hours, minutes, seconds). In this case, the variable 'cost' appended to the returned destination location datset
will store the time figures in decimal numbers, while another variable named 'cost_hms' will store the corresponding value in sexagesimal numbers.
When interpreting the time values stored in the 'cost' variable, the user may want to bear in mind the selected time unit (see right below).
When using cost functions expressing cost in terms of time, the time unit can be selected by the user setting the 'time' parameter to 'h' (hours) or to 'm' (minutes).
In general, the user can also select which type of visualization the function has to produce; this is achieved setting the 'outp' parameter to either 'r' (=raster)
or to 'c' (=contours). The former will produce a raster with a colour scale and contour lines representing the accumulated cost surface; the latter parameter will only
produce contour lines.
The contour lines' interval is set using the parameter 'breaks'; if no value is passed to the parameter, the interval will be set by default to
1/10 of the range of values of the accumulated cost surface.
It is worth reminding the user(s) that all the input layers (i.e., DTM, start location, and destination locations) must use the same projected coordinate system.
Cost surface calculation:
for the cost-surface and LCPs calculation, 'movecost' builds on functions from Jacob van Etten's gdistance package. Under the hood, 'movecost' calculates the slope as rise over run, following the procedure described by van Etten, "R Package gdistance: Distances and Routes on Geographical Grids" in Journal of Statistical Software 76(13), 2017, pp. 14-15. The number of directions in which cells are connected in the cost calculation can be set to 4 (rook's case), 8 (queen's case), or 16 (knight and one-cell queen moves) using the 'move' parameter (see 'Arguments').
Acquiring online elevation data:
if a DTM is not provided, 'movecost()' will download elevation data from online sources. Elevation data will be aquired for the area enclosed by the polygon supplied by the 'studyplot' parameter (SpatialPolygonDataFrame class). To tap online elevation data, 'movecost' internally builds on the
get_elev_raster function from the elevatr package.
The zoom level of the downloaded DTM (i.e., its resolution) is controlled by the parameter 'z', which is
set to 9 by default (a trade off between resolution and download time).
To test this facility, the user may want to try the following code, that will generate a least-cost surface and least-cost paths
in an area close the Mount Etna (Sicily, Italy), whose elevation data are acquired online; the start and end locations, and the
polygon defining the study area, are provided in this same package:
result <- movecost(origin=Etna_start_location, destin=Etna_end_location, studyplot=Etna_boundary)
The LCPs back to the origin can be calculated and plotted setting the parameter 'return.base' to TRUE:
result <- movecost(origin=Etna_start_location, destin=Etna_end_location, studyplot=Etna_boundary, return.base=TRUE)
To know more about what elevation data are tapped from online
sources, visit: https://cran.r-project.org/web/packages/elevatr/vignettes/introduction_to_elevatr.html.
For more information about the elevation data resolution per zoom level, visit
To know what is sourced at what zoom level, visit
Terrain slope and cognitive slope:
when it comes to the terrain slope, the function provides the facility to use the so-called 'cognitive slope', following Pingel TJ (2013), Modeling Slope as a Contributor to Route Selection in Mountainous Areas, in Cartography and Geographic Information Science, 37(2), 137-148. According to Pingel, "Humans tend to overestimate geographic slopes by a surprisingly high margin...This analysis indicates downhill slopes are overestimated at approximately 2.3 times the vertical, while uphill slopes are overestimated at 2 times the vertical.". As a result, if the parameter 'cogn.slp' is set to TRUE, positive slope values are preliminarily multiplied by 1.99, while negative slope values are multiplied by 2.31.
Terrain factor (N):
virtually all the implemented cost functions (with few exceptions) can take into account a 'terrain factor' (N parameter; 1 by default), which represents the easiness/difficulty of moving on different terrain types. According to the type of terrain, the energy spent when walking increases. The same holds true for time, which increases because on a loose terrain (for instance) the walking speed decreases. While a terrain factor is 'natively' part of the Van Leusen's, Pandolf et al.'s, and Bellavia's cost function, it has been integrated into the other cost functions as well (when/if relevant).
Note that the terrain factor has NOT been implemented in the Alberti's, Tobler's off-path, and Irmischer-Clarke's off-path cost function.
As for the latter two, they already natively feature a terrain factor. Therefore, it has been implemented only in their on-path version.
Needless to say, if we use a terrain factor of 1.67 with the Tobler's (on-path) hiking function, the results
will be equal to those obtained using the Tobler's off-path function (the reciprocal of 1.67, i.e. 0.60, is in fact
natively used by the Tobler's function for off-path hiking). In fact, compare the results of the following two runs
of 'movecost()' (first some data are loaded):
volc <- raster::raster(system.file("external/maungawhau.grd", package="gdistance"))
result1 <- movecost(dtm=volc, origin=volc.loc, destin=destin.loc, breaks=0.05, funct="t", N=1.67)
result2 <- movecost(dtm=volc, origin=volc.loc, destin=destin.loc, breaks=0.05, funct="tofp")
The user may want to refer to the following list of terrain factors, which is based on the data collected in Herzog, I. (2020).
Spatial Analysis Based on Cost Functions. In Gillings M, Haciguzeller P, Lock G (eds), "Archaeological Spatial Analysis. A Methodological Guide.",
Routledge: New York, 340 (with previous references). The list is divided into two sections (a and b), the first reporting the terrain
factors to be used for cost functions measuring time, the second for functions measuring cost other than time:
Blacktop roads, improved dirt paths, cement = 1.00
Lawn grass = 1.03
Loose beach sand = 1.19
Disturbed ground (former stone quarry) = 1.24
Horse riding path, flat trails and meadows = 1.25
Tall grassland (with thistle and nettles) = 1.35
Open space above the treeline (i.e., 2000 m asl) = 1.50
Bad trails, stony outcrops and river beds = 1.67
Off-paths = 1.67
Bog = 1.79
Off-path areas below the treeline (pastures, forests, heathland) = 2.00
Rock = 2.50
Swamp, water course = 5.00
Asphalt/blacktop = 1.00
Dirt road or grass = 1.10
Hard-surface road = 1.20
Light brush = 1.20
Ploughed field = 1.30 or 1.50
Heavy brush = 1.50
Hard-packed snow = 1.60
Swampy bog = 1.80
Sand dunes = 1.80
Loose sand = 2.10
Besides citing this package, you may want to refer to the following journal article:
Alberti (2019) <doi:10.1016/j.softx.2019.100331>.
Implemented cost functions:
note that in what follows x[adj] stands for slope as rise/run calculated for adjacent cells:
Tobler's hiking function (on-path) (speed in kmh):
6 * exp(-3.5 * abs(x[adj] + 0.05))
Tobler's hiking function (off-path) (speed in kmh):
(6 * exp(-3.5 * abs(x[adj] + 0.05))) * 0.6
as per Tobler's indication, the off-path walking speed is reduced by 0.6.
Marquez-Perez et al.'s modified Tobler hiking function (speed in kmh):
4.8 * exp(-5.3 * abs((x[adj] * 0.7) + 0.03))
modified version of the Tobler's hiking function as proposed by Joaquin Marquez-Parez, Ismael Vallejo-Villalta & Jose I. Alvarez-Francoso (2017), "Estimated travel time for walking trails in natural areas",
Geografisk Tidsskrift-Danish Journal of Geography, 117:1, 53-62, DOI: 10.1080/00167223.2017.1316212.
Irmischer-Clarke's modified Tobler hiking function (male, on-path; speed in kmh):
(0.11 + exp(-(abs(x[adj])*100 + 5)^2 / (2 * 30)^2)) * 3.6
modified version of the Tobler's function as proposed for (male) on-path hiking by Irmischer, I. J., & Clarke, K. C. (2018). Measuring and modeling the speed of human navigation.
Cartography and Geographic Information Science, 45(2), 177-186. https://doi.org/10.1080/15230406.2017.1292150. Note: all the all the Irmischer-Clarke's functions
originally express speed in m/s; they have been reshaped (multiplied by 3.6) to turn m/s into km/h for consistency with the other Tobler-related cost functions; slope is in percent.
Irmischer-Clarke's modified Tobler hiking function (male, off-path; speed in kmh):
(0.11 + 0.67 * exp(-(abs(x[adj])*100 + 2)^2 / (2 * 30)^2)) * 3.6
Irmischer-Clarke's modified Tobler hiking function (female, on-path; speed in kmh):
(0.95 * (0.11 + exp(-(abs(x[adj]) * 100 + 5)^2/(2 * 30^2)))) * 3.6
Irmischer-Clarke's modified Tobler hiking function (female, off-path; speed in kmh):
(0.95 * (0.11 + 0.67 * exp(-(abs(x[adj]) * 100 + 2)^2/(2 * 30^2)))) * 3.6
Uriarte Gonzalez's slope-dependant walking-time cost function:
1/ (0.0277 * (abs(x[adj])*100) + 0.6115)
proposed by Uriarte Gonzalez;
see: Chapa Brunet, T., Garcia, J., Mayoral Herrera, V., & Uriarte Gonzalez, A. (2008). GIS landscape models for the study of preindustrial settlement patterns in Mediterranean areas.
In Geoinformation Technologies for Geo-Cultural Landscapes (pp. 255-273). CRC Press. https://doi.org/10.1201/9780203881613.ch12.
The cost function is originally expressed in seconds; for the purpose of its implementation in this function, it is the reciprocal of time (1/T) that is used in order to eventually get T/1. Unlike the original cost function, here the pixel resolution is not taken into account since 'gdistance' takes care of the cells' dimension when calculating accumulated costs.
Alberti's Tobler hiking function modified for pastoral foraging excursions (speed in kmh):
(6 * exp(-3.5 * abs(x[adj] + 0.05))) * 0.25
proposed by Gianmarco Alberti; see: Locating potential pastoral foraging routes in Malta through the use of a Geographic Information System. The Tobler’s function has been rescaled to fit animal walking speed during foraging excursions. The distribution of the latter, as empirical data show, turns out to be right-skewed and to vary along a continuum. It ranges from very low speed values (corresponding to slow grazing activities grazing while walking) to comparatively higher values (up to about 4.0 km/h) corresponding to travels without grazing (directional travel toward feeding stations). The function consider 1.5 km/h as the average flock speed, which roughly corresponds to the average speed recorded in some studies, and can be considered the typical speed of flocks during excursions in which grazing takes place while walking (typical form of grazing in most situations). Tobler’s hiking function has been rescaled by a factor of 0.25 to represent the walking pace of a flock instead of humans. The factor corresponds to the ratio between the flock average speed (1.5 km/h) and the maximum human walking speed (about 6.0 km/h) on a favourable slope.
Garmy, Kaddouri, Rozenblat, and Schneider's hiking function (speed in kmh):
4 * exp(-0.008 * ((atan(abs(x[adj]))*180/pi)^2))
slope in degrees;
see: Herzog, I. (2020). Spatial Analysis Based on Cost Functions. In Gillings M, Haciguzeller P, Lock G (eds), "Archaeological Spatial Analysis. A Methodological Guide.",
Routledge: New York, 333-358 (with previous references).
Rees' hiking function (speed in kmh):
(1 / (0.75 + 0.09 * abs(x[adj]) + 14.6 * (abs(x[adj]))^2)) * 3.6
Rees' slope-dependant cost function; it is originally expressed in terms of time (1/v in Rees' publication);
here it is the reciprocal of time (i.e. speed) that is used in order to eventually get the reciprocal of speed (i.e. time).
Slope is dealt with here as originally expressed in Rees' publication (i.e. rise over run). The speed, which is originally expressed in m/s,
has been here transposed to kmh (i.e., multiplied by 3.6) for consistency with other hiking functions.
For this cost function see: Rees, WG (2004). Least-cost paths in mountainous terrain. Computers & Geosciences, 30(3), 203–209. See also: Campbell MJ, Dennison PE, Butler BW, Page WG (2019). Using crowdsourced fitness tracker data to model the relationship between slope and travel rates. Applied Geography 106, 93-107 (with previous references).
Relative energetic expenditure cost function:
1 / (tan((atan(abs(x[adj]))*180/pi)*pi/180) / tan (1*pi/180))
slope-based cost function expressing change in potential energy expenditure;
see Conolly, J., & Lake, M. (2006). Geographic Information Systems in Archaeology. Cambridge: Cambridge University Press, p. 220;
see also Newhard, J. M. L., Levine, N. S., & Phebus, A. D. (2014). The development of integrated terrestrial and marine pathways in the Argo-Saronic region, Greece. Cartography and Geographic Information Science, 41(4), 379-390, with references to studies that use this
function; see also ten Bruggencate, R. E., Stup, J. P., Milne, S. B., Stenton, D. R., Park, R. W., & Fayek, M. (2016). A human-centered GIS approach to modeling mobility on southern Baffin Island, Nunavut,
Canada. Journal of Field Archaeology, 41(6), 684-698. https://doi.org/10.1080/00934690.2016.1234897.
Herzog's metabolic cost function in J/(kg*m):
1 / ((1337.8 * abs(x[adj])^6) + (278.19 * abs(x[adj])^5) - (517.39 * abs(x[adj])^4) - (78.199 * abs(x[adj])^3) + (93.419 * abs(x[adj])^2) + (19.825 * abs(x[adj])) + 1.64)
see Herzog, I. (2016). Potential and Limits of Optimal Path Analysis. In A. Bevan & M. Lake (Eds.), Computational Approaches to Archaeological Spaces (pp. 179-211). New York: Routledge.
Wheeled-vehicle critical slope cost function:
1 / (1 + ((abs(x[adj])*100) / sl.crit)^2)
where sl.crit (=critical slope, in percent) is "the transition where switchbacks become more effective than direct uphill or downhill paths" and typically is in the range 8-16;
see Herzog, I. (2016). Potential and Limits of Optimal Path Analysis. In A. Bevan & M. Lake (Eds.), Computational Approaches to Archaeological Spaces (pp. 179-211). New York: Routledge.
Pandolf et al.'s metabolic energy expenditure cost function (in Watts):
1 / (1.5 * W + 2.0 * (W + L) * (L / W)^2 + N * (W + L) * (1.5 * V^2 + 0.35 * V * (abs(x[adj])*100)))
where W is the walker's body weight (Kg), L is the carried load (in Kg), V is the velocity in m/s, N is a coefficient representing ease of movement on the terrain (see above).
For this cost function, see Pandolf, K. B., Givoni, B., & Goldman, R. F. (1977). Predicting energy expenditure with loads while standing or walking very slowly. Journal of Applied Physiology,
43(4), 577-581. https://doi.org/10.1152/jappl.1918.104.22.1687.
For the use of this cost function in a case study, see Rademaker, K., Reid, D. A., & Bromley, G. R. M. (2012). Connecting the Dots: Least Cost Analysis, Paleogeography, and
the Search for Paleoindian Sites in Southern Highland Peru. In D. A. White & S. L. Surface-Evans (Eds.), Least Cost Analysis of Social Landscapes. Archaeological Case Studies (pp. 32-45).
University of Utah Press;
see also Herzog, I. (2013). Least-cost Paths - Some Methodological Issues, Internet Archaeology 36 (http://intarch.ac.uk/journal/issue36/index.html) with references.
Note: in the returned charts, the cost is transposed from Watts to Megawatts (see, e.g., Rademaker et al 2012 cited above).
Van Leusen's metabolic energy expenditure cost function (in Watts):
1 / (1.5 * W + 2.0 * (W + L) * (L / W)^2 + N * (W + L) * (1.5 * V^2 + 0.35 * V * abs(x[adj])*100) + 10))
which modifies the Pandolf et al.'s equation; see Van Leusen, P. M. (2002). Pattern to process: methodological investigations into the formation and interpretation of spatial patterns in archaeological landscapes. University of Groningen.
Note that, as per Herzog, I. (2013). Least-cost Paths - Some Methodological Issues, Internet Archaeology 36 (http://intarch.ac.uk/journal/issue36/index.html) and unlike Van Leusen (2002), in the above equation slope is expressed in percent and speed in m/s; also, in the last bit of the equantion, 10 replaces the value of 6 used by Van Leusen (as per Herzog 2013).
Note: in the returned charts, the cost is transposed from Watts to Megawatts.
Llobera-Sluckin's metabolic energy expenditure cost function (in KJ/m):
1 / (2.635 + (17.37 * abs(x[adj])) + (42.37 * abs(x[adj])^2) - (21.43 * abs(x[adj])^3) + (14.93 * abs(x[adj])^4))
for which see Llobera M.-Sluckin T.J. (2007). Zigzagging: Theoretical insights on climbing strategies, Journal of Theoretical Biology 249, 206-217.
Bellavia's cost function:
1 / (N * ((atan(abs(x[adj]))*180/pi)+1))
proposed by G. Bellavia, it measures abstract cost. Slope in degrees; N is a terrain factor (see above).
See: Herzog, I. (2020). Spatial Analysis Based on Cost Functions. In Gillings M, Haciguzeller P, Lock G (eds), "Archaeological Spatial Analysis. A Methodological Guide.",
Routledge: New York, 333-358 (with previous references).
Note that the walking-speed-related cost functions listed above are used as they are, while the other functions are reciprocated.
This is done since "gdistance works with conductivity rather than the more usual approach using costs"; therefore
"we need inverse cost functions" (Nakoinz-Knitter (2016). "Modelling Human Behaviour in Landscapes". New York: Springer, p. 183).
As a consequence, if we want to estimate time, we have to use the walking-speed functions as they are since the final accumulated values will correspond to the
reciprocal of speed, i.e. pace. In the other cases, we have to use 1/cost-function to eventually get cost-function/1.
The function returns a list storing the following components
dtm: Digital Terrain Model ('RasterLayer' class)
cost.surface: raster representing the cost-surface ('RasterLayer' class)
accumulated.cost.raster: raster representing the accumualted cost ('RasterLayer' class)
isolines: contour lines derived from the accumulated cost surface ('SpatialLinesDataFrame' class)
LCPs: estimated least-cost paths ('SpatialLines' class)
LCPs.back: estimated least-cost paths back to the origin ('SpatialLines' class)
LCPs$length: length of each least-cost path (units depend on the unit used in the input DTM)
LCPs.back$length: length of each least-cost path back to the origin (units depend on the unit used in the input DTM)
dest.loc.w.cost: copy of the input destination location(s) dataset with a new variable ('cost') added; if the cost is expressed in terms of time, the 'cost' variable will store the time values in decimal numbers, while another variable named 'cost_hms' will store the time values in sexagesimmal numbers (hours, minutes, seconds)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34
# load a sample Digital Terrain Model volc <- raster::raster(system.file("external/maungawhau.grd", package="gdistance")) # load a sample start location on the above DTM data(volc.loc) # load the sample destination locations on the above DTM data(destin.loc) # calculate walking-time isochrones based on the on-path Tobler's hiking function (default), # setting the time unit to hours and the isochrones interval to 0.05 hour; # also, since destination locations are provided, # least-cost paths from the origin to the destination locations will be calculated # and plotted; 8-directions move is used result <- movecost(dtm=volc, origin=volc.loc, destin=destin.loc, move=8, breaks=0.05) # same as above, but using the Irmischer-Clarke's hiking function (male, on-path) result <- movecost(dtm=volc, origin=volc.loc, destin=destin.loc, funct="icmonp", move=8, breaks=0.05) # same as above, but using the 'cognitive slope' result <- movecost(dtm=volc, origin=volc.loc, destin=destin.loc, funct="icmonp", move=8, breaks=0.05, cogn.slp=TRUE) # calculate accumulated cost surface and the least-cost path between the # origin and one destination, and also calculate the LCP back to the origin results <- movecost(dtm=volc, origin=volc.loc, destin=destin.loc[2,], move=8, return.base = TRUE)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.