Add voxel attributes to LAS point cloud

Question

I would like to add data from voxels to a point cloud. Imagine I calculated some voxel statistics with voxel_metrics(), for example:

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

How can I then add the resulting attribute back to the points within their respective voxel without looping through all voxels manually?

I want to then delete points within voxels which have a specific value and keep the other ones. But I guess I will manage to do so with filter_poi() when the attributes are added.

Answer 1

What you want is merge_spatial() but for voxels objets. Sadly it does not exist in lidR. But you can easily use the fact that the point-cloud is stored in a data.table to do it yourself in one shot. It will be more efficient than using voxel_metrics() + an hypothetical merge_spatial() function.

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)

Notice that voxels won't be aligned exactly like in lidR with this code but with 10 cm voxel I guess you don't care.

Answer 2

Since my actual goal was quite specific, my code is probably a bit non-generic, but I solved my issue. I got the voxels containing points, then I altered the voxel values and then I attatched the voxel metric to the point cloud by creating rasters for each vertical layer of voxels and merging them to the point vloud with merge_spatial(). I used 10 cm pixels and did not change the bottom 5 voxel rows because I wanted them to remain unchanged. I will be honest: This code takes quite some time to execute, but it does what I want it to do. I guess it could be improved a lot.

  # 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