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

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 '19 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 '19 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 '19 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 '19 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 '19 at 18:35
@Spacedman summarised correctly the situation. grid_metrics made its job then you have a raster and to processed raster you have the raster package. grid_canopy is different it does no aim to compute metrics but a CHM as nice as possible. — JRR, Jan 24 '19 at 18:45
Got it. Thanks again. Here's a quick fix using the raster package: while (any(is.na(values(int)))) { int <- focal(int, w = matrix(1, 3, 3), fun = function(x){mean(x, na.rm = T)}, pad = TRUE, NAonly = T) } — fgg, Jan 24 '19 at 19:30
Although this works fine for the example in the question, masking should be used if the NAs are outside of the area actually covered by the point cloud data. Also, this is probably not a good solution for larger datasets. — fgg, Jan 24 '19 at 19:39

Andre Silva · Answer 1 · 2019-02-11T16:06:41.403

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}

Interpolating LiDAR intensity images in R?

1 Answers1