3

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?

2
  • Your example only enable to see that LAStools and lidR do not generate the same number of voxels. You must rework your example with skewness and kurtosis to show us what is points and make a valid and reproducible code example. points>1.3m is not a valid R syntax.
    – JRR
    Aug 15, 2020 at 10:12
  • Thanks JRR, I added R code for calculating skewness and kurtosis.
    – Sher
    Aug 15, 2020 at 12:02

1 Answer 1

3

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.

1
  • Thanks a lot JRR for detailed explanation, and sorry for confusion with e1071.
    – Sher
    Aug 15, 2020 at 13:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.