Interpolating LiDAR intensity images in R?

Question

I'd like to generate 0.5 m intensity images using lidR, but the resulting raster is full of pits. Here's a reproducible example:

library(lidR)
LASfile <- system.file("extdata", "MixedConifer.laz", package="lidR")
las <- readLAS(LASfile)
int <- grid_metrics(las, list(imean = mean(Intensity)), 0.5)
plot(int, col = height.colors(50))
summary(values(int))

Can the pit-free algorithm, or some other interpolation solution in lidR, be used to improve the result? It works great in the generation of CHMs with grid_canopy().

Answer 1

As commented, the pits/holes in the final intensity image is due to missing data.

Using OP's example, a resolution of 0.5 units yields 9244 NA cells, while doubling it (squaring in area) yields only 28 NA cells. When gridding LiDAR data, it is a best practice starting the cell size with at least the average point spacing from the point cloud (which in this case is indeed ~ 0.46 units; see from the point cloud header), but be prepared to find an optimum cell resolution around 2 times or 3 times the average point spacing. See for example, LiDAR (.las) to bare-earth DEM: average points spacing and cell size?.

That being said, if interpolation is indeed required, try to investigate which interpolation technique is more suitable for data/statistic being investigated. For example, Ashraf et al. (2017) - "An Investigation of Interpolation Techniques to Generate 2D Intensity Image from LiDAR data" concluded that Inverse Distance Weighting (IDW) and Nearest Neighbour (NN) were the most suitable interpolation methods for their set of conditions.

As explained by user JRR, the pit-free algorithm used for CHMs in lidR package is not available for interpolation of statistic 'intensity' because it has not been proved/tested/supported to be a suitable method to this type of analysis/processing task yet.

Meanwhile, a solution posted by OP in comments was to use a focal filter using a 3 x 3 window (equally weighted) to average values, and applying the focal statistic only to NA cells:

while (any(is.na(values(int)))) { int <- focal(int, w = matrix(1, 3, 3), fun = function(x){mean(x, na.rm = TRUE)}, pad = TRUE, NAonly = TRUE) }

---

_References:

_{Ashraf, I., Hur, S., & Park, Y. (2017). An Investigation of Interpolation Techniques to Generate 2D Intensity Image from LIDAR Data. IEEE Access, 5(June), 8251–8260. doi:10.1109/ACCESS.2017.2699686}