Resultados diferentes em métricas calculadas com LAStools e lidR

Question

Eu tenho um modelo de altura de dossel calculado a partir de dados TLS de alta densidade em um gráfico de 60 por 200 metros. Tentei calcular voxels com LAStools e lidR e obtive resultados significativamente diferentes. Gostaria de saber se alguém pode deixar claro o que está acontecendo. Script Lastools que usei:

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

number of voxels: 189077 

código lidR:

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 

Também calculou a assimetria e curtose da altura do dossel: LAStools:

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

Resultado:

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))
})

Resultado:

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

Como podemos ver, os resultados para assimetria são os mesmos, mas os valores de curtose são muito diferentes. Alguém pode me ajudar a entender por que há uma diferença tão grande?

Answer 1

Em relação ao número de voxels, isso pode ser explicado pelo alinhamento dos voxels. lidRcentraliza o voxel, res/2significando que a parte inferior do voxel terrestre está em 0 e não no centro. Se o LAStools tiver voxels centrados em 0, a mudança pode explicar a diferença. Não tentei ter certeza, mas isso faz sentido.

Sobre assimetria e curtose, sua pergunta é mal colocada. Você não está perguntando por que lidRfornece uma saída diferente de, LAStoolsmas por que e1071fornece uma saída diferente de LAStools. lidRfaz o mesmo do que LAStoolsrealmente.

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

Na lidRcurtose é definida de acordo com as fórmulas da wikipedia . Podemos adivinhar que LAStoolsfaz o mesmo. No lidRcódigo está:

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

No e1071código está:

# 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)

De acordo com a página da wikipedia, parece que corresponde ao excesso de curtose ou algo parecido.