Calculating Point-Based EigenVectors from a LAS file in Lidr R

Question

I am trying to calculate the EigenVector12 metric that the FastPointMetrics function in TreeLS uses. Since this function is made for TLS data, it does not work for the ALS point cloud I am using. The exact wording for the function is: "3D eigenvector coefficients, 1st load of 2nd eigenvector". I then found the point_metrics() function and was able to make my own function.

What I am trying to do is get the 1st loading of the 2nd eigenvector for the nearest neighborhood of 20 points for every point, and then take the mean of the values in a 20m pixel size ending up with a raster of the values.

This is the function I wrote (correct me if it is wrong):

EV12_metrics = function(x,y,z, th1 = 20) {
  xyz <- cbind(x,y,z)         # create XYZ matrix
  cov_m <- cov(xyz)           # compute covariance of the point locations
  e = eigen(cov_m)            # get eigenvalues and eigenvectors from the matrix
  vect12 = e$vectors[4]       # Select the 1st load of the second eigenvector
  return(list(EV12 = vect12)) # return with a list of values 
}

M = point_metrics(laz, ~EV12_metrics(X,Y,Z), k = 20)   # apply the function to the point cloud 
laz = add_attribute(laz, M$EV12, "EV12")               # Add the attribute to the LAS points
gridmets = grid_metrics(laz, mean(EV12), res = 20)     # Create raster of the mean of the values at 20m resolution

This worked well, except it took a while to process. Then I found the point_eigenvalues() function, but that gave me only the largest, medium, and smallest eigenvalues and other metrics that I did not need. I am wondering if this point_eigenvalues() function can do what I need it to or if there is anyway to speed up the function I wrote above. Thanks!

Answer 1

Use C++

The documentation of point_metrics explicitly mention to use C++ code in this case. It gives you an example you can copy and paste changing only the output variable

Rcpp::sourceCpp(code = "
#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]

// [[Rcpp::export]]
SEXP EV12_metrics_cpp(arma::mat A) {
arma::mat coeff;
arma::mat score;
arma::vec latent;
arma::princomp(coeff, score, latent, A);
return(Rcpp::wrap(coeff));
}")

EV12_metrics2 = function(x,y,z, th1 = 25, th2 = 6) {
  xyz <- cbind(x,y,z)
  e <- EV12_metrics_cpp(xyz)
  vect12 = e[4]  
  return(list(EV12 = vect12)) # return with a list of values 
}


M2 = point_metrics(laz, ~EV12_metrics2(X,Y,Z), k = 20)   # apply the function to the point cloud 

This is 10 times faster.

Use lidR 4.1.0

Using the development version of lidR available on github (will be released in Jan 2024 hopefully) the function point_eigenvalue() has a new parameter coeff

M3 = point_eigenvalues(las, 20, coeffs = T)
M3 = M3$coeff01

Benchmark

Using

library(lidR)
LASfile <- system.file("extdata", "Megaplot.laz", package="lidR")
laz = readLAS(LASfile)

• EV12_metrics: 36 s
• EV12_metrics2: 3 s
• point_eigenvalues: 1.2s