从 Lidr R 中的 LAS 文件计算基于点的特征向量

问题描述 投票:0回答:1

我正在尝试计算 TreeLS 中的 FastPointMetrics 函数使用的 EigenVector12 指标。由于此函数是针对 TLS 数据而设计的,因此它不适用于我正在使用的 ALS 点云。该函数的确切措辞是:“3D 特征向量系数,第二特征向量的第一次加载”。然后我找到了 point_metrics() 函数并能够创建我自己的函数。

我想做的是为每个点的 20 个点的最近邻域获取第二个特征向量的第一次加载,然后取 20m 像素大小的值的平均值,最后得到值的栅格。

这是我写的函数(如有错误请指正):

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

这很有效,只是需要一段时间来处理。然后我找到了 point_eigenvalues() 函数,但这只给了我最大、中和最小的特征值以及其他我不需要的指标。我想知道这个 point_eigenvalues() 函数是否可以执行我需要的操作,或者是否有任何方法可以加速我上面编写的函数。谢谢!

point-clouds eigenvalue eigenvector lidar lidr
1个回答
0
投票

使用C++

point_metrics
的文档明确提到在这种情况下使用C++代码。它为您提供了一个示例,您可以复制并粘贴仅更改输出变量

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 

这快了 10 倍。

使用lidR 4.1.0

使用 github 上提供的

lidR
开发版本(预计将于 1 月 24 日发布),函数
point_eigenvalue()
有一个新参数
coeff

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

基准

使用

library(lidR)
LASfile <- system.file("extdata", "Megaplot.laz", package="lidR")
laz = readLAS(LASfile)
  • EV12_metrics
    :36秒
  • EV12_metrics2
    :3秒
  • point_eigenvalues
    :1.2秒
© www.soinside.com 2019 - 2024. All rights reserved.