Different results in metrics calculated with LAStools and lidR

Question

I have canopy height model calculated from high density TLS data on 60 by 200 meter plot. I tried to calculate voxels with LAStools and lidR and got significantly different results. I was wondering if someone can make clear what is happening. Lastools script that I used:

lasvoxel -i infile.laz -drop_class 2 -step 0.5 -o outfile.las

number of voxels: 189077 

lidR code:

las = readLAS("infile.laz", select = "xyzc", filter = "-drop_class 2 -drop_z_below 0 ")
voxels <- voxelize_points(las, res = 0.5)

number of voxels: 196257 

Also calculated canopy height skewness and kurtosis: LAStools:

lascanopy -i infiles\*.laz -kur -ske -height_cutoff 1.3 -files_are_plots -names -o outfile.csv

Result:

plots        ske       kur
72a-4.laz   1.0905  5.58125
11a-4.laz   0.362   2.594
34-2.laz    0.1675  2.00875
63a-1.laz   -0.3115 2.36

lidR:

library(e1071)

files <- list.files(path= "/files", pattern= "*.laz", full.names = TRUE, recursive = FALSE)
O = lapply(files, function(x) {

  las <- readLAS(x, select = "xyzc") 

  z <- las$Z
  z_canopy <- z[z>=1.3]
  skew <- skewness(z_canopy)
  kur <- kurtosis(z_canopy)
  
  return(data.frame(file=x, skewH = skew, kurH=kur))
})

Result:

plots          ske        kur
72a-4.laz   1.090595768  2.58132381
11a-4.laz   0.362007296  -0.40599745
34-2.laz    0.167542141  -0.991227478
63a-1.laz   -0.311523396 -0.640029907

As we can see, the results for skewness are the same, but kurtosis values are very different. Can someone please help me understand why there is such a big difference?

Answer 1

Regarding the number of voxel this may be explained by the alignment of the voxels. lidR centers the voxel at res/2 meaning that the bottom of the ground voxel is at 0 not the center. If LAStools has voxels centered on 0 the shift may explain the difference. I did not try to be sure but this makes sense.

About skewness and kurtosis your question is ill-posed. You are not asking why lidR provides a different output than LAStools but why e1071 provides a different output than LAStools. lidR does the same than LAStools actually.

library(e1071)
library(lidR)

LASfile <- system.file("extdata", "Megaplot.laz", package="lidR")
las <- readLAS(LASfile, filter = "-drop_z_below 1.3")

# e1071
skewness(las$Z) # -0.42
kurtosis(las$Z) # -0.66

# lidR
cloud_metrics(las, .stdmetrics_z)[c("zskew", "zkurt")]
# -0.42
# 2.33

# LAStools
system("lascanopy.exe -i Megaplot.laz -kur -ske -height_cutoff 1.3 -files_are_plots -names -o output.csv")
# -0.42
# 2.33

In lidR the kurtosis is defined according to wikipedia's formulas. We can guess that LAStools does the same. In lidR the code is:

n * sum((z - zmean)^4)/(sum((z - zmean)^2)^2)

In e1071 the code is:

# Here the same formula
r <- n * sum((x-xmean)^4)/(sum((x-xmean)^2)^2)

# Then why output is different
y <- if (type == 1) 
  r - 3
else if (type == 2)
  ((n + 1) * (r - 3) + 6) * (n - 1)/((n - 2) * (n - 3))
else 
  r * (1 - 1/n)^2 - 3

return(y)

According to wikipedia page, it seems it corresponds to the excess kurtosis or something like that.