Markdown document to describe the LCP generation process for a single path

This workflow is for non-rescaled visibility which is weighted at a 50/50 ratio with speed. It is only one path in one direction. Results show descriptive statistics calculated for the generated path and a plot of the path in the study area.

1 - Make sure script and necessary terrain rasters are in designated folder
2 - Load packages

library(gdistance)
library(sf)
library(terra)
library(viridis)

3 - Read and name terrain rasters AS RasterLayer

dtm = raster("dtm.tif")
rgh = raster("roughness.tif")
veg = raster("veg.tif")
vis = raster("visibility.tif")

4 - Define transition functions

tr_diff = function(x){x[2] - x[1]}
tr_mean = function(x){(x[1] + x[2]) / 2}

5 - Generate landscape transition matrices

slo_tr = transition(dtm, tr_diff, directions=16, symm=F)
veg_tr = transition(veg, tr_mean, directions=16)
rgh_tr = transition(rgh, tr_mean, directions=16)
vis_tr = transition(vis, tr_mean, directions=16)

6 - Convert slope from radians -> degrees

slo_tr = atan(slo_tr) * 180/pi

7 - Create adjacency table

adj = adjacent(dtm, cells = 1:ncell(dtm), pairs = TRUE, directions = 16)

8 - Define travel coefficients

a = 0.408
b = 20.953
c = 71.267
d = 0.649
e = -0.008
f = 1.181
g = 6.992

9 - Create cost/conductance matrix using travel rate function and weights of visibility and speed at 50/50

w1=0.5
w2=0.5
all_tr = slo_tr
all_tr[adj] = w1 * ((c * (1/ (pi * b * (1 + ((slo_tr[adj] - a)/b) ^ 2))) + d + e * slo_tr[adj])/ 
                        (f * veg_tr[adj] + g * rgh_tr[adj] + 1)) + w2 * vis_tr[adj]
conductance = geoCorrection(all_tr)

10 - Generate LCPs

#define extent of study area
width <- diff(extent(dtm)[1:2])
height <- diff(extent(dtm)[3:4])
#define start and end points
A = c(extent(dtm)[1] + 0.2 * width, extent(dtm)[4] - 0.1 * height)
B = c(extent(dtm)[2] - 0.1 * width, extent(dtm)[3] + 0.1 * height)
#generate LCP
AtoB <- shortestPath(conductance, A, B, output = "SpatialLines")

11 - Generate metrics

#convert SpatialLines to sf
AtoB <- st_as_sf(AtoB)
#path length
AtoB_len <- st_length(AtoB) |> units::drop_units()
#euclidean distance
AtoB_dst <- sqrt((B[1]-A[1])^2 + (B[2]-A[2])^2)
#sinuosity
AtoB_sin <- AtoB_len / AtoB_dst
#average visibility
vis_rast <- rast(vis)
AtoB_sv <- vect(AtoB)
AtoB_vis_mean <- terra::extract(vis_rast, AtoB_sv, mean)[1,2]
#percent "high visibility" 
vis_rc <- ifel(vis_rast > 0.25, 1, 0)
AtoB_perc_highvis <- terra::extract(vis_rc, AtoB_sv, mean)[1,2]
#average roughness
rgh_rast = rast(rgh)
AtoB_rgh_mean = terra::extract(rgh_rast, AtoB_sv, mean)[1,2]
#average veg dens
veg_rast = rast(veg)
AtoB_veg_mean = terra::extract(veg_rast, AtoB_sv, mean)[1,2]
#save metrics to data frame
df <- data.frame(pathLength = AtoB_len, 
                 startEndDistance = AtoB_dst,
                 sinuosity = AtoB_sin,
                 percentHighVis = AtoB_perc_highvis,
                 avgVis = AtoB_vis_mean,
                 avgRgh = AtoB_rgh_mean,
                 avgVeg = AtoB_veg_mean)

12 - Plot

linecol = rgb(1,0,0.6)
plot(dtm, col = grey(1:100/100), main="DTM", axes = FALSE)
lines(AtoB, col = linecol, lwd = 2)
legend("bottomright", 
       legend="vis weight = 0.5, speed weight = 0.5", 
       title = "LCP", 
       cex = 0.5, 
       lwd = 1,
       col = linecol,
       bty = "o")

print(df, digits=2)
##   pathLength startEndDistance sinuosity percentHighVis avgVis avgRgh avgVeg
## 1      12185             9751       1.2           0.28   0.18  0.035   0.12

Discussion:
Next steps will be to generate paths weighting visibility and speed at different ratios, creating rescaled visibility matrices, and repeating paths at different start and end points in a range of locations and distances apart. Additionally, travel coefficients and the travel rate function is under review and subject to change.