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.

2 Answers 2


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.

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)


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.

  • Unfortunately, this won't work for me, because after calculating the voxel statistics, I change the data in-between. I change the values depending on the neighboring voxels. Also, I actually need them to be aligned exactly. I guess I will then build rasters for each vertical voxel layer so I can use merge_spatial(). Thanks though!
    – Zoe
    Jun 4, 2021 at 8:54
  • 1
    Well you have a good basis with this code. If you have more complex requirements that what you mentioned I cannot guess it. For alignment you can modify the code to align differently. Here voxel centers starts at 0. In lidR corners start at 0. It looks like that round_any(x - 0.5 * res, res) + 0.5 * res
    – JRR
    Jun 4, 2021 at 10:09
  • 1
    Are the coordinates within the object returned by voxel_metrics() the centers of the voxels? Is this also true for the z coordinate? I voxelized a point cloud that had min(z) = 0 and max(z) = 2, but the voxels also had the same z boundaries. Like the highest and lowest voxels were cut in half? Sorry to ask here, just wondering.
    – Zoe
    Jun 4, 2021 at 13:47
  • 2
    Now you are asking I'm not sure myself. The coordinates that are returned are the center of the voxels for sure. I checked and points in Z [-0.5, 0.5[ give a voxel at Z = 0 but for X the alignment is X in [0, 1[ give a voxel at X = 0. I guess it deserves an alignment parameter
    – JRR
    Jun 4, 2021 at 13:58

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.