I'm using the lidR, raster, sp, and gstat libraries in R to create a Digital Terrain Model (DTM) from LiDAR data. I'm encountering unexpected artifacts in my DTM and I'm unsure why.

Here is my code:


las <- readLAS("data.las")

ground_points <- filter_poi(las, Classification == 2)

dtm3 = rasterize_terrain(ground_points, algorithm = kriging(model = gstat::vgm(0.59, "Sph", 874), k = 10L))

plot_dtm3d(dtm3, bg = "white") 

enter image description here

I'm using kriging for the interpolation method, and I've specified a spherical model for the variogram with a nugget of 0.59, a range of 874, and k = 10.

However, when I plot the resultant DTM, I notice strange artifacts, the spikes should not be there.

Does anyone know what could be causing these artifacts? Are there any issues with my kriging parameters or with the overall methodology? How could I correct for this in order to generate a more accurate DTM?

  • Hard to tell without your data. Are points with anomalously high and low values present in the data?
    – Spacedman
    Commented Jun 6, 2023 at 9:03
  • No, these weird artifacts are mostly presented in areas where there are builidngs, Buildings only have rooftops so no point in the ground
    – Purple_Ad
    Commented Jun 6, 2023 at 9:24
  • It might be because rooftops are flat or because a building causes a discrete jump in elevation which the interpolation doesn't handle smoothly. Kriging will estimate the overall smoothness of the whole surface and then use that to interpolate, when that assumption is violated locally then maybe these things can happen... Do the other algorithms mentioned in the docs for rasterize_terrain handle it better?
    – Spacedman
    Commented Jun 6, 2023 at 9:34
  • Yes other algorithm handle the data better. Several paper mentioned that kriging is better for interpolation in city areas. I guess either they are wrong or i am making some mistakes
    – Purple_Ad
    Commented Jun 6, 2023 at 9:40
  • Where do you get your variogram parameters from? Normally you'd fit a variogram model to the data, but here (I think) the data being kriged is something like the average of 10 nearest neighbours and you don't have access to that at this point. What is 874 spatial units compared to the size of your data? How sensitive are these peaks to changes?
    – Spacedman
    Commented Jun 6, 2023 at 15:22


Your Answer

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

Browse other questions tagged or ask your own question.