1 votos

Añadir atributos de voxel a la nube de puntos LAS

Me gustaría añadir datos de voxels a una nube de puntos. Imagina que he calculado algunas estadísticas de los vóxeles con voxel_metrics() por ejemplo:

voxels <- voxel_metrics(las, length(X), res=0.1, all_voxels=TRUE)

¿Cómo puedo entonces añadir el atributo resultante a los puntos dentro de su respectivo vóxel sin recorrer todos los vóxeles manualmente?

Quiero entonces eliminar los puntos dentro de los vóxeles que tienen un valor específico y mantener los otros. Pero supongo que lo conseguiré con filter_poi() cuando se añaden los atributos.

2voto

Andrey Atapin Puntos 384

Lo que quieres es merge_spatial() pero para los objetos de voxels. Lamentablemente no existe en lidR . Pero se puede utilizar fácilmente el hecho de que la nube de puntos se almacena en un data.table para hacerlo tú mismo de una sola vez. Será más eficiente que usar voxel_metrics() + un hipotético merge_spatial() función.

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

res  = 8
xvox = plyr::round_any(las$X, res)
yvox = plyr::round_any(las$Y, res)
zvox = plyr::round_any(las$Z, res)

las@data[, N := length(Z), by = list(xvox, yvox, zvox)]

las2 = filter_poi(las, N > 50)

plot(las2)

Fíjate que los voxels no estarán alineados exactamente como en lidR con este código pero con voxel de 10 cm supongo que no te importa.

0voto

Sterling Duchess Puntos 176

Como mi objetivo real era bastante específico, mi código es probablemente un poco no genérico, pero he resuelto mi problema. Obtuve los vóxeles que contenían puntos, luego alteré los valores de los vóxeles y después adjunté la métrica de los vóxeles a la nube de puntos creando rásters para cada capa vertical de vóxeles y fusionándolos con la nube de puntos con merge_spatial() . Utilicé píxeles de 10 cm y no cambié las 5 filas de vóxeles inferiores porque quería que no se modificaran. Voy a ser honesto: este código tarda bastante tiempo en ejecutarse, pero hace lo que quiero que haga. Supongo que se podría mejorar mucho.

  # voxel metrics: is there a point in the voxel or not?
  voxels <- voxel_metrics(las, length(X), res=0.1, all_voxels=TRUE)
  voxels$V1 <- ifelse(voxels$V1 > 0, 1, 0)
  voxels$V1[is.na(voxels$V1)] <- 0

  # here I filter / change my voxels

  # add empty voxel attributes to the points
  las <- add_lasattribute(las, 1, "V1", "keep voxels with 1")

  # save points which should remain the same
  z_loop_vals <- sort(unique(voxels$Z))[6:length(unique(voxels$Z))]
  unchanged_las <- filter_poi(las, Z <= min(z_loop_vals-5)/100)

  # for each (filtered) vertical voxel layer, create a raster
  for (z_val in z_loop_vals) {

    # create raster
    z_subset <- voxels[voxels$Z==z_val,]
    z_subset <- as.data.frame(z_subset)[,c(1,2,4)]
    new_raster <- rasterFromXYZ(z_subset)
    crs(new_raster) <- CRS("+init=EPSG:25832")

    # add raster values to point cloud
    las_z <- filter_poi(las, Z > (z_val-0.05) & Z <= (z_val+0.05))
    las_z <- merge_spatial(las_z, new_raster, "V1")

    # remove points with V1 == 0
    las_z <- filter_poi(las_z, V1 == 1)
    unchanged_las <- rbind(unchanged_las, las_z)
  } 

  # overwrite old point cloud
  las <- unchanged_las

i-Ciencias.com

I-Ciencias es una comunidad de estudiantes y amantes de la ciencia en la que puedes resolver tus problemas y dudas.
Puedes consultar las preguntas de otros usuarios, hacer tus propias preguntas o resolver las de los demás.

Powered by:

X