# Interpolating LiDAR intensity images in R?

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")
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().

• You are getting what you asked for - 0.5m mean intensity. Its just that lots of the 0.5m cells have no lidar points in them. I see the same thing if I try `grid_canopy` with a small grid cell size, so I don't understand your last comment either. Could you create a 1m intensity mean raster and then subsample that? I would interpolate. It all depends on what you want to do with the data. – Spacedman Jan 23 at 20:39
• The pit-free algorithm is an algorithm designed to compute CHMs. It technically could be computed using `I` of first return instead of `Z` but lidR don't do that because it relies on nothing documented in the state of the art. The best option is maybe to look at raster interpolation of empty cells – JRR Jan 23 at 20:57
• @Spacedman, in grid_canopy you can pass an algorithm to interpolate the data and fill the empty pixels. But this function is specific for the generation o CHMs, as mentioned by JRR, and won't help me with the interpolation of intensity data. – fgg Jan 24 at 17:41
• Thanks, @JRR, I'll have a look at the interpolation options outside of lidR. Do you think it makes sense to allow interpolation in the output of grid_metrics? – fgg Jan 24 at 17:48
• grid_metrics is just the metrics over grids. If you want to do any other processing on the grid then there's things in the raster package or other image analysis packages. If you want to have a continuous model for intensity then there's interpolation algorithms for computing that given a set of (x,y,Z) points. – Spacedman Jan 24 at 18:35

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