Introduction

Here I developed a research methodology using LiDAR data to analyze trees and provided an example of using tree metrics to derive microclimatic from macroclimatic data, showing the potential of LiDAR data in ecological studies.

Needed Packages

Run the following codes to import packages necessary in this project.

library(sf)  # for working with spatial data
library(lidR)  # for working with LiDAR data
library(terra)  # for working with raster data

Get LiDAR data

My LiDAR data used in this case is downloaded on OpenTopography. I want to study the trees in Forest Hill Cemetery in Ann Arbor. This is a screenshot from Google Earth of my research area.

Read and Select LiDAR Data

Use the function lidR::readLAS() to Read the downloaded LiDAR data. It can be in the format of *.las or *.laz. I also set the flag -keep_first to leave only the first returned points so that the superfluous data will be removed at reading time and increase computation speed. We can use the print() function to check the basic information of our data.

fhc_las <- readLAS('../data/lab2/points.laz', filter='-keep_first')  # read LiDAR data
print(fhc_las)  # print information
## class        : LAS (v1.2 format 3)
## memory       : 49 Mb 
## extent       : 274786.5, 275326.3, 4683997, 4684525 (xmin, xmax, ymin, ymax)
## coord. ref.  : WGS 84 / UTM zone 17N 
## area         : 269548 m²
## points       : 675.5 thousand points
## density      : 2.51 points/m²
## density      : 2.51 pulses/m²

Using the function lidR::plot(), we can plot the point clouds easily. To decide which attribute to plot, we can use the View() function and check the attributes in the data.

View(fhc_las)  # show the data table
plot(fhc_las,color='Z',bg='white', axis=T)  # plot the 3D point clouds, set color by elevation

In this plot, there are points within regions I am not interested in. So I will use the function lidR::clip_rectangle() to clip the point clouds. The function range() can be used to decide the range of the coordinates of the LiDAR data.

When doing a clip, it is good practice to include a buffer area to prevent error on the edge of the area in further analyses.

# calculate the original coordinates of the point cloud
x_min <- range(fhc_las$X)[1]  # get the origin of X coordinate
y_min <- range(fhc_las$Y)[1]  # get the origin of Y coordinate
# set the area of interest based on the plotted point cloud
x_aoi <- c(100, 400)  # X range of interest
y_aoi <- c(200, 450)  # Y range of interest
# clip the area of interest
fhc_aoi <- clip_rectangle(fhc_las, x_min+x_aoi[1],y_min+y_aoi[1],x_min+x_aoi[2],y_min+y_aoi[2])
# print the information of the point cloud
print(fhc_aoi)
## class        : LAS (v1.2 format 3)
## memory       : 14.3 Mb 
## extent       : 274886.5, 275186.5, 4684197, 4684447 (xmin, xmax, ymin, ymax)
## coord. ref.  : WGS 84 / UTM zone 17N 
## area         : 76064 m²
## points       : 197.3 thousand points
## density      : 2.59 points/m²
## density      : 2.59 pulses/m²
plot(fhc_aoi,bg='white') # plot the point cloud

Classification

In order to further analyze the trees, classification of ground points is needed. If the original LiDAR data has not been classified, we need to classify the points manually. Common algorithms include PMF, CSF, and MCC, which are included in the lidR package. Luckily, my downloaded LiDAR data has been classified. So there is no need to classify again.

plot(fhc_aoi,color='Classification',bg='white' )  # plot the 3D point cloud, set color by class

As shown in the plot, ground points have been correctly classified.

Get a Digital Terrain model

Get a Digital Terrain Model is usually the second step in processing that follows classification of ground points. Common algorithms include TIN, IDW, and Kriging. Kriging is the slowest but can get a better DTM. Since my data volume is not so huge, I can use Kriging to generate a DTM.

To generate a DTM model with the kriging algorithm we use lidR::rasterize_terrain() where algorithm = kriging()

# get DTM
dtm_kriging <- rasterize_terrain(fhc_aoi, algorithm=kriging(k=20))
# get terrain characteristics for generating a shaded model
dtm_prod <- terrain(dtm_kriging, v = c("slope", "aspect"), unit = "radians")
# compute hill shade
dtm_hillshade <- shade(slope = dtm_prod$slope, aspect = dtm_prod$aspect)
# plot shaded DTM
plot(dtm_hillshade, col =gray(0:30/30), 
     legend =F, xlab='X (m)', ylab='Y (m)', main='Shaded Digital Terrain Model')

The DTM can also be plotted 3d, using the function lidR::plot_dtm3d().

Height Normalization

In order to get a Canopy Height Model (CHM) instead of a Digital Surface Model (DSM), we need to do the height normalization first, so that the derived surface will be representing the relative canopy height instead of the absolute elevations.

Height normalization removes the influence of terrain on above ground measurements. This can be realized either by point cloud based (use the point cloud) or raster based (use the DTM) method. The point cloud based is more accurate since DTM is a discretized raster and the locations of the pixels do not always match the locations of the ground points. Therefore the result ground points will always be \(Z=0\) using the point cloud based method. This method is however computationally intensive.

Since my data volume is not so huge, I will use the point cloud based method.

# do height normalization
nlas <- normalize_height(fhc_aoi, kriging())
# plot the histogram of the ground points
hist(filter_ground(nlas)$Z, breaks=seq(-0.5,0.5,0.01),
     main='Ground Points Height Distribution', xlab='Ground Point Height (m)')

As shown in the histogram, all of the ground points of my LiDAR data have been normalized to zero.

Canopy Height Model

The rasterize_canopy() function uses the ‘local maximum’ or the highest points. The “pit-free” algorithm can be used to develop the CHM. It avoids the empty pixels in the result CHM raster if the grid resolution is set too small when developing the CHM.

# get CHM
chm <- rasterize_canopy(nlas, res = 0.5, 
                        pitfree(thresholds = c(0, 5, 10, 15, 20),  # height thresholds
                                max_edge = c(0, 1)))  # max edge length of the triangles
# plot CHM
plot(chm,col=height.colors(25),
     xlab='X (m)', ylab='Y (m)', main='Canopy Height Model (m)')

Individual Tree Detection and Segmentation

Now we have got the CHM, we can deploy individual tree detection based on that. The simplest way is to use a fix-sized window to detect the highest point on the CHM.

Here I use the function lidR::locate_trees(), where the window size ws is 10 meters. This means I think the radius of the tree canopy within my area is close to 5m.

ttops <- locate_trees(nlas, lmf(ws = 10))  # tree detection using a fixed window size
plot(chm, col = height.colors(50), main='CHM with Tree Tops (m)')  # plot CHM
plot(sf::st_geometry(ttops), add = TRUE, pch = 3)  # add tree centers to the CHM

x <- plot(nlas,bg='white',size=4)  # plot point cloud
add_treetops3d(x,ttops)  # add tree tops 3d

While the sizes of tree canopy can vary based on the ages of trees, it is better to use a variable size window to detect individual trees. Here I developed a simple linear relationship. For research purposes, a more sophisticated model is needed.

# define a window size function
f <- function(x) {
  y <- 7/18*x+20/9
  y[x < 2] <- 3
  y[x > 20] <- 10
  return(y)
}
# plot the function
heights <- seq(-5, 30, 0.5)
ws <- f(heights)
plot(heights, ws, type='l', ylim=c(0,12),
     xlab='Tree Height (m)', ylab='Diameter (m)',
     main='Relationship between Tree Heights and Diameter')

ttops <- locate_trees(nlas,lmf(f))  # get tree tops
plot(chm,col=height.colors(50),main='CHM with Tree Tops (m)')  # plot CHM
plot(sf::st_geometry(ttops),add=T,pch=3)  # plot tree tops

x <- plot(nlas,bg='white',size=4)  # plot point cloud
add_treetops3d(x,ttops)  # plot tree tops

In order to derive tree metrics for further analyses, we need the point cloud data to include tree segmentation information. Here I use the algorithm silva2016 to segment trees because it can get the best result on my data. To use other algorithms, do ?segment_trees to check the document.

algo <- silva2016(chm,treetops=ttops)  # set algorithm
las_seg = segment_trees(nlas,algo)  # deploy tree segmentation
plot(las_seg,bg='white',color='treeID')  # plot point cloud, set color by tree segmentation result

Derive Tree Metrics

Analyses of point cloud data are often based on metrics calculations. Metrics are scalar summaries of point distributions that can be computed using varying neighborhood definitions and varying reference locations. The notion of metrics is at the core of the lidR package, which enables the computation of standard or custom user-defined metrics at varying levels of regularization.

# get tree metrics
metrics <- crown_metrics(las_seg,
                         ~list(z_max = max(Z), z_mean = mean(Z)),  # set function
                         geom='concave')  # set the crowns' shape
head(metrics)
## Simple feature collection with 6 features and 3 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 274886.5 ymin: 4684258 xmax: 274912 ymax: 4684387
## Projected CRS: WGS 84 / UTM zone 17N
##   treeID z_max   z_mean                       geometry
## 1      1 25.18 17.65179 POLYGON ((274893.1 4684348,...
## 2      2 21.89 14.09292 POLYGON ((274897.2 4684259,...
## 3      3 12.87  6.87518 POLYGON ((274903.7 4684287,...
## 4      4 24.19 15.08976 POLYGON ((274909.7 4684343,...
## 5      5 22.24 14.72880 POLYGON ((274906.2 4684355,...
## 6      7 24.50 16.91201 POLYGON ((274911.8 4684382,...
plot(metrics['z_max'],pal=hcl.colors,pch=19,main='Tree Height Map (m)')

Application

Using the tree metrics, we can do a lot of things including selecting the trees or doing tree based inventory. One interesting application I recently read is using the tree metrics to derive microclimate from macrocliamte, since trees have an effect on meterology indexes like temperature.

For example, I have a land surface temperature raster data which can be computed from Landsat imagery. Here I generate a raster data from CHM to represent the temperature data which can be calculated from Landsat imagery. Everywhere in my study area is assigned 35 degree Celsius.

temp=setValues(chm,35)  # 35 degrees Celsius

Let’s say the trees have a buffering effect on temperature, and it is decided by the height and intensity of the cloud points with a relationship:

\(CE=0.01 \times Z_{max} + 0.0002 \times I_{mean}\)

Where \(CE\) stands for the cooling effect in degrees Celsius (how many degrees Celsius the temperature will decrease), \(Z_{max}\) is the max height of a tree point cloud, and \(I_{mean}\) is the mean intensity of the tree point cloud.

ce_metrics <- crown_metrics(las_seg, ~list(zmax=max(Z), imean=mean(Intensity)))  # get tree metrics
ce_metrics$cooling_effect <- 0.01*ce_metrics$zmax + 0.0002*ce_metrics$imean  # get cooling effect index
# rasterize the cooling effect map
r <- rast(ext(las_seg),resolution=10)
v <- vect(ce_metrics["cooling_effect"])
ce_map <- terra::rasterize(v, r, field = "cooling_effect", fun = mean)
# replace the NA cells with 0
ce_map[is.na(ce_map[])] <- 0
# set the map Coordinate System
crs(ce_map) <- crs(temp)
# plot the cooling effect map
plot(ce_map, col = hcl.colors(15),
     xlab='X (m)',ylab='Y (m)',main='Cooling Effect (degrees Celsius)') # plot the cooling effect map

Then we can do the raster computation and get the microclimate.

# resample the temperature map to the same resolution of the cooling effect map
temp <- resample(temp,ce_map,method='bilinear')
# do raster computation
micro_climate <- temp-ce_map
plot(micro_climate,col = hcl.colors(10),
     xlab='X (m)',ylab='Y(m)',main='Temperature Map (degrees Celsius)')

Cool! However, this is a rather simplified model. For research purposes, a more sophisticated model is needed.