1. Introducing leastcostpath

This vignette describes leastcostpath, a package written for use in the R environment (R Core Team, 2018) and built on classes and functions provided in the R package gdistance (Van Etten, 2017). leastcostpath provides the functionality to calculate Least Cost Paths (LCPs) using numerous cost functions (time- and energy-based) that approximate the difficulty of moving across a landscapee. Furthermore, additional cost surfaces can be incorporated into the analysis via create_barrier_cs() or create_feature_cs().

Secondly, leastcostpath provides the functionality to calculate Stochastic Least Cost Paths (Pinto and Keitt, 2009); and Probabilistic Least Cost Paths (Lewis, 2020).

Thirdly, leastcostpath provides functionality to calculate movement potential within a landscape through the implementation of From-Everywhere-to-Everywhere (White and Barber, 2012); Cumulative Cost Paths (Verhagen 2013), and Least Cost Path calculation within specified distance bands (Llobera, 2015).

Lastly, leastcostpath provides the functionality to validate the accuracy of the computed Least Cost Path relative to another path via validate_lcp() (Goodchild and Hunter, 1997) and PDI_validation() (Jan et al. 1999).

2. using leastcostpath

2.1 Setup

library(rgdal)
library(rgeos)
library(sp)
library(raster)
library(gdistance)
library(leastcostpath)
dem <- raster::raster(system.file('external/maungawhau.grd', package = 'gdistance'))

dem_extent <- as(raster::extent(dem), 'SpatialPolygons')

raster::crs(dem_extent) <- raster::crs(dem)

raster::plot(dem)
raster::plot(dem_extent, add = T, border = "red")

2.2 Creating Slope-based Cost Surfaces

# cost functions currently implemented within leastcostpath

cfs <- c("tobler", "tobler offpath", "irmischer-clarke male", "irmischer-clarke offpath male", "irmischer-clarke female", "irmischer-clarke offpath female",
        "modified tobler", "wheeled transport", "herzog", "llobera-sluckin", "campbell 2019")

# neighbours can be 4, 8, 16, 32, or 48. 8 used for illustration purposes, but greater the number the better the cost surface / LCP approximates reality. 
neigh <- 4

slope_cs <- leastcostpath::create_slope_cs(dem = dem, cost_function = "tobler", neighbours = neigh)

plot(raster(slope_cs), col = grey.colors(100))

2.3 Creating and Incorporating Barriers within Cost Surfaces

# create barrier of altitude values above 180

altitude <- dem > 180
# convert value 0 to NA - this ensures that the areas below 180 are not viewed as barriers. That is, if the raster value is NA then create_barrier_cs() will assume that the cell is background and will therefore assign the background argument value 
altitude[altitude == 0] <- NA

# the values NOT NA will be assigned the field argument value. If 0 (default) then movement within the area will be completely inhibited.
altitude_cs <- leastcostpath::create_barrier_cs(raster = altitude, barrier = altitude, neighbours = neigh, field = 0, background = 1)

# multiplying the two cost surfaces ensures that barriers continue to completely inhibit movement (i.e. slope_cs values * 0 = 0)
slope_altitude_cs <- slope_cs * altitude_cs

plot(raster(slope_altitude_cs), col = grey.colors(100))

2.3 Calculating Least Cost Paths

A <- sp::SpatialPoints(cbind(2667930, 6479466))
B <- sp::SpatialPoints(cbind(2667668, 6478818))

# if cost function is anisotropic (e.g. Tobler's Hiking Function) then LCP from A-to-B and B-to-A may be different. To create LCP in both directions set the directional argument to FALSE (default)
lcp <- leastcostpath::create_lcp(cost_surface = slope_altitude_cs, origin = A, destination = B, directional = FALSE)

plot(raster(slope_altitude_cs), col = grey.colors(100))
plot(A, add = T, col = "black")
plot(B, add = T, col = "black")
# LCP from A-to-B
plot(lcp[1,], add = T, col = "red")
# LCP from B-to-A
plot(lcp[2,], add = T, col = "blue")

2.4 Calculating Cost Corridors

# calculates accumulated cost surface from A and from B. These are then averaged. rescale argument TRUE (not default) rescales the accumulated cost surface to 1
cc <- leastcostpath::create_cost_corridor(cost_surface = slope_altitude_cs, origin = A, destination = B, rescale = TRUE)

plot(cc, col = heat.colors(100))
plot(A, add = T, col = "white", pch = 16)
plot(B, add = T, col = "white", pch = 16)

2.5 Calculating From-Everywhere-To-Everywhere Least Cost Paths

random_locs <- sp::spsample(x = dem_extent, n = 25, type = 'random')

fete_lcp <- leastcostpath::create_FETE_lcps(cost_surface = slope_altitude_cs, locations = random_locs, cost_distance = FALSE, parallel = FALSE)

plot(raster(slope_altitude_cs), col = grey.colors(100))
plot(fete_lcp, add = T, col = "red")
plot(random_locs, add = T, col = "black", pch = 16)

2.6 Least Cost Path Density

# rasterises the LCPs and calculates the number of times an LCP crosses a raster cell. If rescale argument TRUE (not default) then LCP density scaled to 1. If rasterize_as_points argument TRUE (default) then LCP vertices points are rasterized. This is quicker than rasterizing lines (rasterize_as_lines argument FALSE) but will contain 'gaps' where no LCP vertex point is present 
fete_lcp_density <- leastcostpath::create_lcp_density(lcps = fete_lcp, raster = dem, rescale = TRUE, rasterize_as_points = TRUE)

fete_lcp_density[fete_lcp_density == 0] <- NA

plot(fete_lcp_density, col = heat.colors(100))

2.7 Probabilistic Least Cost Paths

# least cost paths are calculated using input data that contains errors. The most notable is the Digital Elevation Model (DEM) which is used to model movement difficulty within the landscape. By incorporating DEM error (in this case vertical error), the error is propagated throughout the LCP modelling process. This results in the input data better approximating the real world, with the LCP being more reliable as a result  

# for illustration purposes, we will assume that the vertical error (root-mean-square-error) is +/- 0.5m.
RMSE <- 0.5
n <- 100
lcps <- list()

for (i in 1:n) {
  
  print(i)

lcps[[i]] <- leastcostpath::create_lcp(cost_surface = leastcostpath::create_slope_cs(dem = leastcostpath::add_dem_error(dem = dem, rmse = RMSE, type = "autocorrelated"), cost_function = "tobler", neighbours = neigh) * altitude_cs, origin = A, destination = B, directional = FALSE)

}
## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] 6
## [1] 7
## [1] 8
## [1] 9
## [1] 10
## [1] 11
## [1] 12
## [1] 13
## [1] 14
## [1] 15
## [1] 16
## [1] 17
## [1] 18
## [1] 19
## [1] 20
## [1] 21
## [1] 22
## [1] 23
## [1] 24
## [1] 25
## [1] 26
## [1] 27
## [1] 28
## [1] 29
## [1] 30
## [1] 31
## [1] 32
## [1] 33
## [1] 34
## [1] 35
## [1] 36
## [1] 37
## [1] 38
## [1] 39
## [1] 40
## [1] 41
## [1] 42
## [1] 43
## [1] 44
## [1] 45
## [1] 46
## [1] 47
## [1] 48
## [1] 49
## [1] 50
## [1] 51
## [1] 52
## [1] 53
## [1] 54
## [1] 55
## [1] 56
## [1] 57
## [1] 58
## [1] 59
## [1] 60
## [1] 61
## [1] 62
## [1] 63
## [1] 64
## [1] 65
## [1] 66
## [1] 67
## [1] 68
## [1] 69
## [1] 70
## [1] 71
## [1] 72
## [1] 73
## [1] 74
## [1] 75
## [1] 76
## [1] 77
## [1] 78
## [1] 79
## [1] 80
## [1] 81
## [1] 82
## [1] 83
## [1] 84
## [1] 85
## [1] 86
## [1] 87
## [1] 88
## [1] 89
## [1] 90
## [1] 91
## [1] 92
## [1] 93
## [1] 94
## [1] 95
## [1] 96
## [1] 97
## [1] 98
## [1] 99
## [1] 100
lcps <- do.call(rbind, lcps)

lcps_A_B <- leastcostpath::create_lcp_density(lcps = lcps[lcps$direction == "A to B",], raster = dem, rescale = TRUE, rasterize_as_points = TRUE)

lcps_A_B[lcps_A_B == 0] <- NA

lcps_A_B <- (lcps_A_B/n)

lcps_B_A <- leastcostpath::create_lcp_density(lcps = lcps[lcps$direction == "B to A",], raster = dem, rescale = FALSE, rasterize_as_points = TRUE)

lcps_B_A[lcps_B_A == 0] <- NA

lcps_B_A <- (lcps_B_A/n)

plot(raster(slope_altitude_cs), legend = FALSE, col = grey.colors(100))
plot(lcps_A_B, add = T, col = heat.colors(100))
plot(A, add = T, col = "white", pch = 16)
plot(B, add = T, col = "white", pch = 16)

plot(raster(slope_altitude_cs), legend = FALSE, col = grey.colors(100))
plot(lcps_B_A, add = T, col = heat.colors(100))
plot(A, add = T, col = "white", pch = 16)
plot(B, add = T, col = "white", pch = 16)

References

Goodchild, M. F., & Hunter, G. J. (1997). A simple positional accuracy measure for linear features. International Journal of Geographical Information Science, 11(3), 299–306. https://doi.org/10.1080/136588197242419

Jan, O., Horowitz, A., & Peng, Z. (1999). Using GPS Data to Understand Variations in Path Choice.

Lewis, J. (2020). Probabilistic Modelling using Monte Carlo Simulation for Incorporating Uncertainty in Least Cost Path Results: a Postdictive Roman Road Case Study (preprint). SocArXiv. https://doi.org/10.31235/osf.io/mxas2

Llobera, M. (2015). Working the digital: some thoughts from landscape archaeology. Material evidence: learning from archaeological practice. Routledge, Abingdon, 173–188.

Pinto, N., & Keitt, T. H. (2009). Beyond the least-cost path: evaluating corridor redundancy using a graph-theoretic approach. Landscape Ecology, 24(2), 253–266. https://doi.org/10.1007/s10980-008-9303-y

R Core Team. (2018). R: A language and environment for statistical computing. Viena, Austria: R Foundation for Statistical Computing. https://www.R-project.org/

van Etten, J. (2017). R Package gdistance: Distances and Routes on Geographical Grids. Journal of Statistical Software, 76(13). https://doi.org/10.18637/jss.v076.i13

Verhagen, P. (2013). On the Road to Nowhere?: Least Cost Paths, Accessibility and the Predictive Modelling Perspective (pp. 383–389). Presented at the CAA2010: fusion of cultures: proceedings of the 38th Annual Conference on Computer Applications and Quantitative Methods in Archaeology, Granada, Spain, April 2010, Archaeopress.

White, D., & Barber, S. L. (2012). Geospatial modeling of pedestrian transportation networks: a case study from precolumbian Oaxaca, Mexico. Journal of Archaeological Science, 39(8), 2684–2696. https://doi.org/10.1016/j.jas.2012.04.017