4

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

8
  • 1
    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
    Commented Jan 23, 2019 at 20:39
  • 1
    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
    Commented Jan 23, 2019 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
    Commented Jan 24, 2019 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
    Commented Jan 24, 2019 at 17:48
  • 1
    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
    Commented Jan 24, 2019 at 18:35

1 Answer 1

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

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.