💾 Archived View for republic.circumlunar.space › users › johngodlee › posts › 2020-12-30-height.gmi captured on 2023-03-20 at 18:58:29. Gemini links have been rewritten to link to archived content

-=-=-=-=-=-=-

Extracting a vertical height profile from terrestrial LiDAR

DATE: 2020-12-30

AUTHOR: John L. Godlee

Previously, measuring canopy cover with hemispherical photography only provided a 2D representation of the canopy, but with LiDAR it's possible to measure variation in canopy cover over the height of the canopy to create a canopy height profile. Here I want to describe how I used R to process the XYZ point cloud data to create a canopy height profile. I have already described in a previous post how I voxelise, clean and crop the point cloud, using PDAL[1].

1: https://pdal.io/

# Packages
library(ggplot2)
library(dplyr)
library(data.table)
library(scico)
library(zoo)

I used data.table::fread() to read the XYZ point cloud .csv files into R, as they are very large, about 500 MB, and fread() seems to do a better job at reading large files into memory.

For each file, I rounded the elevation (Z) coordinates to the nearest cm, then for each cm height bin I calculated the volume of space occupied by voxels.

I created a height foliage density profile with ggplot().

I calculated the effective number of layers according to Ehbrecht et al. 2016 (Forest Ecology and Management), which basically splits the height profile into 1 m bins and calculates the Shannon diversity index of the foliage volume occupied in each layer. Here is the function for it:

#' Effective number of layers in a point cloud distribution
#'
#' @param x vector of Z (elevation) coordinates 
#' @param binwidth width of vertical bins in units of x
#'
#' @return atomic vector of length one describing the effective number of layers
#'     in the canopy
#'
#' @details Uses the Shannon diversity index (Entropy) to estimate the 
#'     "Effective Number of Layers" in the vertical profile of a point cloud 
#'     distribution.
#'
#' @references
#' Martin Ehbrecht, Peter Schall, Julia Juchheim, Christian Ammer, & 
#'     Dominik Seidel (2016). Effective number of layers: A new measure for 
#'     quantifying three-dimensional stand structure based on sampling with 
#'     terrestrial LiDARForest Ecology and Management, 380, 212–223.
#' 
#' @examples 
#' x <- rnorm(10000)
#' enl(x)
#' 
# Calculate effective number of layers in canopy
## Assign to Z slices
## Count number of points within each slice
## Calculate shannon diversity index (entropy) on vertical layer occupancy
enl <- function(x, binwidth) { 
    binz <- cut(x, include.lowest = TRUE, labels = FALSE,
        breaks = seq(floor(min(x)), ceiling(max(x)), by = binwidth))

    n <- unlist(lapply(split(x, binz), length))

    entropy <- exp(-sum(n / sum(n) * log(n / sum(n))))

    return(entropy)
}

I calculated the area under the curve of the foliage density profile using density() then zoo::rollmean(), a method I stole of Stack Overflow.

I also calculated the height above the ground of the peak of foliage density.

Here is the script in its entirety:

# Import data
file_list <- list.files(path = "../dat/tls/height_profile", pattern = "*.csv", 
  full.names = TRUE)

# Check for output directories
hist_dir <- "../img/foliage_profile"
if (!dir.exists(hist_dir)) {
  dir.create(hist_dir, recursive = TRUE)
}

out_dir <- "../dat/subplot_profile"
if (!dir.exists(out_dir)) {
  dir.create(out_dir, recursive = TRUE)
}

# Define parameters 
voxel_dim <- 0.01
z_width <- 1
cylinder_radius <- 10

# Calculate maximum 1 voxel layer volume
layer_vol <- pi * cylinder_radius^2 * voxel_dim

# For each subplot:
profile_stat_list <- lapply(file_list, function(x) {

  # Get names of subplots from filenames
  subplot_id <- gsub("_.*.csv", "", basename(x))
  plot_id <- gsub("(^[A-Z][0-9]+).*", "\\1", subplot_id)
  subplot <- gsub("^[A-Z][0-9]+(.*)", "\\1", subplot_id)

  # Read file
  dat <- fread(x)

  # Round Z coords to cm
  dat$z_round <- round(dat$Z, digits = 2)

  # Calculate volume and gap fraction
  bin_tally <- dat %>% 
    group_by(z_round) %>%
    filter(z_round > 0) %>%
    tally() %>% 
    as.data.frame() %>%
    mutate(vol = n * voxel_dim,
      gap_frac = vol / layer_vol)

  # Plot gap fraction density plot 
  pdf(file = paste0(hist_dir, "/", subplot_id, "_foliage_profile.pdf"), 
    width = 8, height = 6)
    print(
      ggplot(bin_tally, aes(x = z_round, y = gap_frac)) +
        geom_line() +
        theme_bw() + 
        labs(x = "Elevation (m)", y = "Gap fraction") + 
        coord_flip()
    )
  dev.off()

  # Calculate effective number of layers
  layer_div <- enl(dat$Z, z_width)

  # Calculate area under curve 
  den <- density(dat$z_round)

  den_df <- data.frame(x = den$x, y = den$y)

  auc_canopy <- sum(diff(den_df$x) * rollmean(den_df$y, 2))

  # Calculate height of max peak
  dens_peak_height <- den_df[den_df$y == max(den_df$y), "x"]

  # Create dataframe from stats
  out <- data.frame(plot_id, subplot, layer_div, auc_canopy, dens_peak_height)

  # Write to file
  write.csv(out,
    file.path(out_dir, 
      paste0(paste(plot_id, subplot, sep = "_"), "_summ.csv")),
    row.names = FALSE)

  return(out)
})

Foliage height density profile