Forestry Analysis Using R and LiDAR Data
Tutorial: Use R Language in RStudio and LAS Data to Complete Various Forestry Analyses
Using R language in RStudio, and LiDAR point clouds, you are able to complete three distinct types of forestry analyses: Individual Tree Detection, Individual Tree Segmentation, and the creation of Canopy Height Models. These analyses are important to make data-driven decisions in forestry management, they provide tools to support sustainable forest management by providing precise information on forest composition, structure, and dynamics, they help us in making environmental decisions to conserve biodiversity and protect ecosystems, and enhance the economic efficiency of timber harvesting and other forestry operations.
Individual Tree Detection (ITD) is the process of identifying and locating individual trees within a forest using LiDAR. It involves detecting the position of each tree’s highest point, also known as the apex, within the forest canopy. ITD is important for forest inventories essential for sustainable forest management, biodiversity assessments to monitor the species diversity and health of specific trees, resource management to help estimate timber volume and assess economic value of forest resources, and environmental monitoring where forest structure can be tracked overtime to understand the ecological pressures and impact of environmental changes.
Individual Tree Segmentation (ITS) involves not only detecting individual trees but also delineating their crowns from the LiDAR point cloud data. This is important for detail analysis of the tree characteristics (i.e. the crown shape, size, and structure); forest health monitoring to detect any signs of disease, pest infections or damage; species identification to help with biodiversity studies and conservation efforts, and precision forestry which enables targeted interventions and optimizes forestry management activities.
Canopy Height Models (CHM) are a digital representation of the height of the forest canopy above the ground, and is created by subtracting the Digital Terrain Model (the bare Earth surface) from the Digital Surface Model (elevation of both ground and any above ground features). The result is a map of the tree heights. CHM’s are important for forest structure analysis, biomass and carbon storage estimation which are important for understanding carbon cycles and climate change mitigation, habitat mapping as various wildlife species depend on specific forest structures, and fire risk assessment to identify areas with tall, dense canopies that are more prone to wildfires.
This tutorial goes through the following steps:
Step 1: Source LiDAR Data
If you already have a LiDAR point cloud, jump down to Step 3. If you need to source a LiDAR point cloud, visit Equator and follow the steps below to download your LAZ file.
- Sign in or create an account with Equator.
- Search for your site, then select the +NEW SITE button at the bottom of the screen.
- Edit your site bounds and select okay.
- Navigate to the Data tab and choose ‘Point Cloud’. Choose your parameter and select ‘Generate’.
- Your point cloud is stored in the Layers tab under your site.
- To download, select the download button (down arrow). Choose your parameters and select ‘Process’.
- The .laz file of your site will appear in your Downloads folder (or whichever folder you have directed your downloads to go).
RStudio & LiDAR: LiDAR Point Cloud from Equator
Step 2: Convert LAZ to LAS
If you downloaded your data from Equator, it downloaded as an LAZ file. An LAZ file is a compressed version of the LAS File. Some programs will accept LAZ files; however, other programs such as RStudio will not. You will need to convert your LAZ file to an LAS file. There are many ways you can do this. We prefer Cloud Compare and have a whole tutorial on how do do this here.
If you already had LiDAR point cloud data and it was in the form of an .las extension, jump down to Step 3.
Step 3: Set up your Project in RStudio
The first step in setting up your project is installing and loading the R Packages that are useful for working with LAS data. We recommend the lidR package, which is commonly used for LiDAR data processing in R. To download the lidR package, click here.
Next, read the LAS dataset that you have sourced. Use the readLAS function from the lidR package to read the LAS data into R.
The code in RStudio should read:
install.packages(“lidR“))
library(lidR)
las <- readLAS(“C:/Users/…/file.las”, select = “xyzr“, filter = “-drop_z_below 0“)
chm <- rasterize canopy(las, 0.5, pitfree(subcircle = 0.2)
plot(las, bg = “white”, size = 4)
plot(sf::st_geometry(ttops), add =TRUE, pch = 3)
add_treetops3d(x, ttops)
las <- readLAS(“C:/Users/…/file.las”)
ttops <- locate_trees(las, lmf(ws = 30))
plot(chm, col = height.colors(50))
x <- plot(las, bg = “white”, size = 4)
Step 4: Individual Tree Detection (ITD) and Segmentation (ITS)
In this section, the application of Individual Tree Detection (ITD) and Individual Tree Segmentation (ITS) will be carried out. The LMF (Local Maximum Filter) can be applied with a constant size windows of ws = 30 meters meaning that for a given point the algorithm looks to the neigbourhood points within a 2.5 radius circle to figure out if the point is the local highest.
The code is as follows:
Part 1: Read LAS file
las <- readLAS(“C:/Downloads/USGS_LPC_IN_Central_Hamilton_2017_LAS_2019.las”, filter = “-drop_z_below 0”)
las <- segment_trees(las, li2012())
metrics <- crown_metrics(las, ~list(z_max = max(Z), z_mean = mean(Z)))
head(metrics)
plot(metrics[“z_max”], pal = hcl.colors, pch = 19)
colnames(metrics)
RStudio & LiDAR: Plotting the maximum canopy height (z_max) using color coding
Part 2: Define the Custom Metrics Function
A custom function is defined to calculate the maximum height (z_max), standard deviation of heights (z_sd), mean intensity (i_mean), and maximum intensity (i_max). The custom crown metrics (ccm) calculates the mean intensity (imean) for each tree crown.
custom_crown_metrics <- function(z, i) { # user-defined function
metrics <- list(
z_max = max(z), # max height
z_sd = sd(z), # vertical variability of points
i_mean = mean(i), # mean intensity
i_max = max(i) # max intensity
)
return(metrics) # output
}
ccm = ~custom_crown_metrics(z = Z, i = Intensity)
metrics <- crown_metrics(las, func = ccm, geom = “convex”)
plot(metrics[“z_max”], pal = hcl.colors)
Part 3: Calculate the metrics
This step calculates, then filters out the tree crowns with mean intensities below 10000.
metrics <- crown_metrics(las, ~list(imean = mean(Intensity))) # calculate tree intensity metrics
metrics <- metrics[metrics$imean > 10000,] # filter intensity
subset <- filter_poi(las, treeID %in% metrics$treeID)
x <- plot(las, bg = “white”, size = 4)
plot(subset, add = x + c(-100, 0), size = 5) # some plotting
Part 4: Calculates the metrics data for the tree based inventory
This step converts the metrics data to a raster format and plots it, summing all the G values over a 15m resolution grid.
metrics <- crown_metrics(las, func = ~custom_crown_metrics(Z, Intensity)) # calculate intensity metrics
metrics$G <- 1.7 * metrics$z_max + 1.1 * metrics$i_mean # set value of interest
plot(metrics[“G”], pal = hcl.colors, pch = 19) # some plotting
r <- terra::rast(ext(las), resolution = 15)
v <- terra::vect(metrics[“G”])
map <- terra::rasterize(v, r, field = “G”, fun = sum) # extract sum of G at 15m
plot(map, col = hcl.colors(15)) # some plotting
RStudio & LiDAR: Convert metric data to raster format
RStudio & LiDAR: Summing the G values over a 15m resolution grid
Part 5: Tree Based Segmentation & Canopy Height Model
This step creates a canopy height model with a resolution of 0.5m. It locates tree tops using a local maximum filter with a window size of 45, then plots the canopy height model and tree tops.
chm <- rasterize_canopy(las, 0.5, pitfree(subcircle = 0.2))
plot(las, bg = “white”, size = 4)
ttops <- locate_trees(las, lmf(ws = 45))
plot(chm, col = height.colors(50))
plot(sf::st_geometry(ttops), add = TRUE, pch = 3)
x <- plot(las, bg = “white”, size = 4)
add_treetops3d(x, ttops)
RStudio & LiDAR: Create a Canopy Height Model (CHM) with a
resolution of 0.5m
RStudio & LiDAR: Create a plot of the representative trees