Introduction

Lately, several works have been conducted aiming to investigate possibilities of extracting forests’ structural metrics from Digital Aerial Photogrammetry (DAP) point clouds (PC).
This markdown demonstrates a simple workflow on lidR for DAP-PC processing and tree metrics deriving.

Loading Data

roi <- st_read(path_roi, quiet =T)
point_cloud <- readLAS(path_point_cloud) %>% clip_roi(roi)
print(point_cloud)
#> class        : LAS (v1.2 format 2)
#> memory       : 115.7 Mb 
#> extent       : 306694.5, 306842.1, 8353199, 8353330 (xmin, xmax, ymin, ymax)
#> coord. ref.  : SIRGAS 2000 / UTM zone 24S 
#> area         : 9960.452 m²
#> points       : 1.44 million points
#> density      : 144.9 points/m²

Point Cloud Visual Inspection

plot(point_cloud,color="RGB",size = 3, bg = "white",  axis = TRUE)

Creating a Digital Surface Model (DSM)

dsm <- grid_canopy(point_cloud,res=0.5,algorithm = dsmtin())
plot(dsm,col = height.colors(50))

Ground Classification

ws <- seq(3, 12, 3)
th <- seq(0.1, 1.5, length.out = length(ws))
point_cloud_class <- classify_ground(point_cloud, algorithm = pmf(ws = ws, th = th))

Noise Classification

point_cloud_class <- classify_noise(point_cloud_class, ivf(res = 0.5, n = 50))
plot(point_cloud_class %>% filter_poi(Classification!=LASNOISE),
     color = "Classification", size = 3, bg = "white", axis = TRUE)

Creating a Digital Terrain Model (DTM)

dtm_tin <- grid_terrain(filter_ground(point_cloud_class), 0.5, tin(),keep_lowest = T)
plot_dtm3d(dtm_tin, bg = "white")

Normalizing Point Cloud

point_cloud_norm <- normalize_height(point_cloud_class, knnidw())%>% filter_poi(Z>= 0)
#> Warning: There were 6681 degenerated ground points. Some X Y Z coordinates were
#> repeated. They were removed.
plot(point_cloud_norm, color = "Z", size = 3, bg = "white", axis = TRUE)

Analysing a Vertical Transect

p1 <- c(306775, 8353310)
p2 <- c(306732, 8353250)

las_tr <- clip_transect(point_cloud_norm, p1, p2, width = 4, xz = TRUE)

ggplot(las_tr@data, aes(X,Z, color = Z)) + 
  geom_point(size = 0.5) + 
  coord_equal() + 
  theme_minimal() +
  scale_color_gradientn(colours = height.colors(50))

Create a Canopy Height Model (CHM)

chm <- grid_canopy(point_cloud_norm,res=0.5,algorithm = dsmtin(max_edge = 3))
plot(chm,col = height.colors(50))

Individual Tree Detection

ttops <- find_trees(point_cloud_norm, lmf( ws = 5))
plot(chm, col = height.colors(50))
plot(ttops, add = TRUE)

x <- plot(point_cloud_norm, bg = "white", size = 4, axis = TRUE, legend = TRUE)
add_treetops3d(x, ttops)

Individual Tree Segmentation

algo <- dalponte2016(chm, ttops)
point_cloud_norm <- segment_trees(point_cloud_norm, algo) 
plot(point_cloud_norm, bg = "white", size = 3, color = "treeID") 

plot(filter_poi(point_cloud_norm, treeID == 5), size = 3, bg = "white",  axis = TRUE, legend = TRUE)

Individual Tree Metrics

ttops <- tree_metrics(point_cloud_norm, func = .stdmetrics) %>% st_as_sf()
kable(summary(ttops[,c( "zmax", "zmean","zsd","zq95","area")]%>% st_drop_geometry()))
zmax zmean zsd zq95 area
Min. : 5.004 Min. :2.376 Min. :0.6396 Min. : 4.463 Min. : 0.1772
1st Qu.: 8.679 1st Qu.:4.472 1st Qu.:2.6744 1st Qu.: 8.214 1st Qu.:19.3972
Median : 9.528 Median :5.250 Median :2.9570 Median : 8.890 Median :35.2686
Mean : 9.440 Mean :5.137 Mean :2.9920 Mean : 8.784 Mean :31.9503
3rd Qu.:10.157 3rd Qu.:5.795 3rd Qu.:3.4776 3rd Qu.: 9.663 3rd Qu.:46.1242
Max. :11.792 Max. :7.338 Max. :4.1238 Max. :10.751 Max. :59.4550 
hist(ttops$zmax,plot = T, main = "Height Distribruition (m)",xlab =   "Height (m)",ylab = "Frequency")