grid_terrain with tin - triangulation outside of point cloud "bounds"

[lyonsdata]

Hello Jean-Romain,

I have recently been learning to use the lidR package and have had a great time with it so far. However, I am having a rather specific issue when producing a dtm using the tin algorithm from one or more .las files. Triangles are being formed that connect points on the outer bounds of the point clouds thus creating large triangles over "empty" space. This is likely due to the irregular way the point cloud has been cut. Is there any way to restrict triangles from being formed outside of the point cloud bounds? I suspect this can be fixed somehow with clipping or buffering, but I am just not sure how to go about it.

Thank you,
Mike

ctg <- readLAScatalog("/path")
dtm <- grid_terrain(ctg, res = 0.5, algorithm = tin(extrapolate = knnidw(1, 1, 50))
plot_dtm3d(dtm, bg = "white")

[woolleyt]

Hi Mike - I have experienced the same issue creating a dtm and chm using various lidR methods. This only occurs when, as you say, the point cloud is an irregular shape. I haven't figured out any way to deal with it other than clipping the output raster using a shape file of the original extent of the data. Thanks, Travis

[Jean-Romain]

There is no native way in lidR to deal with that yet in tin(). But I can easily enable it because the code do exist it is simply disabled in tin() ( but enabled in dsmtin()). The reason for that is because if I allows for removing large triangles (because this is the way I deal with this kind of shape) it has other consequences such as removing also the DTM in lakes when points are missing. Also with such feature it is easy to generate a DTM that do not encompass all the points. This is why I chose to do not provide the feature.

[lyonsdata]

I see. This would be an excellent feature if you would be willing to implement it. I understand the issue with lakes and regions with sparse points, but it would be nice to at least have the option to try. Perhaps there is a more specific implementation that would only apply this to large triangles on the outer bounds of a point cloud. In either case, thanks for the great package!

[Jean-Romain]

Perhaps there is a more specific implementation that would only apply this to large triangles on the outer bounds of a point cloud.

This is definitively complex to implement. I already thought about it. I already have some code that might be useful for such feature. But it is not that simple and again prone to unexpected troubles.

I'm currently in vacation. I'll look at that later. If you want to have a look you can probably enable it yourself by passing a value to trim here

https://github.com/Jean-Romain/lidR/blob/37b81e713270c0697f5b42c439518502b4b06693/R/algorithm-spi.R#L63

You must also pay attention to the behavior of the extrapolation from L65 that will fill the missing values and invalidate the removal of large triangle. This is the most complex part to modify actually.

[Jean-Romain]

Also please share your files. That will help me to work with a good example of datasets that presents the issue you are experiencing

[lyonsdata]

Okay, I will try this out. I also sent you a link to the .las files.

[Jean-Romain]

I added and argument is_concave in grid_terrain(). If TRUE it computes a concave hull instead of a convex hull to dertermine where are the limits of the dataset. However it is much computationally demanding. Your dataset being huge (25 pt/m², 4km²) I only made a test on a 10^th of the ground points.

shade = function(dtm) {
  na.proj <- is.na(projection(dtm))
  dummy.proj <- CRS("+init=epsg:26917")
  if (na.proj) projection(dtm) <- dummy.proj
  dtm_prod <- terrain(dtm, opt = c("slope", "aspect"))
  dtm_hillshade <- hillShade(slope = dtm_prod$slope, aspect = dtm_prod$aspect)
  if (na.proj) projection(dtm) <- CRS()
  return(dtm_hillshade)
}

library(lidR)
las = readLAS("pt000001.las", select = 'c', filter = "-keep_class 2 -keep_random_fraction 0.1")

dtm1 = grid_terrain(las, 1, tin())
sdtm1 = shade(dtm1)

dtm2 = grid_terrain(las, 1, tin(), is_concave = TRUE)
sdtm2 = shade(dtm2)

sdtm = addLayer(sdtm1, sdtm2)
plot(sdtm, col = gray.colors(50,0,1))

[lyonsdata]

Great, thank you!