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