Provides the facility to calculate non-isotropic accumulated cost surface, least-cost paths, and least-cost corridors using a number of human-movement-related cost functions that can be selected by the user. It just requires a Digital Terrain Model, a start location and (optionally) destination locations. See Alberti (2019)
movecost provides the facility to calculate accumulated cost surface and least-cost paths using a number of human-movement-related cost functions that can be selected by the user. It just requires a Digital Terrain model, a start location and (optionally) destination locations.
The package comes with some toy datasets:
destin.loc: SpatialPointsDataFrame representing spots on the volcano Maunga Whau (Auckland, New Zealand), to be used as destination locations for least-cost paths calculation.
volc.loc: SpatialPointsDataFrame representing a spot on the volcano Maunga Whau (Auckland, New Zealand).
movecost(): function for calculating accumulated cost of movement across the terrain and least-cost paths from an origin.
moveCost(): provides the facility to calculate the accumulated cost of movement around a starting location and to optionally calculate least-cost paths 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.
If the parameter
destin is fed with a dataset representing destination location(s) (
SpatialPointsDataFrame class), the function will also calculate least-cost path(s) plotted on the input DTM; the length of each path will be saved under the variable
length stored in the
LCPs dataset (
SpatialLines class) returned by the function. The red dot(s) representing the destination location(s) will be labelled with numeric values representing the cost value at the location(s). The cost value will be also appended to the updated destination dataset returned by the function and storing a new variable named
The function builds on functions out of Jacob van Etten's
gdistance package, and by default uses a 16-directions movement in calculating the accumulated cost-surface. The number of movements can be optionally set by the user via the
The slope is internally calculated by the function using the
terrain() function out of the
raster package, using 8 neighbors. The slope is calculated in degrees and, for the sake of its use within the implemented cost functions, it is internally modified (i.e., turned into either gradient or percent) according to the user-selected cost function. The user can input a slope dataset (
RasterLayer class) using the parameter
slope; this can prove usefull in cases when the slope has to be preliminarily modified, for instance to mask out areas that cannot be crossed or to weight the slope values according to a user-defined weighting scheme.
The following cost functions are implemented (x stands for slope):
6 * exp(-3.5 * abs(tan(x*pi/180) + 0.05))
(6 * exp(-3.5 * abs(tan(x*pi/180) + 0.05))) * 0.6
as per Tobler's indication, the off-path walking speed is reduced by 0.6.
4.8 * exp(-5.3 * abs((tan(x*pi/180) * 0.7) + 0.03))
modified version of the Tobler's hiking function as proposed by Joaquín Márquez-Pérez, Ismael Vallejo-Villalta & José I. Álvarez-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.
(0.11 + exp(-(tan(x*pi/180)*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: the function originally expresses speed in m/s; it has been is reshaped (multiplied by 3.6) to turn it into kmh for consistency with the other Tobler-related cost functions; also, originally the slope is in percent;
tan(x*pi/180)*100 turns the slope from degrees to percent (=rise/run*100).
(0.11 + 0.67 * exp(-(tan(x*pi/180)*100 + 2)^2 / (2 * 30)^2)) * 3.6
1/ (0.0277 * (tan(x*pi/180)*100) + 0.6115)
proposed by Uriarte González; see: Chapa Brunet, T., García, J., Mayoral Herrera, V., & Uriarte González, 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. Also, originally the slope is in percent:
tan(x*pi/180)*100 turns the slope from degrees to percent (=rise/run*100). 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.
1 / (tan(x*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.
1 / ((1337.8 * tan(x*pi/180)^6 + 278.19 * tan(x*pi/180)^5 - 517.39 * tan(x*pi/180)^4 - 78.199 * tan(x*pi/180)^3 + 93.419 * tan(x*pi/180)^2 + 19.825 * tan(x*pi/180) + 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.
1 / (1 + ((tan(x*pi/180)*100) / sl.crit)^2)
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.
1 / (1.5 * W + 2.0 * (W + L) * (L / W)^2 + N * (W + L) * (1.5 * V^2 + 0.35 * V * (tan(x*pi/180)*100)))
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. As for the latter, suggested values available in literature are:
Ploughed field=1.3; see de Gruchy, M., Caswell, E., & Edwards, J. (2017). Velocity-Based Terrain Coefficients for Time-Based Models of Human Movement. Internet Archaeology, 45(45). https://doi.org/10.11141/ia.45.4.
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.19184.108.40.2067. 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).
1 / (1.5 * W + 2.0 * (W + L) * (L / W)^2 + N * (W + L) * (1.5 * V^2 + 0.35 * V * (tan(x*pi/180)*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.
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.
When using the Tobler-related cost functions, the time unit can be selected by the user setting the
time parameter to
h (hour) or to
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 image 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.
The function returns a list storing:
$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$length: length of each least-cost path (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.
first release to CRAN.