Diferentes resultados en métricas calculadas con LAStools y lidR
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
Con respecto al número de vóxeles, esto puede explicarse por la alineación de los vóxeles. lidR
centra el vóxel en res/2
lo 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é lidR
proporciona una salida diferente a, LAStools
sino por qué e1071
proporciona una salida diferente a LAStools
. lidR
hace lo mismo que en LAStools
realidad.
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 lidR
la curtosis se define según las fórmulas de wikipedia . Podemos adivinar que LAStools
hace lo mismo. En lidR
el código está:
n * sum((z - zmean)^4)/(sum((z - zmean)^2)^2)
En e1071
el 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í.