Generating Landscape Rasters

---

Introduction

In this step of my thesis project, we create landscape metrics from the high resolution lidar downloaded from the USGS. These landscape metrics are saved as rasters for use later in generating and analyzing least cost paths.

Requirements include the packages listed below, and USGS 3DEP lidar data.

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.

Set the working directory and packages

setwd(wd)

library(lidR)
library(terra)
library(future)
library(RColorBrewer)

Create a catalog of lidar data

The lidR package is used frequently in this script to work with lidar files using the various tools and functions it provides. Please view the lidR book as a reference for the methods used in this process.

First we create a catalog of our lidar files downloaded from the USGS.

#name and read in catalog
ctg = readLAScatalog(lakemtns_lidar)
#plot catalog as tiles
plot(ctg)

Other catalogs are created in a few different sections in an attempt to preserve the original data downloaded. Those catalogs follow the same order of operations described above.

---

Generate DTM, DSM, CHM

In this section, commonly derived landscape models including a digital terrain model (DTM), digital surface model (DSM), and canopy height model (CHM) are constructed. I will share the code I used for these steps and plot the results.

Digital Terrain Model (DTM)

We use the default parameters to process our DTM using the catalog of lidar specified above. Because we use the catalog engine, everything is merged into a single raster output.

#start parallel processing
plan(multisession, workers = 5L)
#name and create dtm by tin method
dtm = lidR::rasterize_terrain(ctg, algorithm = tin())
plan(sequential)

plot(dtm)

Digital Surface Model (DSM)

Our DSM raster represents the highest elevation of ALS returns. We use the same TIN method as our DTM.

#start parallel processing
plan(multisession, workers = 5L)
#name and create dsm by dsmtin method
dsm = rasterize_canopy(ctg, algorithm = dsmtin())
plan(sequential)

plot(dsm)

Looks pretty much the same as the DTM.. good!

Canopy Height Model (CHM)

The lidR book includes a workflow for generating a CHM, but a faster method which satisfies our needs is to find the difference between the DSM and the DTM.

chm = dsm - dtm

plot(chm, breaks=cuts, col=pal(8))

Knowledge of the area combined with the CHM results lets me know that the vegetation in this area is low-lying and shrub-dominant. Therefore, canopy is generally no greater than 15m.

---

Generate Terrain Roughness and Vegetation Density

These landscape metrics are more novel and computationally expensive. See Campbell et al. (2017) for metrics reference.

Terrain Roughness

Terrain roughness is calculated as the absolute value of the difference between a fine-resolution DTM and a “smoothed” DTM.

Start parallel processing and create lidar, DTM, and working directories for the terrain roughness computation process.

#turn on parallel processing
plan(multisession, workers = 10L)

#define roughness lidar directory
roughness_dir_laz = paste0(roughness_dir, "\\rgh_lidar")
#create new directory to copy lidar files into 
dir.create(roughness_dir_laz, showWarnings = T, recursive = T)

#define roughness dtm directory
roughness_dir_dtm = paste0(roughness_dir, "\\dtms_0_25m")
#create new directory for individual DTMs
dir.create(roughnees_dir_dtm, showWarnings = T, recursive = T)

Copy lidar files into the new directory, create a list of the roughness lidar files, and read the files into a catalog

#create list of lidar tile files
lidar_files = list.files(lakemtns_lidar)
#copy list to new roughness lidar directory
file.copy(from = lidar_files, 
          to = paste0(roughness_dir_laz, lidar_files),
          overwrite = TRUE,
          recursive = TRUE)
#create new list of roughness lidar tile files
#specify only the .laz files (just in case), and the full file name
roughness_lidar = list.files(roughness_dir_laz, 
                             pattern = "*.laz", 
                             full.names = TRUE)
#create roughness lidar catalog
ctg_roughness = readLAScatalog(roughness_lidar)

plot(ctg_roughness)

Before performing operations on the data, perform processing efficiency operations

#prevent stopping early, ignoring errors
opt_stop_early(ctg_roughness) = F
#perform operations on disk
opt_output_files(ctg_roughness) = paste0(roughness_dir_dtm, "/tile_{XLEFT}_{YBOTTOM}")
#reduce chunk size to reduce processing speed
opt_chunk_size(ctg_roughness) = 250
#reduce chunk buffer to reduce processing speed
opt_chunk_buffer(ctg_roughness) = 5

Generate fine-resolution DTM, save, read DTM back in, and turn off parallel processing

#generate dtm at 0.25m resolution for each lidar tile
dtm_r = lidR::rasterize_terrain(ctg_roughness, 
                                res = 0.25, 
                                algorithm = tin())
#define file path and name for output dtm
out_file = file.path(roughness_dir_dtm, "dtm_mosaic_0_25m.tif")
#write output dtm to file
writeRaster(dtm_r, out_file, overwrite = T)
#call back in dtm
dtm_r = rast(out_file)
#turn off parallel processing
plan(sequential)

Create a ‘smoothed’ DTM by creating a circular focal matrix, then calculating focal values using the new fine-resolution DTM, the circular focal matrix as weights, and specifying the smoothing function as mean.

#generate focal matrix
circular_focal_matrix = focalMat(dtm_r, 2.5, "circle")
#create smooth dtm
smooth_dtm = terra::focal(dtm_r, circular_focal_matrix, fun=mean)

Perform terrain roughness calculations using new DTMs

#get absolute value of smooth/rough dtm difference as roughness measure
terr_rgh = abs(smooth_dtm - dtm_r)
#generate 10m res roughness raster and save to file
terr_rgh = aggregate(terr_rgh,
                     fact = 40,
                     fun = mean,
                     filename = paste0(roughness_dir, "roughness.tif"),
                     overwrite = TRUE)

Delete roughness lidar and DTM directories to save space in database

dir_delete(roughness_dir_laz)
dir_delete(roughness_dir_dtm)

Plot roughness raster!

plot(terrain_roughness)

Vegetation Density

Vegetation density is calculated using height normalized lidar point clouds, and finding a normalized relative density (NRD) value for all cells in the study area raster.

NRD is defined as follows: \[ NRD_{ij} = \frac{\sum_{i}^{j} n}{\sum_{0}^{j} n} \]

First we perform many of the same operations we did in the terrain roughness generation when setting up the calculation environment (creating directories, performing efficiency steps, etc.).

#start parallel processing
plan(multisession, workers = 10L)

#define vegetation density lidar directory
veg_dir_laz = paste0(veg_dir, "/veg_lidar")
#create new directory to copy lidar files into
dir.create(veg_dir_laz, showWarnings = T, recursive = T)

#define normalized height directory
veg_dir_norm = paste0(veg_dir, "/norm_ht")
#create new directory for normalized height files
dir.create(veg_dir_norm, showWarnings = T, recursive = T)

#copy lidar into new vegetation density lidar directory
file.copy(from = lidar_files,
          to = paste0(veg_dir_laz, lidar_files),
          overwrite = T,
          recursive = T)
#create new list of vegetation density lidar tile files
#specify only the .laz files and the full file name
laz_files = list.files(veg_dir_laz, 
                       pattern = "*.laz", 
                       full.names = T)
#create catalog of lidar to perform height normalizing operations
ctg_n = readLAScatalog(laz_files)

#direct operations to perform and save in new normalized height directory
opt_output_files(ctg_n) = paste0(veg_dir_norm, "/tile_{XLEFT}_{YBOTTOM}")
#reduce chunk size to reduce processing speed
opt_chunk_size(ctg_n) = 500
#reduce chunk buffer to reduce processing speed
opt_chunk_buffer(ctg_n) = 10

Now we perform the height normalizing operation as defined by the lidR package

# height normalizing command using new DTM
ctg_n = normalize_height(ctg_n, dtm)

Calculate vegetatIon density by calculating NRD

# set up NRD function
nrdfun = function(z){
  y = z[z >= 0.15 & z <= 2.75]
  x = z[z <= 2.75]
  nrd = length(y)/length(x)
  return(nrd)
}
# use pixel metrics to calculate NRD for every pixel in the height normalized point cloud
veg = pixel_metrics(ctg_n, func = ~nrdfun(Z), 10)

#turn off parallel processing
plan(sequential)

Plot results!

plot(vegetation_density)

---

Visibility

Visibility is another important metric used in this analysis. We generate visibility according to Mistick et al. (2023), using the VisiMod package in R.

The all-encompassing VisiMod() function is a single function which reads the dtm and dsm, defines the number of test points and the distance of radius for the visibility index (VI) calculation, determines if visibility is in a single direction or omnidirectional, and produces a visibility raster.

visibility = VisiMod(dtm = dtm,
                     dsm = dsm,
                     num_pts = 1000,
                     dist = 500,
                     vi_type = "omnidirectional",
                     save_dir = paste0(vis_dir, "/vis_files"),
                     cores = 4)

plot(vis)

Okay, I’m tired now… will return later to make prettier plots.