Purpose & Requirements

This script generates least cost paths (LCPs) and their descriptive metrics.

This step is part of a larger project on using LCPs to effectively model wildland firefighter escape routes through the use of landscape metrics including slope, vegetation density, and terrain roughness. This novel approach includes the use of a modeled visibility index in an attempt to estimate how visibility impacts situational awareness for firefighters along the escape routes and ultimately improve safety in the field.

The workflow is completed entirely in R, and required packages include raster, terra, and sf for spatial modeling, and gdistance by van Etten et al (2023) for the LCP generation. Landscape metrics in the format of rasters are the result of lidar data courtesy of USGS 3DEP. The travel cost or escape route algorithm description is provided by Mickey Campbell.

This script was created as part of my master’s thesis research requirements through the University of Utah Department of Geography. By visiting the Utah Remote Sensing Applications (URSA) website, you can find more information on my research lab’s current projects, published papers, and collaborators.


Workflow Summary


Define Workspace and Packages

Set the working directory for wherever scripts and data will be saved.
Load packages including raster and terra due to gdistance’s use of raster; the “older version” of the updated terra package. We also use sf when computing LCP metrics, and gdistance for the transition matrices, cost matrix, and other LCP generation requirement steps.

library(terra)
library(raster)
library(sf)
library(gdistance)

Define and Check Landscape Rasters

Landscape rasters include DTM, vegetation density, terrain roughness, and visibility. These rasters are the result of a separate workflow described here.

Check and plot the landscape rasters we defined.

dtm = raster(dtm)
plot(dtm)

veg = raster(veg)
plot(veg)

rgh = raster(rgh)
plot(rgh)

vis = raster(vis)
plot(vis)

Clipping raster extent

To avoid edge effects which have been shown to result in NA values when generating LCPs, we clip down the extent of our area of interest. This still gives us a large area to analyze, while avoiding those NA values.

# current raster extents: 
#      min(x) = 417000, max(x) = 423000, min(y) = 4450000, max(y) = 4461000
# setting 'e' to equal desired extent, removing 1000 meter buffer 
e = as(extent(418000, 422000, 4451000, 4460000), "Spatial Polygons")
# matching coordinate reference system
crs(e) = crs(vis)
# apply the buffer (aka 'crop') to all the landscape rasters
vis = crop(vis, e, 'in', overwrite = T)
dtm = crop(dtm, e, 'in', overwrite = T)
veg = crop(veg, e, 'in', overwrite = T)
rgh = crop(rgh, e, 'in', overwrite = T)

Define Transition Functions and Create Matrices

First, we define our transition functions which will help calculate the cost of movement between nodes in our transition matrices. Difference and mean are the two functions used in this workflow.

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

Next, we create the transition matrices from the landscape metrics rasters.

veg_tr = transition(veg, tr_mean, directions=8)
ter_tr = transition(rgh, tr_mean, directions=8)
slo_tr = transition(dtm, tr_diff, directions=8, symm=F)

Cells are connected using 8 directions, both orthagonally and diagonally through the Moore neighborhood, resulting in a king’s graph matrix. This is a commonly used method used in GIS software.

We also convert our slope transition matrix from radians to degrees.

slo_tr = atan(slo_tr) * 180/pi

Define Visibility Rescaling Method and Factor

To further explore the effect that visibility has on selecting a safe path, we will create multiple paths with multiple values of visibility. Not only do we change the weights of importance of a non-rescaled visibility-to-speed iteration, we also conduct logistically, logarithmically, and exponentially rescaled visibility values tests. For an additional explanation of this process including algorithmic information, see here.

First we define the rescaling method; we create a list of values, with 0 resulting in a logarithmic rescale, 0.5 resulting in logistic growth rescaling, and 1.0 resulting in exponentially rescaled visibility. In addition to the documentation above, see here for more information on rescaling visibility for this application.

rescale_method = c(0, 0.5, 1.0)

Next, the rescaling factor defines the magnitude of the visibility rescaling operation.

rescale_factor = c(0, 5, 10, 20)

Define Cost Matrix Variables

A few more steps are needed to generate our final cost matrix. For a more in-depth explanation of these steps, please view the gdistance vignette.

Create adjacency table

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

Define travel rate coefficients

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

Define any other serial needs for nested loop

# boundary extent for selecting start and end points
bdy = extent(dtm)
# start counter for progress and time keeping
counter = 0
# set seed to be able to reference previous start and end points
set.seed(1)

LCP Generation + Calculation Loop

Now we jump into our LCP Generation Loop.
This is a four-part nested loop to satisfy our variables for LCP generation.
* The outer loop is to establish start and end points; we choose 100 start and end pairs for this process.
* The second loop in is to call our rescaling methods which we established above.
* The third loop is to call our rescaling factors which we also established above.
* Our fourth and final inner loop establishes weights to apply to the values of maximizing visibility versus maximizing travel rate.

When we create a cost matrix which the LCPs are based on, we use a travel rate modeling function which estimates travel rate based on our landscape metrics. In this project, we analyze the tradeoffs of maximizing visibility versus travel rate, so by applying weights from 0 to 1 by 0.1 increments as metrics of “importance”, we can analyze how the LCPs change course. View the equation below which we use in the LCP loop:

\[ cost = \left[w_1(\frac{(c(\frac{1}{(\pi b (1 + (\frac{(slope - a)}{b})^{2}))}) + d + e(slope))}{f(density) + g(roughness) + 1})\right] + \left[w_2(visibility)\right] \]

The first bracket which contains weight 1 (\(w_1\)) defines the travel rate equation, and the second bracket with \(w_2\) defines visibility.

In the remainder of the loop after generating each individual LCP, we also calculate the LCP’s metrics, including path length, straight-line distance from the start to end points, sinuosity, average values of visbility, roughness, and vegetation density, the percent of the path that is determined to have high visibility (above 0.25), and travel time along the LCP. We save these values and the LCP shape to a data frame and write to file in the next step.

#set start timer
t1 = Sys.time()
#start and end points
for (i in seq(1,100)){
  pt1x <- runif(1, bdy[1], bdy[2])
  pt1y <- runif(1, bdy[3], bdy[4])
  pt2x <- runif(1, bdy[1], bdy[2])
  pt2y <- runif(1, bdy[3], bdy[4])
#rescaling methods
  for (v0 in rescale_method){
#rescaling factors
    for (k in rescale_factor){
      log.fun <- function(v, v0, L, k){
        vr <- L/ (1 + exp(-k * (v - v0)))
        return(vr)
      }
      vis.rs <- log.fun(vis, v0, 1, k)
      vis_tr = transition(vis.rs, tr_mean, directions=8)
#weights
      for (w1 in seq(0,1,0.1)){
        w2 = 1 - w1
        all_tr = veg_tr
#make cost transition matrix
        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 * ter_tr[adj] + 1))) + (w2 * vis_tr[adj])
        conductance = geoCorrection(all_tr)
        start = c(pt1x, pt1y)
        finish = c(pt2x, pt2y)
        lcp = gdistance::shortestPath(conductance, start, finish, output = "SpatialLines")
        
#generate metrics
        #convert line from SpatialLine to simple features object 
        lcp_sf <- st_as_sf(lcp)
        #path length
        lcp_len <- st_length(lcp_sf) |> units::drop_units()
        #euclidean distance
        lcp_dst <- sqrt((pt2x-pt1x)^2 + (pt2y-pt1y)^2)
        #sinuosity
        lcp_sin <- lcp_len / lcp_dst
        #average visibility
        vis_rast <- rast(vis)
        lcp_sv <- vect(lcp_sf)
        lcp_vis_mean <- terra::extract(vis_rast, lcp_sv, mean)[1,2]
        #percent "high visibility" 
        vis_rc <- ifel(vis_rast > 0.25, 1, 0)
        lcp_perc_highvis <- terra::extract(vis_rc, lcp_sv, mean)[1,2]
        #average roughness
        rgh_rast = rast(rgh)
        lcp_rgh_mean = terra::extract(rgh_rast, lcp_sv, mean)[1,2]
        #average veg dens
        veg_rast = rast(veg)
        lcp_veg_mean = terra::extract(veg_rast, lcp_sv, mean)[1,2]
        
        #travel time
        # creating copy of dtm for resolution/extent/etc.
        AtoB_rast <- dtm 
        # setting all of the values of the new raster to 1
        values(AtoB_rast) <- 1 
        # setting all pixels outside of the AtoB path to zero
        AtoB_rast <- mask(AtoB_rast, lcp, updatevalue = 0) 
        # creating a transition matrix where all cell-to-cell movement *on* the path has conductance of 1, and all movement *off* the path has conductance of 0 (acts as barrier), forcing movement to be on the path
        AtoB_tr <- transition(AtoB_rast, min, 8) 
        # creating a copy of the transition matrix specifically to calculate speed
        speed_tr <- slo_tr 
        # creates speed transition matrix, but only on the AtoB path
        speed_tr[adj] = ((c * (1/ (pi * b * (1 + ((slo_tr[adj] - a)/b) ^ 2))) + d + e * slo_tr[adj])/(f * veg_tr[adj] + g * ter_tr[adj] + 1)) * AtoB_tr[adj] 
        # geocorrection to adjust for angular distances
        speed_conductance = geoCorrection(speed_tr) 
        # calculate time along that path
        time_AtoB <- costDistance(speed_conductance, start, finish) 
        #end travel time calc
        
        #attach metrics as attributes to lcp
        lcp_sf$pt1x = pt1x
        lcp_sf$pt1y = pt1y
        lcp_sf$pt2x = pt2x
        lcp_sf$pt2y = pt2y
        lcp_sf$speedwt = w1
        lcp_sf$viswt = w2
        lcp_sf$rscmeth = v0
        lcp_sf$rscfact = k
        lcp_sf$lcpLen = lcp_len
        lcp_sf$eucDist = lcp_dst
        lcp_sf$sinuosity = lcp_sin
        lcp_sf$percHV = lcp_perc_highvis
        lcp_sf$avgVis = lcp_vis_mean
        lcp_sf$avgRgh = lcp_rgh_mean
        lcp_sf$avgVeg = lcp_veg_mean
        
        # summarize results to data frame
        if (i == 1 &
            v0 == 0 & 
            k == 0 &
            w1 == 0){
          lcps = lcp_sf
          df <- data.frame(pt1x = pt1x,
                           pt1y = pt1y,
                           pt2x = pt2x,
                           pt2y = pt2y,
                           v0 = v0,
                           k = k,
                           w1 = w1,
                           pathLength = lcp_len,
                           startEndDistance = lcp_dst,
                           sinuosity = lcp_sin,
                           percentHighVis = lcp_perc_highvis,
                           avgVis = lcp_vis_mean,
                           avgRgh = lcp_rgh_mean,
                           avgVeg = lcp_veg_mean,
                           traveltime = time_AtoB
                           )
        } else {
          lcps = rbind(lcps, lcp_sf)
          df <- rbind(df, data.frame(pt1x = pt1x,
                                     pt1y = pt1y,
                                     pt2x = pt2x,
                                     pt2y = pt2y,
                                     v0 = v0,
                                     k = k,
                                     w1 = w1,
                                     pathLength = lcp_len,
                                     startEndDistance = lcp_dst,
                                     sinuosity = lcp_sin,
                                     percentHighVis = lcp_perc_highvis,
                                     avgVis = lcp_vis_mean,
                                     avgRgh = lcp_rgh_mean,
                                     avgVeg = lcp_veg_mean,
                                     traveltime = time_AtoB
                                     ))
        }
        
        #calculate time of line creation + calculations
        t2 = Sys.time()
        time = t2-t1
        counter = counter + 1
        print(paste0("lines + metrics created: ", counter, ", time: ", time))
      }
    }
  }
}

Set system to run sequentially

plan(sequential)

Write Results to File and Plot

We make sure to remember to save our data frame generated from the loop above and the metrics as a csv, and the LCPs and their attributes as a shapefile.

write.csv(df, "lcp_metrics/lcp_metrics_100loop.csv")
st_write(lcps, "lcp_metrics/lcps_100loop.shp", overwrite = T)

Plot to see the resulting lines! Warning: there are a lot of them

plot(lcps)