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.This is an updated version.

Getting Started

Installation

``` r

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

```{r pressure, echo=FALSE} 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

```{r CHM pitfree 0.5 m, echo=TRUE, message=FALSE, warning=FALSE}

LIDARdir <- file.path(system.file("extdata", package = "LadderFuelsR")) lidarfile<- lidR::readLAS(file.path(LIDARdir, "Eglinzone1clipped000000.las"), filter = "-dropzbelow 0")

chmpitfree<- gridcanopy(lidarfile, res=0.5,pitfree( c(0,2,5,10,15,20,25,30,35,40), c(0,1.5), subcircle=0.15)) chmpitfree[chmpitfree > 40] <- NA chmpitfree[chmpitfree < 0] <- 0 chmpitfree1 <- projectRaster(chm_pitfree, crs=26916)

col <- height.colors(25) plot(chm_pitfree1,col=col)

```

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

```{r Tree tops detection, echo=TRUE, message=FALSE, warning=FALSE}

# parameters ws= 2.5 hmin = 2 res=0.5 ttopsmultichm = findtrees(lidarfile, multichm(res = res, dist2d = 2,ws= ws, layerthickness = 0.3,dist3d = 1, hmin = hmin)) proj4string(ttops_multichm) <- CRS('+init=EPSG:26916')

Create an rgl point cloud

x<-addtreetops3d(plot(lidarfile, 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)

```{r Crowns Silva, echo=TRUE, message=FALSE, warning=FALSE}

algosilva1 <-silva2016(chmpitfree1, ttopsmultichm, maxcrfactor = 0.6, exclusion = 0.3, ID = "treeID") crownssilvalas1 <-segmenttrees(lidarfile, algosilva1, attribute = "treeID", uniqueness = "incremental") crownssilvalas2<-filterpoi(crownssilva_las1, !is.na(treeID))

mypalette <- colorRampPalette(col) x1<-plot(crownssilvalas2, color = "treeID", pal = mypalette, 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

```{r tree metrics function, echo=TRUE}

customcrownmetrics <- function(z, i) { # user-defined function metrics <- list( dz = 1, th = 1, zmax = max(z),# max height zmin = min(z),# min height zmean = mean(z),# mean height zsd = sd(z), # vertical variability of points zq1=quantile(z, probs = 0.01), zq5=quantile(z, probs = 0.05), zq25=quantile(z, probs = 0.25), zq50=quantile(z, probs = 0.50), zq75=quantile(z, probs = 0.75), zq95=quantile(z, probs = 0.95), crr=(mean(z)-min(z))/(max(z)-min(z)) ) return(metrics) # output } ccm = ~customcrownmetrics(z = Z, i = Intensity)

```

5.Computing crown level standard metrics within all trees detected

```{r tree and crown standard and own metrics, echo=TRUE, message=FALSE, warning=FALSE} crownssilvafilter<-filterpoi(crownssilva_las2, Z >= 1)

metrics1 = crownmetrics(crownssilvafilter,func = .stdtreemetrics, geom = "convex") crowndiam<-data.frame(sqrt(metrics1$convhullarea/ pi) * 2) names(crowndiam)<-"crowndiam" metrics2 = crownmetrics(crownssilvafilter,func = ccm, geom = "convex") #concave metricsall <- dplyr::bindcols(list(metrics1,crowndiam,metrics2)) metricsall1 <- metricsall[,c(1:4,6,10:21)] names(metricsall1)<-c("treeID", "Z", "npoints", "convhullarea", "crowndiam", "zmax", "zmin", "zmean","zsd", "zq1","zq5", "zq25","zq50","zq75", "zq95", "crr", "geometry" )

treecrowns <- stassf(metricsall1)

ttops1<-stassf(ttopsmultichm) crowns1<-stassf(treecrowns) ttopswithincrowns <- 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(stgeometry(crowns1), pch = 16, col = "green") plot(ttopswithin_crowns, add = TRUE, pch= 16, col = "darkblue", main = "Tree tops over the crowns") ```

6.Crop las files with crown polygons

```{r cropLAS files with no overlapping crowns, echo=TRUE, message=FALSE, warning=FALSE}

treesID <- treecrowns %>% dplyr::select(treeID) n <- nrow(trees_ID)

crown_cort <- vector("list", length=n)

for (i in 1:n) { kk <- treesID[i,] crowncort[[i]] = cliproi(crownssilva_las2, kk) }

mypalette <- colorRampPalette(col) x2<-plot(crowncort[[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

```{r LAI and LAD tree metrics, echo=TRUE, message=FALSE, warning=FALSE}

LIDARdir <- file.path(system.file("extdata", package = "LadderFuelsR")) laslist1 <- list.files(LIDARdir, 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

fileswithmorethan10_points <- c()

loop through each file

for (file in laslist1) { lasdata <- lidR::readLAS(file) lasdata1<-filterpoi(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 (lasdata1@header$Number of point records > 10) { fileswithmorethan10points <- c(fileswithmorethan10_points, file) } }

Creates a data frame of the 3D voxels information (xyz) with Leaf Area Density values

shortname1<-NULL profilelist<-NULL lidarlailist<-NULL understorylailist<-NULL LAHVmetriclist<-NULL

for (j in seqalong(fileswithmorethan10points)){

shortname<-strisub(fileswithmorethan10points[j], 1, -5) shortname1<-gsub(".*/","",short_name)

normlasfile<-fileswithmorethan10points[[j]]

VOXELSLAD = lad.voxels(normlasfile, grain.size = 2)

ladprofile = lad.profile(VOXELSLAD, relative = F) laitot = lai(ladprofile) understorylai <- lai(ladprofile, min = 0.3, max = 2.5) LAHVmetric<- LAHV(ladprofile, LAI.weighting = FALSE, height.weighting = FALSE)

ladprofile1 = data.frame(ladprofile, treeID = shortname1) laitot1 = data.frame(laitot, treeID = shortname1) understorylai1 = data.frame(understorylai, treeID = shortname1) LAHVmetric1 = data.frame(LAHVmetric, treeID = shortname1)

profilelist<-rbind(profilelist, ladprofile1) lidarlailist<-rbind(lidarlailist,laitot1) understorylailist <-rbind(understorylailist,understorylai1) LAHVmetriclist<-rbind(LAHVmetriclist,LAHVmetric1) }

head(profile_list,10) ```

8.Depurating Tree LAD profiles

```{r depurating LAD databases, echo=TRUE, message=FALSE, warning=FALSE}

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

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

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

```

9.Gaps and Fuel Layers Base Height (FBH)

```{r Gaps and Fuel layers Base Height (fbh), echo=TRUE, message=FALSE, warning=FALSE}

LAD profiles derived from normalized ALS data after applying [lad.profile()] function from leafR package

profilelist2$treeID <- factor(profilelist2$treeID)

treesname1 <- as.character(profilelist2$treeID) treesname2 <- factor(unique(treesname1))

gapsfbhslist<-list() for (i in levels(treesname2)) { tree2 <- profilelist2 |> dplyr::filter(treeID == i) gapsfbhs <- getgapsfbhs(tree2, step=1, minheight=1.5, percgap= 25,percbase= 25, verbose=TRUE) gapsfbhslist[[i]] <- gaps_fbhs }

gapsfbhslist1 <- dplyr::bindrows(gapsfbhslist) gapsfbhslist1$treeID <- factor(gapsfbhs_list1$treeID)

Remove the row with all NA values from the original data frame

First remove "treeID" and "treeID1" columns

gapsfbhslist1notreeID <- gapsfbhslist1[, -which(names(gapsfbhslist1) == c("treeID","treeID1"))]

Check if any row has all NA values

rowswithallNAorzero <- apply(gapsfbhslist1no_treeID, 1, function(row) all(is.na(row) | row == 0))

Get the row index with all NA values

rowindex <- which(rowswithallNAorzero)

Remove the row with all NA values from the original data frame

if (length(rowindex) > 0) { gapsfbhsmetrics <- gapsfbhslist1[-rowindex, ] } else { gapsfbhsmetrics <- gapsfbhslist1 } rownames(gapsfbhsmetrics) <- NULL head(gapsfbhsmetrics) ```

10.LAD percentile of each height bin

```{r LAD percentile of each height bin, echo=TRUE, message=FALSE, warning=FALSE}

LAD profiles derived from normalized ALS data after applying [lad.profile()] function from leafR package

profilelist2$treeID <- factor(profilelist2$treeID)

treesname1 <- as.character(profilelist2$treeID) treesname2 <- factor(unique(treesname1))

gapsperclist <- 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

```{r Distances (and their heights) between fuel layers, echo=TRUE, message=FALSE, warning=FALSE}

Tree metrics derived from getgapsfbhs() function

numericvars <- setdiff(names(gapsfbhsmetrics), c("treeID", "treeID1")) gapsfbhsmetrics[numericvars] <- lapply(gapsfbhsmetrics[numericvars], function(x) as.numeric(ifelse(x == "NA", NA, x))) gapsfbhsmetrics$treeID <- factor(gapsfbhs_metrics$treeID)

Tree metrics derived from calculategapsperc() function

gapsperc$treeID <- factor(gapsperc$treeID)

treesname1 <- as.character(gapsfbhsmetrics$treeID) treesname2 <- factor(unique(trees_name1))

metricsdistancelist <- list()

for (i in levels(trees_name2)) {

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

# Get distance metrics for each tree metricsdistance <- getdistance(tree1,tree2,step=1, minheight=1.5) metricsdistancelist[[i]] <- metricsdistance }

Combine the individual data frames

distancemetrics <- dplyr::bindrows(metricsdistancelist) distancemetrics <- distancemetrics[, order(names(distancemetrics))] rownames(distancemetrics) <- NULL head(distance_metrics) ```

12.Fuel Layers Depth

```{r Distane between fuel layers, echo=TRUE, message=FALSE, warning=FALSE}

library(dplyr) library(magrittr)

LAD profiles derived from normalized ALS data after applying [lad.profile()] function

profilelist2$treeID <- factor(profilelist2$treeID)

Tree metrics derived from get_distance() function

distancemetrics$treeID <- factor(distancemetrics$treeID)

metricsdepthlist <- list()

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

tree1 <- profilelist2 |> dplyr::filter(treeID == i) tree2 <- distancemetrics |> dplyr::filter(treeID == i)

# Get depths for each tree metricsdepth <- getdepths(tree1, tree2,step= 1,minheight= 1.5) metricsdepthlist[[i]] <- metricsdepth }

Combine the individual data frames

depthmetrics <- dplyr::bindrows(metricsdepthlist)

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

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

```{r Plots Gaps and Fuel layers Base Height (fbh), echo=TRUE, message=FALSE, warning=FALSE}

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

LAD profiles derived from normalized ALS data after applying [lad.profile()] function

profilelist2$treeID <- factor(profilelist2$treeID)

Tree metrics derived from getgapsfbhs() function

gapsfbhsmetrics$treeID <- factor(gapsfbhsmetrics$treeID)

Generate plots for gaps and fbhs

plotsgapsfbhs <- getplotsgapfbh(profilelist2, gapsfbhsmetrics,min_height=1.5)

par(mfrow = c(2, 2))

Plot in RED are the GAPS and in GREEN the FBHs

plot(plotsgapsfbhs[[1]]) plot(plotsgapsfbhs[[2]]) plot(plotsgapsfbhs[[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

```{r Fuels base height recalculation after after considering distances greater than any number of height bin steps, echo=TRUE, message=FALSE, warning=FALSE}

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

Tree metrics derived from get_depths() function

depthmetrics$treeID <- factor(depthmetrics$treeID)

treesname1 <- as.character(depthmetrics$treeID) treesname2 <- factor(unique(treesname1))

fbhcorrlist <- 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 fbhcorr <- getrealfbh(tree3,step= 1, numbersteps = 1, min_height=1.5)

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

Combine fbh values for all trees

fbhmetricscorr <- dplyr::bindrows(fbhcorrlist) fbhmetricscorr$treeID <- factor(fbhmetrics_corr$treeID)

Reorder columns

Get original column names

originalcolumnnames <- colnames(fbhmetricscorr)

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 matchingcolumns <- grep(paste0("^", prefix), originalcolumnnames, value = TRUE) # Append to the new order neworder <- c(neworder, matchingcolumns) }

Reorder values

fbhmetricscorr <- fbhmetricscorr[, neworder] rownames(fbhmetricscorr) <- NULL head(fbhmetrics_corr) ```

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

```{r Fuel layers depth after considering distances greater than the actual height bin step, echo=TRUE, message=FALSE, warning=FALSE}

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

Tree metrics derived from getrealfbh() function

fbhmetricscorr$treeID <- factor(fbhmetricscorr$treeID)

treesname1 <- as.character(fbhmetricscorr$treeID) treesname2 <- factor(unique(trees_name1))

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

Combine depth values for all trees

depthmetricscorr <- dplyr::bindrows(depthmetricscorrlist) rownames(depthmetricscorr) <- NULL head(depthmetricscorr) ```

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

```{r Fuel layers distances after considering distances greater than any number of height bin steps, echo=TRUE, message=FALSE, warning=FALSE}

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

Tree metrics derived from getrealdepths() function

depthmetricscorr$treeID <- factor(depthmetricscorr$treeID)

treesname1 <- as.character(depthmetricscorr$treeID) treesname2 <- factor(unique(trees_name1))

distancemetricscorrlist <- lapply(levels(treesname2), function(i) { # Filter data for each tree tree2 <- depthmetricscorr |> dplyr::filter(treeID == "1CROWN") # Get effective gap for each tree geteffectivegap(tree2,numbersteps = 1, min_height= 1.5) })

Combine the individual data frames

distancesmetricscorr <- dplyr::bindrows(distancemetricscorrlist)

=======================================================================

REORDER COLUMNS:

=======================================================================

Get original column names

originalcolumnnames <- colnames(distancesmetricscorr)

Specify prefixes

prefixes <- c("treeID", "Hcbh", "dptf", "Hdptf", "effdist", "dist", "Hdist", "maxHcbh", "maxdptf", "maxHdptf", "lastHcbh", "lastdptf", "lastHdptf", "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 matchingcolumns <- grep(paste0("^", prefix), originalcolumn_names, value = TRUE)

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

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

Reorder values

distancesmetricscorr1 <- distancesmetricscorr[, new_order]

Unlist the data frame

distancesmetricscorr2 <- as.data.frame(lapply(distancesmetricscorr1, function(x) unlist(x))) rownames(distancesmetricscorr2) <- NULL head(distancesmetricscorr2) ```

17.Fuels LAD percentage (greater than a threshold)

```{r Fuels LAD percentage for fuel layers with a LAD percentage above a threshold (10 %), echo=TRUE, message=FALSE, warning=FALSE}

library(dplyr) library(magrittr)

LAD profiles derived from normalized ALS data after applying [lad.profile()] function

profilelist2$treeID <- factor(profilelist2$treeID)

Tree metrics derived from geteffectivegap() function

distancesmetricscorr2$treeID <- factor(distancesmetricscorr2$treeID)

treesname1 <- as.character(distancesmetricscorr2$treeID) treesname2 <- factor(unique(trees_name1))

LADmetrics1 <- list() LADmetrics2 <- list()

for (i in levels(treesname2)) { # Filter data for each tree tree1 <- profilelist2 |> dplyr::filter(treeID == i) tree2 <- distancesmetricscorr2 |> dplyr::filter(treeID ==i)

# Get LAD metrics for each tree LADmetrics <- getlayerslad(tree1, tree2, threshold=10, step = 1,minheight= 1.5) LADmetrics1[[i]] <- LADmetrics$df1 LADmetrics2[[i]] <- LADmetrics$df2 }

LADmetricsall1 <- dplyr::bindrows(LADmetrics1) LADmetricsall2 <- dplyr::bindrows(LADmetrics2)

List of data frames

LADmetricslist <- list(LADmetricsall1, LADmetricsall2)

Initialize an empty list to store reordered data frames

fuelsLADmetrics <- list()

Specify prefixes (adjust accordingly)

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

Loop over each data frame

for (i in seqalong(LADmetrics_list)) {

LADmetricsall <- LADmetricslist[[i]]

# Get original column names originalcolumnnames <- colnames(LADmetricsall)

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

# Loop over prefixes for (prefix in prefixes) { # Find column names matching the current prefix matchingcolumns <- grep(paste0("^", prefix), originalcolumn_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 LADmetricsall <- LADmetricsall[, neworder] # Store the reordered data frame in the list fuelsLADmetrics[[i]] <- LADmetricsall } rownames(fuelsLADmetrics[[1]]) <- NULL rownames(fuelsLAD_metrics[[2]]) <- NULL

head(fuelsLADmetrics[[2]]) ```

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

```{r Plots of fuel layers with LAD percentage greater than a threshold, echo=TRUE, message=FALSE, warning=FALSE}

library(ggplot2)

LAD profiles derived from normalized ALS data after applying [lad.profile()] function

profilelist2$treeID <- factor(profilelist2$treeID)

Tree metrics derived from getlayerslad() function

LADgt10p <- fuelsLAD_metrics[[2]]

treesname1 <- as.character(LADgt10p$treeID) treesname2 <- factor(unique(treesname1))

Generate plots for fuels LAD metrics

plotstreesLAD <- getplotseffective(profilelist2, LADgt10p, min_height=1.5)

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

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

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

```{r CBH based on different criteria: maximum LAD, maximum and last distance, echo=TRUE, message=FALSE, warning=FALSE}

library(dplyr) library(magrittr)

Tree metrics derived from getlayerslad() function

LADgt10p <- fuelsLAD_metrics[[2]]

treesname1 <- as.character(LADgt10p$treeID) treesname2 <- factor(unique(treesname1))

cbhmetricslist <- list()

for (j in levels(trees_name2)){

# Filter data for each tree tree1 <- LADgt10p |> dplyr::filter(treeID == j) cbhmetrics <- getcbhmetrics(tree1,minheight= 1.5) cbhmetricslist[[j]] <- cbhmetrics }

Combine depth values for all trees

cbhmetricsall <- dplyr::bindrows(cbhmetrics_list)

Get original column names

originalcolumnnames <- colnames(cbhmetricsall)

# Specify prefixes desiredorder <- c("treeID", "Hcbh", "dptf","effdist","dist", "Hdist", "Hdptf","maxlad","max","last","nlayers")

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

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

# Reorder columns cbhmetricsall <- cbhmetricsall[, new_order]

```

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

```{r Plots of CBH based on different criteria: maximum LAD, maximum and last distance, echo=TRUE, message=FALSE, warning=FALSE} library(ggplot2)

LAD profiles derived from normalized ALS data after applying [lad.profile()] function

profilelist2$treeID <- factor(profilelist2$treeID)

Tree metrics derived from getcbhmetrics() function

cbhmetricsall$treeID <- factor(cbhmetricsall$treeID)

treesname1 <- as.character(cbhmetricsall$treeID) treesname2 <- factor(unique(trees_name1))

Generate plots for fuels LAD metrics

plotscbhmaxlad <- getplotscbhLAD(profilelist2, cbhmetricsall,minheight=1.5) plotscbhmaxdist <- getplotscbhmaxdist(profilelist2, cbhmetricsall,minheight=1.5) plotscbhlastdist <- getplotscbhlastdist(profilelist2, cbhmetricsall,min_height=1.5)

par(mfrow = c(2, 2)) plot(plotscbhmaxlad[[1]]) plot(plotscbhmaxdist[[2]]) plot(plotscbhlastdist[[3]]) ```

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

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

```{r CBH and the LAD percentage below and above the CBH using the breaking point method, echo=TRUE, message=FALSE, warning=FALSE}

library(dplyr) library(magrittr)

LAD profiles derived from normalized ALS data after applying [lad.profile()] function

profilelist2$treeID <- factor(profilelist2$treeID)

Tree metrics derived from getcbhmetrics() function

cbhmetricsall$treeID <- factor(cbhmetricsall$treeID)

treesname1 <- as.character(cbhmetricsall$treeID) treesname2 <- factor(unique(trees_name1))

cumLADmetrics_list <- list()

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

# Get cumulative LAD metrics for each tree cumLADmetricsall <- getcumbreak(tree1, tree2, threshold=75, minheight= 1.5, verbose=TRUE) cumLADmetricslist[[i]] <- cumLADmetricsall }

Combine the individual data frames

cumLADmetrics <- dplyr::bindrows(cumLADmetricslist)

=======================================================================

REORDER COLUMNS

=======================================================================

Get original column names

originalcolumnnames <- colnames(cumLADmetrics)

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 matchingcolumns <- grep(paste0("^", prefix), originalcolumn_names, value = TRUE)

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

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

Reorder columns

cumLADmetrics <- cumLADmetrics[, new_order]

when % LAD is < 75 % below or above the BP (Breaking Point), Hcbh1 is derived from CBH maximum LAD criterium

rownames(cumLADmetrics) <- NULL head(cumLADmetrics) ```

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

```{r Plots of the CBH and the LAD percentage below and above the CBH using the breaking point method, echo=TRUE, message=FALSE, warning=FALSE}

library(ggplot2)

LAD profiles derived from normalized ALS data after applying [lad.profile()] function

profilelist2$treeID <- factor(profilelist2$treeID)

Tree metrics derived from getcumbreak() function

cumLADmetrics$treeID <- factor(cumLADmetrics$treeID)

treesname1 <- as.character(cumLADmetricsall$treeID) treesname2 <- factor(unique(treesname1))

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(plotlist[[1]]) plot(plotlist[[2]]) plot(plot_list[[3]]) ```

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

23. Joinning Fuel ladder properties with Crown polygons

```{r Joining crown polygons and ladder fuels metrics, echo=TRUE, message=FALSE, warning=FALSE}

Tree metrics derived from getlayerslad() function

cbhmetricsall$treeID1 <- factor(cbhmetricsall$treeID1)

crown polygons (output from step 4)

treecrowns$treeID1 <- factor(treecrowns$treeID)

crownsproperties<-merge (treecrowns,cbhmetricsall, by="treeID1") crownsproperties$maxladHcbhfactor <- cut(crownsproperties$maxlad_Hcbh, breaks = 5)

Plotting with a discrete legend

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

ggplot() + geomsf(data = crownsproperties, aes(fill = maxladHcbhfactor)) + scalefillmanual(values = palette(5)) + thememinimal() + labs(title = "Tree Crowns", fill = "maxladHcbh")

```

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 (olga.viedma\@uclm.es{.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.7, accessed on november 2024, 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.