Diferentes resultados en métricas calculadas con LAStools y lidR

Aug 15 2020

Tengo un modelo de altura del dosel calculado a partir de datos TLS de alta densidad en una parcela de 60 por 200 metros. Intenté calcular vóxeles con LAStools y lidR y obtuve resultados significativamente diferentes. Me preguntaba si alguien puede aclarar lo que está sucediendo. Script de Lastools que utilicé:

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 

También se calculó la inclinación y la curtosis de la altura del dosel: 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, los resultados para la asimetría son los mismos, pero los valores de curtosis son muy diferentes. ¿Alguien puede ayudarme a comprender por qué hay una diferencia tan grande?

Respuestas

3 JRR Aug 15 2020 at 19:34

Con respecto al número de vóxeles, esto puede explicarse por la alineación de los vóxeles. lidRcentra el vóxel en res/2lo que significa que la parte inferior del vóxel de tierra está en 0, no en el centro. Si LAStools tiene vóxeles centrados en 0, el cambio puede explicar la diferencia. No intenté estar seguro, pero esto tiene sentido.

Sobre la asimetría y la curtosis, su pregunta está mal planteada. No se pregunta por qué lidRproporciona una salida diferente a, LAStoolssino por qué e1071proporciona una salida diferente a LAStools. lidRhace lo mismo que en LAStoolsrealidad.

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

En lidRla curtosis se define según las fórmulas de wikipedia . Podemos adivinar que LAStoolshace lo mismo. En lidRel código está:

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

En e1071el có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)

Según la página de wikipedia, parece que corresponde al exceso de curtosis o algo así.