Classification of snag points using lassnags

Document created on February 6th 2017 by @bi0m3trics, updated on February 13th 2020 and up-to-date with lidR v3.0.0

The following example shows how to classify points in a lidar dataset (using segment_snags) based on their intensity (e.g., 8-bit, 0-255 values) using the algorithm developed by Wing et al. (2014).

Read in and filter the file

First, we need to read the file (and mixed conifer example that comes with the lidR package). For this we only need the spatial coordinates and intensity and in accordance with Wing et al. (2014) we should only need the first and single returns but here we’ll just keep the first returns. We’re assuming the data are normalized, and thus don’t include the classification attribute; however, if the Z values are in absolute elevation you’ll need to normalize them first.

library(lidR)
LASfile <- system.file("extdata", "MixedConifer.laz", package="lidR")
las = readLAS(LASfile, select = "xyzi", filter="-keep_first") 
plot(las)

Which results in the following las data being loaded

Algorithm overview

Figure 2 from Wing et al (2014) outlines his approach, which is described in full in the article's text. In short, the algorithm uses the four following steps:

1. Calculate plot-level lidar variables which are used to train the sensitivity of the algorithm. These are arguments passed to the segment_snags function call.
2. Calculate individual point neighborhood intensity and point density statistics for three separate 3D neighborhood variables (figure below).
3. Use these variables in a conditional assessment framework (decision tree) to determine whether an individual point is associated with a live tree or snag. Note: This example stops here!
4. Lastly, Wing et al’s approach takes all overstory points located within a 1m radius cylinder of the snag classified points and classifies them as snag points. As an additional step, they then take all points not receiving a snag classification (i.e., live tree points) and eliminate them by replacing their height values with a zero value to create the final snag-filtered point cloud.

The following is Figure 4 from Wing et al (2014) describing the point neighborhood intensities used in step two of the algorithm.

Preparing and running segment_snags with Wing's algorithm

The following are the parameters that the wing2015 algorithm needs to run:

• neigh_radii – this is a vector of three radii used in quantifying local-area centered neighborhoods. The defaults are 1.5, 1, and 2 for the sphere, small cylinder and large cylinder neighborhoods, respectively.
• low_int_thrsh – This is the lower intensity threshold filtering value. As mentioned above, these should be based on apriori knowledge (plot-level lidar data) The default is 50.
• uppr_int_thrsh – This is the upper intensity threshold filtering value. The default is 170.
• pt_den_req – This is the point density requirement based on plot-level point density defined classes. See Wing et al. (2015) page 172 fir it’s description and the segment_snags default is 3.
• BBPRthresh_mat – This is a 3x4 matrix which provides the four average BBPR (branch and bole point ratio) values for each of the three neighborhoods (sphere, small cylinder and large cylinder) to be used for conditional assessments and classification into the following four snag classes: (1) general snag (2) small snag (3) live crown edge snag (4) high canopy cover snag. See Wing et al. (2015) page 172 and Table 2. This matrix must be provided by the user.

Here's the code for our example to parameterize and call segment_snags:

BBPRthrsh_mat <- matrix(
c(0.80, 0.80, 0.70, 0.85, 
  0.85, 0.60, 0.80, 0.80, 
  0.60, 0.90, 0.90, 0.55),
  nrow = 3L, ncol = 4L)

algorithm <- wing2015(neigh_radii = c(1.5, 1, 2), low_int_thrsh = 50, uppr_int_thrsh = 170, pt_den_req = 3, BBPRthrsh_mat = BBPRthrsh_mat)
las <- segment_snags(las, algorithm)

Which results in the the point being classifies into non-snag and snag points.

Note: that this algorithm strictly performs a classification based on user input while the original publication's methods also included a segmentation step and some pre- (filtering for first and single returns only) and post-process (filtering for only the snag classified points prior to segmentation) tasks which are now expected to be performed by the user.

This algorithm attributes each point in the point cloud (snagCls column) into the following five snag classes:

0. live tree - not a snag
1. general snag - the broadest range of snag point situations
2. small snag - isolated snags with lower point densities
3. live crown edge snag - snags located directly adjacent or intermixing with live trees crowns
4. high canopy cover snag - snags protruding above the live canopy in dense conditions (e.g., canopy cover >= 55%).

Examining the results

Once we’ve classified the points, we can plot the results using the newly defined snagCls field in the lidar data using our own predefined color palette…

colorPal <- c("white", "red", "orange", "yellow", "light green")
plot(las, color = "snagCls", colorPalette = colorPal) 

Lastly, we can create a new las object of only snags points. This new object could then be segmented using segment_trees and/or analyzed using a custom metrics function to extract snag attributes using tree_metrics.

snags <- filter_poi(las, snagCls > 0)
colorPal <- c("red", "orange", "yellow", "light green")
plot(snags, color = "snagCls", colorPalette = colorPal)