lidR (version 2.2.5)

point_metrics: Point-based metrics

Description

Computes a series of user-defined descriptive statistics for a LiDAR dataset for each point. This function is very similar to grid_metrics but computes metrics for each point based on its k-nearest neighbours or its sphere neighbourhood.

Usage

point_metrics(las, func, k, r, xyz = TRUE, filter = NULL)

Arguments

las

An object of class LAS

func

formula. An expression to be applied to each cell (see section "Parameter func").

k

integer. k-nearest neighbours

r

numeric. radius of the neighborhood sphere (not supported yet).

xyz

logical. Coordinates of each point are returned in addition to each metric. If filter = NULL coordinates are references to the original coordinates and do not occupy additional memory. If filter != NULL it obviously takes memory.

filter

formula of logical predicates. Enables the function to run only on points of interest in an optimized way. See also examples.

Parameter <code>func</code>

The function to be applied to each cell is a classical function (see examples) that returns a labeled list of metrics. For example, the following function f is correctly formed.

f = function(x) {list(mean = mean(x), max = max(x))}

And could be applied either on the Z coordinates or on the intensities. These two statements are valid:

point_metrics(las, ~f(Z), k = 8)
point_metrics(las, ~f(Intensity), k = 5)

Everything that works in grid_metrics should also work in point_metrics but might be meaningless. For example, computing the quantile of elevation does not really makes sense here.

Details

It is important to bear in mind that this function is very fast for the feature it provides i.e. mapping a user-defined function at the point level using optimized memory management. However, it is still computationally demanding. To help users to get an idea of how computationally demanding this function is, let's compare it to grid_metrics. Assuming we want to apply mean(Z) on a 1 km<U+00B2> tile with 1 point/m<U+00B2> with a resolution of 20 m (400 m<U+00B2> cells), then the function mean is called roughly 2500 times (once per cell). On the contrary, with point_metrics, mean is called 1000000 times (once per point). So the function is expected to be more than 400 times slower in this specific case (but it does not provide the same feature). This is why the user-defined function is expected to be well optimized, otherwise it might drastically slow down this already heavy computation. See examples. Last but not least, grid_metrics() relies on the data.table package to compute a user-defined function in each pixel. point_metrics() relies on a similar method but with a major difference: is does not rely on data.table and thus has not been tested over many years by thousands of people. Please report bugs, if any.

See Also

Other metrics: cloud_metrics(), grid_metrics(), hexbin_metrics(), tree_metrics(), voxel_metrics()

Examples

# NOT RUN {
LASfile <- system.file("extdata", "Topography.laz", package="lidR")

# Read only 0.5 points/m^2 for the purposes of this example
las = readLAS(LASfile, filter = "-thin_with_grid 2")

# Computes the eigenvalues of the covariance matrix of the neighbouring
# points and applies a test on these values. This function simulates the
# 'shp_plane()' algorithm from 'lasdetectshape()'
plane_metrics1 = function(x,y,z, th1 = 25, th2 = 6) {
  xyz <- cbind(x,y,z)
  cov_m <- cov(xyz)
  eigen_m <- eigen(cov_m)$value
  is_planar <- eigen_m[2] > (th1*eigen_m[3]) && (th2*eigen_m[2]) > eigen_m[1]
  return(list(planar = is_planar))
}

# Apply a user-defined function
M <- point_metrics(las, ~plane_metrics1(X,Y,Z), k = 25)
#> Computed in 6.3 seconds

# We can verify that it returns the same as 'shp_plane'
las <- lasdetectshape(las, shp_plane(k = 25), "planar")
#> Computed in 0.1 second

all.equal(M$planar, las$planar)

# At this stage we can be clever and find that the bottleneck is
# the eigenvalue computation. Let's write a C++ version of it with
# Rcpp and RcppArmadillo
Rcpp::sourceCpp(code = "
#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]

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

plane_metrics2 = function(x,y,z, th1 = 25, th2 = 6) {
  xyz <- cbind(x,y,z)
  eigen_m <- eigen_values(xyz)
  is_planar <- eigen_m[2] > (th1*eigen_m[3]) && (th2*eigen_m[2]) > eigen_m[1]
  return(list(planar = is_planar))
}

M <- point_metrics(las, ~plane_metrics2(X,Y,Z), k = 25)
#> Computed in 0.5 seconds

all.equal(M$planar, las$planar)
# Here we can see that the optimized version is way better but is still 5 times slower
# because of the overhead of calling R functions and switching back and forth from R to C++.


# Use the filter argument to process only first returns
M1 <- point_metrics(las, ~plane_metrics2(X,Y,Z), k = 25, filter = ~ReturnNumber == 1)
dim(M1) # 13894 instead of 17182 previously.

# is a memory-optimized equivalent to:
first = lasfilterfirst(las)
M2 <- point_metrics(first, ~plane_metrics2(X,Y,Z), k = 25)
all.equal(M1, M2)
# }