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.

  1. Sign in or create an account with Equator.
  2. Search for your site, then select the +NEW SITE button at the bottom of the screen.
  3. Edit your site bounds and select okay.
  4. Navigate to the Data tab and choose ‘Point Cloud’. Choose your parameter and select ‘Generate’.
  5. Your point cloud is stored in the Layers tab under your site.
  6. To download, select the download button (down arrow). Choose your parameters and select ‘Process’.
  7. The .laz file of your site will appear in your Downloads folder (or whichever folder you have directed your downloads to go).
RStudio & LiDAR: Equator Point Cloud

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 height (z_max) using colour coding

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: Converting metric data to raster format

RStudio & LiDAR: Convert metric data to raster format

RStudio & LiDAR: Plotting G value over a 15m resolution grid

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 canopy height model

RStudio & LiDAR: Create a Canopy Height Model (CHM) with a
resolution of 0.5m

RStudio & LiDAR: Plot CHM of representative trees

RStudio & LiDAR: Create a plot of the representative trees

Frequently Asked Questions

R is a programming language and open-source software environment commonly used for statistical computing and data analysis. It provides a wide variety of statistical and graphical techniques, and it is extensible, allowing users to add new functions through packages. R was created by Ross Ihaka and Robert Gentleman at the University of Auckland, New Zealand, and it is part of the GNU project.

There are a number of key applications the R language can be used for. These include:

  1. Statistical Analysis: R provides a comprehensive set of statistical and mathematical functions, making it popular among statisticians and data analysts.
  2. Data Visualization: R has powerful tools for creating various types of graphs and visualizations, such as scatter plots, bar charts, histograms, and more.
  3. Data Manipulation: R has extensive capabilities for manipulating and transforming data. It supports data frames, which are similar to tables in databases, making it easy to work with structured data
  4. Extensibility: R is highly extensible through packages. There is a vast repository of packages available for various purposes, and users can create their own packages to share functionality
  5. Open Source: R is free and open-source, which means that anyone can access, use, and modify the source code.
  6. Community Support: R has a large and active community of users and developers. This community contributes to the development of packages, provides support through forums, and shares resources.

RStudio is an integrated development environment (IDE) specifically designed for working with the R programming language. It provides a user-friendly interface and a set of powerful tools to enhance the workflow of R users. RStudio is widely used by statisticians, data scientists, and analysts for data analysis, visualization, and report generation.

  1. Script Editor: RStudio includes a script editor with syntax highlighting, code completion, and other features that make it easier to write and edit R code.
  2. Console: RStudio provides an interactive R console where users can directly execute R code and see the output in real-time.
  3. Workspace: It includes a workspace pane that displays the current R environment, allowing users to view and manage variables, data frames, and other objects.
  4. Plots and Graphs: RStudio has integrated support for creating and viewing plots and graphs directly within the IDE. This feature is useful for exploratory data analysis.
  5. Packages: RStudio provides tools for managing and installing R packages, which extend the functionality of R. Users can easily install, update, and load packages from the IDE.
  6. Integrated Help: The IDE includes an integrated help pane that provides documentation and assistance on R functions and packages, making it easier for users to find information.
  7. Version Control: RStudio has built-in support for version control systems like Git and SVN, making it convenient for users to manage their code repositories.
  8. Markdown Support: RStudio supports the creation of dynamic documents using R Markdown, allowing users to integrate R code with narrative text, create interactive documents, and generate reports.

R Studio can be downloaded from the Posit website for free.

Check out this video on how to download LAZ data from Equator.

The .las is a file extension used for the Lidar Data Exchange Format, and is the file format used to store point cloud data. It contains information about each point collected, including x, y, and z coordinates, the classification of the point and the projection.

When you download a point cloud from Equator, it exports as an LAZ file. An LAZ file is just a compressed version of the LAS file. When bringing a point cloud into R Studio, as well as many other programs, an LAS file is required, which means you will need to convert the file. We have a whole written and video tutorial on how to use Cloud Compare to convert an LAZ file to an LAS file here.