3

I am using the lidR package in R to manage my ALS data. I was having problems generating DTMs on the edge files of my catalog because the tool lidR::grid_terrain() was doing some interpolation between points that define the boundary of the point cloud. I attach an image presenting the error: plotting of edge file DTM and Slope(degrees) of edge file

I came up with what I think is just a small fix for my case but maybe is helpfull for a small analysis:

las <- readLAS(file)
#projection(las) <- CRS("+init=epsg:REF") #In case it lacks projection

dtm <- grid_terrain(las, res = 1, tin())

p <- raster::rasterToPolygons(lidR::grid_density(las)*0+1, dissolve = TRUE) #To create a mask
m <- raster::mask(raster::crop(dtm, raster::extent(p)),p) #To mask the DTM

The result is a DTM perfectly masked to where there are points. Result of masking

I was wondering if there is another method for big areas as I am afraid of this could take a long time to process al the files in a catalog.

1

1 Answer 1

4

There is no simple and fast solution. Your DTM is indeed incorrect, but only where there are no points. Consequently this problem that is a known issue was considered not big problem.

grid_terrain() has a parameter is_concave = FALSE. If your turns it to TRUE it computes a concave hull and interpolates only in the hull. Sadly it is slow to compute and will increase the computation time. You have only ~30 tiles so it can do the job if the point density is not super high. Try it and see if it can works for you.

Your solution is perfectly fine and probably faster. Yet I would try with a larger resolution to avoid the small "holes" I can see in your image. Or maybe fill the holes before to mask.

The lidR package has a function catalog_boundaries() that computes the concave hull of each tile exactly like grid_terrain() but that is much faster. You could use the polygons to mask the pixels. You can compute the hull on a thinned point-clouds to go faster.

5
  • I did not know about this package. Thanks for the head up. I have tried increasing the resolution of ´´´grid_density´´´ and got better results with less holes. I am trying to use the ´´´is_concave = TRUE´´´ parameter but it indicates that this argument is unused and returns Error. Feb 19, 2021 at 11:55
  • 1
    ìs_concave was introduced in lidR 3.1.0
    – JRR
    Feb 19, 2021 at 12:07
  • Just a note that catalog_boundaries() is slow for me even with concavity = Inf, as in a 3+ hour ETA from lidR 4.0.1 on <15 GB of .laz tiles with a Zen 3 core at 4.8 GHz. Since lasboundary.exe -convex (from LAStools) completes the same task in less than two minutes it may be more efficient to read its shapefiles from R.
    – Todd West
    Sep 28, 2022 at 15:15
  • When I have done ctg2 <- catalog_boundaries(ctg, concavity=1, length_threshold=15), and display it with plot(ctg2, add=TRUE), how can I convert the polygon to for example a shapefile
    – BERA
    May 17, 2023 at 8:26
  • @BERA you can convert a 'catalog' to a polygon layer by using 'sf::st_as_sf(ctg2)' May 26, 2023 at 7:07

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.