Automated Tool for Vertical Fuel Continuity Analysis using Airborne Laser Scanning Data.

LadderFuelsR: An R Package for vertical fuel continuity analysis using Airborne Laser Scanning data

Authors: Olga Viedma, Carlos Silva, JM Moreno and A.T. Hudak

Automated tool for vertical fuel continuity analysis using Airborne Laser Scanning data that can be applied on multiple tree species and for large-scale studies.The workflow consisted of 1) calculating the Leaf Area Density (LAD) profiles of each segmented tree; 2) identifying gaps and fuel layers; 3) estimating the distance between fuel layers; and 4) retrieving the fuel layers base height (FBH) and depth. Additionally, other functions recalculate previous metrics after considering distances greater than certain threshold and calculate the CBH based on three criteria: maximum LAD, and the largest- and the last-distance. Moreover, the package calculates: i) the percentage of LAD comprised in each fuel layer and remove fuel layers below a specified threshold (default 10 % LAD) recalculating the distances among the reminder ones. On the other hand, when the LAD profiles showed only one fuel layer with CBH at the minimum base height, it identifies the CBH performing a segmented linear regression (breaking point) on the cumulative sum of LAD as a function of height. Finally, a collection of plotting functions is developed to represent all previous metrics.

Getting Started

Installation

#The CRAN version:
install.packages("LadderFuelsR")

# The development version:
#install.packages("remotes")
library(remotes)
install_github("https://github.com/olgaviedma/LadderFuelsR", dependencies = TRUE)

# loading LadderFuelsR package
library(LadderFuelsR)

Required libraries

if (!require("pacman")) install.packages("pacman")
pacman::p_load(plyr,
               dplyr,
               tidyr,
               stringr,
               stringi,
               purrr,
               rlang,
               tidyverse,
               sf,
               terra,
               raster,
               data.table,
               rgdal,
               lidR,
               leafR,
               segmented,
               lidRplugins,
               ggplot2,
               gt,
               gridExtra,
               patchwork,
               SSBtools,
               tibble,
               rgl,
               rglwidget,
               LadderFuelsR,
               magrittr,
               gdata)

1. Computing Canopy height model (CHM) using lidR package

 LIDAR_dir <- file.path(system.file("extdata", package = "LadderFuelsR"))
 lidar_file<- lidR::readLAS(file.path(LIDAR_dir, "Eglin_zone1_clipped_000000.las"), filter = "-drop_z_below 0")

  chm_pitfree<- grid_canopy(lidar_file, res=0.5,pitfree( c(0,2,5,10,15,20,25,30,35,40), c(0,1.5), subcircle=0.15))
  chm_pitfree[chm_pitfree > 40] <- NA
  chm_pitfree[chm_pitfree < 0] <- 0
  chm_pitfree1 <- projectRaster(chm_pitfree, crs=26916)
  
col <- height.colors(25)
plot(chm_pitfree1,col=col)

2.Detecting individual tree top from the lidar-derived CHM

  # parameters
  ws= 2.5
  hmin = 2
  res=0.5
  ttops_multichm = find_trees(lidar_file, multichm(res = res, dist_2d = 2,ws= ws, layer_thickness = 0.3,dist_3d = 1, hmin = hmin))
  proj4string(ttops_multichm) <- CRS('+init=EPSG:26916')

# Create an rgl point cloud
x<-add_treetops3d(plot(lidar_file, bg = "white", size = 4), ttops_multichm)
# Customize the plot orientation
rgl.viewpoint(theta = 0, phi = 0, fov = 60, zoom = 0.75)
# Convert the rgl scene to an HTML widget
rglwidget(elementId = "x", width = 800, height = 600)

3. Individual tree crown deliniation (Silva et al. 2016)

  algo_silva1 <-silva2016(chm_pitfree1, ttops_multichm, max_cr_factor = 0.6, exclusion = 0.3, ID = "treeID")
  crowns_silva_las1 <-segment_trees(lidar_file, algo_silva1, attribute = "treeID", uniqueness = "incremental")
  crowns_silva_las2<-filter_poi(crowns_silva_las1, !is.na(treeID))
  
my_palette <- colorRampPalette(col)
x1<-plot(crowns_silva_las2, color = "treeID", pal = my_palette, bg = "white")

# Customize the plot orientation
rgl.viewpoint(theta = 0, phi = 0, fov = 10, zoom = 0.75)

# Convert the rgl scene to an HTML widget
rglwidget(elementId = "x1", width = 800, height = 600)

4. Defining function for computing crown-level metrics

custom_crown_metrics <- function(z, i) { # user-defined function
  metrics <- list(
     dz = 1,
     th = 1,
     z_max = max(z),# max height
     z_min = min(z),# min height
     z_mean = mean(z),# mean height
     z_sd = sd(z), # vertical variability of points
     z_q1=quantile(z, probs = 0.01),
     z_q5=quantile(z, probs = 0.05),
    z_q25=quantile(z, probs = 0.25),
    z_q50=quantile(z, probs = 0.50),
    z_q75=quantile(z, probs = 0.75),
    z_q95=quantile(z, probs = 0.95),
     crr=(mean(z)-min(z))/(max(z)-min(z))
   )
   return(metrics) # output
}
ccm = ~custom_crown_metrics(z = Z, i = Intensity)

5.Computing crown level standard metrics within all trees detected

 crowns_silva_filter<-filter_poi(crowns_silva_las2, Z >= 1)
  
  metrics1 = crown_metrics(crowns_silva_filter,func = .stdtreemetrics, geom = "convex")
  crown_diam<-data.frame(sqrt(metrics1$convhull_area/ pi) * 2)
  names(crown_diam)<-"crown_diam"
  metrics2 = crown_metrics(crowns_silva_filter,func = ccm, geom = "convex") #concave
  metrics_all <- dplyr::bind_cols(list(metrics1,crown_diam,metrics2))
  metrics_all1 <- metrics_all[,c(1:4,6,10:21)]
  names(metrics_all1)<-c("treeID", "Z", "npoints", "convhull_area", "crown_diam", "z_max", "z_min", "z_mean","z_sd", "z_q1","z_q5", "z_q25","z_q50","z_q75", "z_q95", "crr", "geometry" )
  
  tree_crowns <- st_as_sf(metrics_all1)
  
ttops1<-st_as_sf(ttops_multichm)
crowns1<-st_as_sf(tree_crowns)
ttops_within_crowns <- st_intersection(ttops1, crowns1)

# Set the size of the plotting device
par(mfrow = c(1, 1), mar = c(1, 1, 1, 1), pin = c(5, 4))
plot(st_geometry(crowns1), pch = 16, col = "green")
plot(ttops_within_crowns, add = TRUE, pch= 16, col = "darkblue", main = "Tree tops over the crowns")

6.Crop las files with crown polygons

  trees_ID <- tree_crowns %>% dplyr::select(treeID)
  n <- nrow(trees_ID)
  
  crown_cort <- vector("list", length=n)
  
  for (i in 1:n) {
    kk <- trees_ID[i,]
    crown_cort[[i]] = clip_roi(crowns_silva_las2, kk)
  }

my_palette <- colorRampPalette(col)
x2<-plot(crown_cort[[1]], color = "Z", pal = my_palette, bg = "black", size = 2.5)

# Customize the plot orientation
rgl.viewpoint(theta = 0, phi = 0, fov = 60, zoom = 0.75)

# Convert the rgl scene to an HTML widget
rglwidget(elementId = "x2", width = 400, height = 600)

7.LAI-LAD metrics by Trees

LIDAR_dir <- file.path(system.file("extdata", package = "LadderFuelsR"))
las_list1 <- list.files(LIDAR_dir, pattern = "*_CROWN.las", full.names = TRUE, ignore.case = TRUE)

# create a vector to hold the file names of .las files with more than 10 points
files_with_more_than_10_points <- c()

# loop through each file
for (file in las_list1) {
  las_data <- lidR::readLAS(file)
  las_data1<-filter_poi(las_data, Z >= 1)
  
  # skip to next file if there was a problem reading
  if (is.null(las_data1)) next
  
  # check if it contains more than three points
  if (las_data1@header$`Number of point records` > 10) {
  files_with_more_than_10_points <- c(files_with_more_than_10_points, file)
  }
}

# Creates a data frame of the 3D voxels information (xyz) with Leaf Area Density values
short_name1<-NULL
profile_list<-NULL
lidar_lai_list<-NULL
understory_lai_list<-NULL
LAHV_metric_list<-NULL

for (j in seq_along(files_with_more_than_10_points)){
  
  short_name<-stri_sub(files_with_more_than_10_points[j], 1, -5)
  short_name1<-gsub(".*/","",short_name)
  
  normlas_file<-files_with_more_than_10_points[[j]]

  VOXELS_LAD = lad.voxels(normlas_file, grain.size = 2)
  
  lad_profile = lad.profile(VOXELS_LAD, relative = F)
  lai_tot = lai(lad_profile)
  understory_lai <- lai(lad_profile, min = 0.3, max = 2.5)
  LAHV_metric<- LAHV(lad_profile, LAI.weighting = FALSE, height.weighting = FALSE)
  
  lad_profile1 = data.frame(lad_profile, treeID = short_name1)
  lai_tot1 = data.frame(lai_tot, treeID = short_name1)
  understory_lai1 = data.frame(understory_lai, treeID = short_name1)
  LAHV_metric1 = data.frame(LAHV_metric, treeID = short_name1)
  
  profile_list<-rbind(profile_list, lad_profile1)
  lidar_lai_list<-rbind(lidar_lai_list,lai_tot1)
  understory_lai_list <-rbind(understory_lai_list,understory_lai1)
  LAHV_metric_list<-rbind(LAHV_metric_list,LAHV_metric1)
}

head(profile_list,10)

8.Depurating Tree LAD profiles

cols <- c('treeID')
profile_list[cols] <- lapply(profile_list[cols], function (x) as.factor(x))
profile_list$lad<-round(profile_list$lad,digits = 4)

cases <- data.frame(table(profile_list$treeID))
cases1 <-cases[cases$Freq > 5, ]
names(cases1)<-c("treeID", "Freq")

profile_list1 <- profile_list[profile_list$treeID %in% cases1$treeID, ]
profile_list2 <- data.frame(profile_list1 %>% replace(is.na(.), 0.01))

9.Gaps and Fuel Layers Base Height (FBH)

# LAD profiles derived from normalized ALS data after applying [lad.profile()] function from leafR package
profile_list2$treeID <- factor(profile_list2$treeID)

trees_name1 <- as.character(profile_list2$treeID)
trees_name2 <- factor(unique(trees_name1))

gaps_fbhs_list<-list()
for (i in levels(trees_name2)) {
  tree2 <- profile_list2 |> dplyr::filter(treeID == i)
  gaps_fbhs <- get_gaps_fbhs(tree2, step=1,
                          min_height=1.5,
                          perc_gap= 25,perc_base= 25,
                          verbose=TRUE)
  gaps_fbhs_list[[i]] <- gaps_fbhs
}

gaps_fbhs_list1 <- dplyr::bind_rows(gaps_fbhs_list)
gaps_fbhs_list1$treeID <- factor(gaps_fbhs_list1$treeID)

# Remove the row with all NA values from the original data frame
# First remove "treeID" and "treeID1" columns
gaps_fbhs_list1_no_treeID <- gaps_fbhs_list1[, -which(names(gaps_fbhs_list1) == c("treeID","treeID1"))]
# Check if any row has all NA values
rows_with_all_NA_or_zero <- apply(gaps_fbhs_list1_no_treeID, 1, function(row) all(is.na(row) | row == 0))
# Get the row index with all NA values
row_index <- which(rows_with_all_NA_or_zero)

# Remove the row with all NA values from the original data frame
if (length(row_index) > 0) {
  gaps_fbhs_metrics <- gaps_fbhs_list1[-row_index, ]
} else {
  gaps_fbhs_metrics <- gaps_fbhs_list1
}
rownames(gaps_fbhs_metrics) <- NULL
head(gaps_fbhs_metrics)

10.LAD percentile of each height bin

# LAD profiles derived from normalized ALS data after applying [lad.profile()] function from leafR package
profile_list2$treeID <- factor(profile_list2$treeID)

trees_name1 <- as.character(profile_list2$treeID)
trees_name2 <- factor(unique(trees_name1))

gaps_perc_list <- list()  # Initialize outside the loop

    for (i in levels(trees_name2)) {
      tree1 <- profile_list2 |> dplyr::filter(treeID == i)
      percentiles <- calculate_gaps_perc(tree1, min_height=1.5)
      gaps_perc_list[[i]] <- percentiles
    }

    gaps_perc <- dplyr::bind_rows(gaps_perc_list)
 
head(gaps_perc)

11.Distance between Fuel Layers

# Tree metrics derived from get_gaps_fbhs() function
numeric_vars <- setdiff(names(gaps_fbhs_metrics), c("treeID", "treeID1"))
gaps_fbhs_metrics[numeric_vars] <- lapply(gaps_fbhs_metrics[numeric_vars], function(x) as.numeric(ifelse(x == "NA", NA, x)))
gaps_fbhs_metrics$treeID <- factor(gaps_fbhs_metrics$treeID)

# Tree metrics derived from calculate_gaps_perc() function
gaps_perc$treeID <- factor(gaps_perc$treeID)

trees_name1 <- as.character(gaps_fbhs_metrics$treeID)
trees_name2 <- factor(unique(trees_name1))

metrics_distance_list <- list()

for (i in levels(trees_name2)) {

  # Filter data for each tree
  tree1 <- gaps_fbhs_metrics |> dplyr::filter(treeID == i)
  tree2 <- gaps_perc |> dplyr::filter(treeID == i)

  # Get distance metrics for each tree
  metrics_distance <- get_distance(tree1,tree2,step=1, min_height=1.5)
  metrics_distance_list[[i]] <- metrics_distance
}

# Combine the individual data frames
distance_metrics <- dplyr::bind_rows(metrics_distance_list)
distance_metrics <- distance_metrics[, order(names(distance_metrics))]
rownames(distance_metrics) <- NULL
head(distance_metrics)

12.Fuel Layers Depth

library(dplyr)
library(magrittr)

# LAD profiles derived from normalized ALS data after applying [lad.profile()] function
profile_list2$treeID <- factor(profile_list2$treeID)
# Tree metrics derived from get_distance() function
distance_metrics$treeID <- factor(distance_metrics$treeID)

metrics_depth_list <- list()

for (i in levels(profile_list2$treeID)){

  tree1 <- profile_list2 |> dplyr::filter(treeID == i)
  tree2 <- distance_metrics |> dplyr::filter(treeID == i)

  # Get depths for each tree
  metrics_depth <- get_depths(tree1, tree2,step= 1,min_height= 1.5)
  metrics_depth_list[[i]] <- metrics_depth
}

# Combine the individual data frames
depth_metrics <- dplyr::bind_rows(metrics_depth_list)

depth_metrics <- depth_metrics[, order(names(depth_metrics))]
rownames(depth_metrics) <- NULL
head(depth_metrics)

13.Plot Gaps and Fuel Layers Base Height (FBH)

library(LadderFuelsR)
library(ggplot2)
library(lattice)

# LAD profiles derived from normalized ALS data after applying [lad.profile()] function
profile_list2$treeID <- factor(profile_list2$treeID)
# Tree metrics derived from get_gaps_fbhs() function
gaps_fbhs_metrics$treeID <- factor(gaps_fbhs_metrics$treeID)

# Generate plots for gaps and fbhs
plots_gaps_fbhs <- get_plots_gap_fbh(profile_list2, gaps_fbhs_metrics,min_height=1.5)

par(mfrow = c(2, 2))
# Plot in RED are the GAPS and in GREEN the FBHs
plot(plots_gaps_fbhs[[1]])
plot(plots_gaps_fbhs[[2]])
plot(plots_gaps_fbhs[[3]])

{alt="Plot 1"}	{alt="Plot 2"}
{alt="Plot 3"}

14.Fuel Layers Base Height (FBH) after after considering distances greater than any number of height bin steps

library(SSBtools)
library(dplyr)
library(magrittr)

# Tree metrics derived from get_depths() function
depth_metrics$treeID <- factor(depth_metrics$treeID)

trees_name1 <- as.character(depth_metrics$treeID)
trees_name2 <- factor(unique(trees_name1))

fbh_corr_list <- list()

for (i in levels(trees_name2)){

  # Filter data for each tree
  tree3 <- depth_metrics |> dplyr::filter(treeID ==i)

  # Get real fbh for each tree
  fbh_corr <- get_real_fbh(tree3,step= 1, number_steps = 1, min_height=1.5)

  # Store fbh values in a list
  fbh_corr_list[[i]] <- fbh_corr
}

# Combine fbh values for all trees
fbh_metrics_corr <- dplyr::bind_rows(fbh_corr_list)
fbh_metrics_corr$treeID <- factor(fbh_metrics_corr$treeID)

# Reorder columns
# Get original column names
original_column_names <- colnames(fbh_metrics_corr)

# Specify prefixes
prefixes <- c("treeID", "Hdist", "Hcbh", "Hdepth", "dist", "depth", "max_height")

# Initialize vector to store new order
new_order <- c()

# Loop over prefixes
for (prefix in prefixes) {
  # Find column names matching the current prefix
  matching_columns <- grep(paste0("^", prefix), original_column_names, value = TRUE)
  # Append to the new order
  new_order <- c(new_order, matching_columns)
}

# Reorder values
fbh_metrics_corr <- fbh_metrics_corr[, new_order]
rownames(fbh_metrics_corr) <- NULL
head(fbh_metrics_corr)

15.Fuel Layers Depth after removing distances greater than the actual height bin step

library(dplyr)
library(magrittr)
library(tidyr)

# Tree metrics derived from get_real_fbh() function
fbh_metrics_corr$treeID <- factor(fbh_metrics_corr$treeID)

trees_name1 <- as.character(fbh_metrics_corr$treeID)
trees_name2 <- factor(unique(trees_name1))

depth_metrics_corr_list <- lapply(levels(trees_name2), function(i) {
  # Filter data for each tree
  tree2 <- fbh_metrics_corr |> dplyr::filter(treeID == i)
  # Get real depths for each tree
  get_real_depths(tree2,step=1, min_height=1.5)
})

# Combine depth values for all trees
depth_metrics_corr <- dplyr::bind_rows(depth_metrics_corr_list)
rownames(depth_metrics_corr) <- NULL
head(depth_metrics_corr)

16.Distance between Fuel Layers greater than any number of height bin steps

library(dplyr)
library(magrittr)
library(stringr)

# Tree metrics derived from get_real_depths() function
depth_metrics_corr$treeID <- factor(depth_metrics_corr$treeID)

trees_name1 <- as.character(depth_metrics_corr$treeID)
trees_name2 <- factor(unique(trees_name1))

distance_metrics_corr_list <- lapply(levels(trees_name2), function(i) {
  # Filter data for each tree
  tree2 <- depth_metrics_corr |> dplyr::filter(treeID == "1_CROWN")
  # Get effective gap for each tree
  get_effective_gap(tree2,number_steps = 1, min_height= 1.5)
})

# Combine the individual data frames
distances_metrics_corr <- dplyr::bind_rows(distance_metrics_corr_list)

# =======================================================================#
# REORDER COLUMNS:
# =======================================================================#
# Get original column names
original_column_names <- colnames(distances_metrics_corr)

# Specify prefixes
prefixes <- c("treeID", "Hcbh", "dptf", "Hdptf", "effdist", "dist", "Hdist", "max_Hcbh", "max_dptf", "max_Hdptf", "last_Hcbh", "last_dptf", "last_Hdptf", "max_height")

# Initialize vector to store new order
new_order <- c()

# Loop over prefixes
for (prefix in prefixes) {
  # Find column names matching the current prefix
  matching_columns <- grep(paste0("^", prefix), original_column_names, value = TRUE)

  # Extract numeric suffixes and order the columns based on these suffixes
  numeric_suffixes <- as.numeric(gsub(paste0("^", prefix), "", matching_columns))
  matching_columns <- matching_columns[order(numeric_suffixes)]

  # Append to new order
  new_order <- c(new_order, matching_columns)
}

# Reorder values
distances_metrics_corr1 <- distances_metrics_corr[, new_order]
# Unlist the data frame
distances_metrics_corr2 <- as.data.frame(lapply(distances_metrics_corr1, function(x) unlist(x)))
rownames(distances_metrics_corr2) <- NULL
head(distances_metrics_corr2)

17.Fuels LAD percentage (greater than a threshold)

library(dplyr)
library(magrittr)

# LAD profiles derived from normalized ALS data after applying [lad.profile()] function
profile_list2$treeID <- factor(profile_list2$treeID)

# Tree metrics derived from get_effective_gap() function
distances_metrics_corr2$treeID <- factor(distances_metrics_corr2$treeID)

trees_name1 <- as.character(distances_metrics_corr2$treeID)
trees_name2 <- factor(unique(trees_name1))

LAD_metrics1 <- list()
LAD_metrics2 <- list()

for (i in levels(trees_name2)) {
  # Filter data for each tree
  tree1 <- profile_list2 |> dplyr::filter(treeID == i)
  tree2 <- distances_metrics_corr2 |> dplyr::filter(treeID ==i)

  # Get LAD metrics for each tree
  LAD_metrics <- get_layers_lad(tree1, tree2, threshold=10,
                           step = 1,min_height= 1.5)
  LAD_metrics1[[i]] <- LAD_metrics$df1
  LAD_metrics2[[i]] <- LAD_metrics$df2
}

LAD_metrics_all1 <- dplyr::bind_rows(LAD_metrics1)
LAD_metrics_all2 <- dplyr::bind_rows(LAD_metrics2)

# List of data frames
LAD_metrics_list <- list(LAD_metrics_all1, LAD_metrics_all2)

# Initialize an empty list to store reordered data frames
fuels_LAD_metrics <- list()

# Specify prefixes (adjust accordingly)
prefixes <- c("treeID", "Hdist", "Hcbh", "effdist", "dptf", "Hdptf", "max", "last", "nlayers")

# Loop over each data frame
for (i in seq_along(LAD_metrics_list)) {

  LAD_metrics_all <- LAD_metrics_list[[i]]

  # Get original column names
  original_column_names <- colnames(LAD_metrics_all)

  # Initialize vector to store new order
  new_order <- c()

  # Loop over prefixes
  for (prefix in prefixes) {
    # Find column names matching the current prefix
    matching_columns <- grep(paste0("^", prefix), original_column_names, value = TRUE)

    # Extract numeric suffixes and order the columns based on these suffixes
    numeric_suffixes <- as.numeric(gsub(paste0("^", prefix), "", matching_columns))
    
       # Order the columns based on numeric suffixes
    matching_columns <- matching_columns[order(numeric_suffixes)]

    # Append to new order
    new_order <- c(new_order, matching_columns)
  }
  # Reorder columns
  LAD_metrics_all <- LAD_metrics_all[, new_order]
  # Store the reordered data frame in the list
  fuels_LAD_metrics[[i]] <- LAD_metrics_all
 }
rownames(fuels_LAD_metrics[[1]]) <- NULL
rownames(fuels_LAD_metrics[[2]]) <- NULL

head(fuels_LAD_metrics[[2]])

18.Plot Effective Fuel Layers with LAD percentage greater than a threshold

library(ggplot2)

# LAD profiles derived from normalized ALS data after applying [lad.profile()] function
profile_list2$treeID <- factor(profile_list2$treeID)
# Tree metrics derived from get_layers_lad() function
LAD_gt10p <- fuels_LAD_metrics[[2]]

trees_name1 <- as.character(LAD_gt10p$treeID)
trees_name2 <- factor(unique(trees_name1))

# Generate plots for fuels LAD metrics
plots_trees_LAD <- get_plots_effective(profile_list2, LAD_gt10p, min_height=1.5)

par(mfrow = c(2, 2))
plot(plots_trees_LAD[[1]])
plot(plots_trees_LAD[[2]])
plot(plots_trees_LAD[[3]])

{alt="Plot 1"}	{alt="Plot 2"}
{alt="Plot 3"}

19.CBH based on different criteria: maximum LAD, maximum and last distance

library(dplyr)
library(magrittr)

# Tree metrics derived from get_layers_lad() function
LAD_gt10p <- fuels_LAD_metrics[[2]]

trees_name1 <- as.character(LAD_gt10p$treeID)
trees_name2 <- factor(unique(trees_name1))

cbh_metrics_list <- list()

for (j in levels(trees_name2)){

  # Filter data for each tree
  tree1 <- LAD_gt10p |> dplyr::filter(treeID == j)
  cbh_metrics <- get_cbh_metrics(tree1,min_height= 1.5)
  cbh_metrics_list[[j]] <- cbh_metrics
}

# Combine depth values for all trees
cbh_metrics_all <- dplyr::bind_rows(cbh_metrics_list)

# Get original column names
  original_column_names <- colnames(cbh_metrics_all)

  # Specify prefixes
desired_order <- c("treeID", "Hcbh", "dptf","effdist","dist", "Hdist", "Hdptf","maxlad_","max_","last_","nlayers")

  # Identify unique prefixes
  prefixes <- unique(sub("^([a-zA-Z]+).*", "\\1", original_column_names))
  # Initialize vector to store new order
  new_order <- c()

  # Loop over desired order of prefixes
  for (prefix in desired_order) {
    # Find column names matching the current prefix
    matching_columns <- grep(paste0("^", prefix), original_column_names, value = TRUE)
    # Append to the new order
    new_order <- c(new_order, matching_columns)
  }
  
  # Reorder columns
  cbh_metrics_all <- cbh_metrics_all[, new_order]

20.Plots CBH based on different criteria: maximum LAD, maximum and last distance

library(ggplot2)

# LAD profiles derived from normalized ALS data after applying [lad.profile()] function
profile_list2$treeID <- factor(profile_list2$treeID)
# Tree metrics derived from get_cbh_metrics() function
cbh_metrics_all$treeID <- factor(cbh_metrics_all$treeID)

trees_name1 <- as.character(cbh_metrics_all$treeID)
trees_name2 <- factor(unique(trees_name1))

# Generate plots for fuels LAD metrics
plots_cbh_maxlad <- get_plots_cbh_LAD(profile_list2, cbh_metrics_all,min_height=1.5)
plots_cbh_maxdist <- get_plots_cbh_maxdist(profile_list2, cbh_metrics_all,min_height=1.5)
plots_cbh_lastdist <- get_plots_cbh_lastdist(profile_list2, cbh_metrics_all,min_height=1.5)

par(mfrow = c(2, 2))
plot(plots_cbh_maxlad[[1]])
plot(plots_cbh_maxdist[[2]])
plot(plots_cbh_lastdist[[3]])

{alt="Plot 1"}	{alt="Plot 2"}
{alt="Plot 3"}

21.CBH based on the Breaking Point method and LAD percentage

library(dplyr)
library(magrittr)

# LAD profiles derived from normalized ALS data after applying [lad.profile()] function
profile_list2$treeID <- factor(profile_list2$treeID)

# Tree metrics derived from get_cbh_metrics() function
cbh_metrics_all$treeID <- factor(cbh_metrics_all$treeID)

trees_name1 <- as.character(cbh_metrics_all$treeID)
trees_name2 <- factor(unique(trees_name1))

cum_LAD_metrics_list <- list()

for (i in levels(trees_name2)) {
  # Filter data for each tree
  tree1 <- profile_list2 |> dplyr::filter(treeID == i)
  tree2 <- cbh_metrics_all |> dplyr::filter(treeID == "1_CROWN")

  # Get cumulative LAD metrics for each tree
  cum_LAD_metrics_all <- get_cum_break(tree1, tree2, threshold=75, min_height= 1.5, verbose=TRUE)
  cum_LAD_metrics_list[[i]] <- cum_LAD_metrics_all
}

# Combine the individual data frames
cum_LAD_metrics <- dplyr::bind_rows(cum_LAD_metrics_list)

# =======================================================================#
# REORDER COLUMNS
# =======================================================================#

# Get original column names
original_column_names <- colnames(cum_LAD_metrics)

# Specify prefixes (adjust accordingly)
prefixes <- c("treeID", "Hcbh", "below", "above", "bp", "max", "cumlad")

# Initialize vector to store new order
new_order <- c()

# Loop over prefixes
for (prefix in prefixes) {
  # Find column names matching the current prefix
  matching_columns <- grep(paste0("^", prefix), original_column_names, value = TRUE)

  # Extract numeric suffixes and order the columns based on these suffixes
  numeric_suffixes <- as.numeric(gsub(paste0("^", prefix), "", matching_columns))
  matching_columns <- matching_columns[order(numeric_suffixes)]

  # Append to new order
  new_order <- c(new_order, matching_columns)
}

# Reorder columns
cum_LAD_metrics <- cum_LAD_metrics[, new_order]

## when % LAD is < 75 % below or above the BP (Breaking Point), Hcbh1 is derived from CBH maximum LAD criterium
rownames(cum_LAD_metrics) <- NULL
head(cum_LAD_metrics)

22.Plot CBH based on the Breaking Point method and LAD percentage

library(ggplot2)

# LAD profiles derived from normalized ALS data after applying [lad.profile()] function
profile_list2$treeID <- factor(profile_list2$treeID)

# Tree metrics derived from get_cum_break() function
cum_LAD_metrics$treeID <- factor(cum_LAD_metrics$treeID)

trees_name1 <- as.character(cum_LAD_metrics_all$treeID)
trees_name2 <- factor(unique(trees_name1))

plot_list <- list()

  for (j in levels(trees_name2)) {
    
    tree1 <- profile_list2 |> dplyr::filter(treeID ==j)
    tree2 <- cum_LAD_metrics |> dplyr::filter(treeID == j)
    
    plots_cbh_bp <- get_plots_cbh_bp(tree1, tree2,min_height = 1.5)
    plot_list[[j]] <- plots_cbh_bp
    
  }

par(mfrow = c(2, 2))
plot(plot_list[[1]])
plot(plot_list[[2]])
plot(plot_list[[3]])

{alt="Plot 1"}	{alt="Plot 2"}
{alt="Plot 3"}

23. Joinning Fuel ladder properties with Crown polygons

# Tree metrics derived from get_layers_lad() function
cbh_metrics_all$treeID1 <- factor(cbh_metrics_all$treeID1)

# crown polygons (output from step 4)
tree_crowns$treeID1 <- factor(tree_crowns$treeID)

crowns_properties<-merge (tree_crowns,cbh_metrics_all, by="treeID1")
crowns_properties$maxlad_Hcbh_factor <- cut(crowns_properties$maxlad_Hcbh, breaks = 5)

# Plotting with a discrete legend

palette <- colorRampPalette(c("orange", "dark green"))

ggplot() +
  geom_sf(data = crowns_properties, aes(fill = maxlad_Hcbh_factor)) +
  scale_fill_manual(values = palette(5)) +
  theme_minimal() +
  labs(title = "Tree Crowns", fill = "maxlad_Hcbh")

Acknowledgements

We gratefully acknowledge funding from project INFORICAM (PID2020-119402RB-I00), funded by the Spanish MCIN/AEI/ 10.13039/501100011033 and by the "European Union NextGenerationEU/PRTR". Carlos Silva was supported by the NASA's Carbon Monitoring System funding (CMS, grant 22-CMS22-0015).

Reporting Issues

Please report any issue regarding the LadderFuelsR package to Dr. Olga Viedma ([email protected]{.email})

Citing LadderFuelsR

Viedma,O.;Silva, C; Moreno, JM & Hudak, AT: LadderFuelsR: An R Package for vertical fuel continuity analysis using LiDAR data.version 0.0.1, accessed on November. 22 2023, available at: https://CRAN.R-project.org/package=LadderFuelsR

Disclaimer

LadderFuelsR package comes with no guarantee, expressed or implied, and the authors hold no responsibility for its use or reliability of its outputs.