1

I have several tiles in a catalog covering bays along the shore. There are some hits in the water but they are sparse. Ideally I think I should use the is_concave = TRUE as an option in grid_terrain(), but it was taking so long I ended up terminating R.

I looked at the resulting terrain rasters that had the "outside the convex hull" warning - they had a small area of extreme negative values due to the NAs in the water area. Is it possible to remove the NAs?

Data may be downloaded from NS elevation explorer (search for 284_4955).

library(lidR)

myLas <- readLAS("284_4955_201901.laz", filter = "-drop_withheld -drop_overlap -drop_z_below -1")
x <- grid_terrain(myLas, 1, tin())
#> Interpolation of 89 points outside the convex hull defined by ground points (outside the triangulation) failed and returned NAs.

x
#> class      : RasterLayer 
#> dimensions : 1000, 1000, 1e+06  (nrow, ncol, ncell)
#> resolution : 1, 1  (x, y)
#> extent     : 284000, 285000, 4955000, 4956000  (xmin, xmax, ymin, ymax)
#> crs        : +proj=utm +zone=20 +ellps=GRS80 +units=m +vunits=m +no_defs 
#> source     : memory
#> names      : Z 
#> values     : -21474836, 106.02  (min, max)

2 diagonal black lines are negative raster values which I think come from the NAs. 2 diagonal black lines are negative raster values

2
  • You problem is likely to be reproducible with a LAS object. Please simplify your question to make it minimally reproducible. Moreover I'm not sure to understand. normalize_height is likely to complain about NAs but not grid_terrain(). This is why we need a minimal reproducible example and ideally a link to the point cloud or at least an image.
    – JRR
    Commented Mar 29, 2021 at 17:31
  • Thank you, I've simplified the code and included an image in the edits of my question.
    – Ray J
    Commented Mar 29, 2021 at 20:16

2 Answers 2

1

Well, those cases of fractal-ish costal region with sparse points in water are nasty. And your files is special: points classified as water are also flagged withheld. So using drop_withheld discards all water points. In addition I think you found an edge case bug (when you see -21474836 you are in presence of an infinite value). So here are some options available for you.

Keep water points

myLas <- readLAS("284_4955_201901.laz", filter = "-drop_z_below -1 -keep_class 2 9")
x <- grid_terrain(myLas, 1, tin())
plot(x)

enter image description here

Increase extrapolation outide convex hull limits

See documentation of tin() and knnidw()

myLas <- readLAS("284_4955_201901.laz", filter = "-drop_withheld -drop_overlap -drop_z_below -1 -keep_class 2 9")
x <- grid_terrain(myLas, 1, tin(extrapolate = knnidw(3, 1, 75)))
plot(x)

enter image description here

Erase infinites

Erase infinite values that actually are in water at the very edge of the convex hull.

myLas <- readLAS("284_4955_201901.laz", filter = "-drop_withheld -drop_overlap -drop_z_below -1 -keep_class 2 9")
x <- grid_terrain(myLas, 1, tin())
x[x < 0] <- 0
plot(x)

enter image description here

2
  • perfect, thank you. I think in my case extending the extrapolation with knnidw is ideal and safer than my initial request. This should reduce/eliminate the infinite value but can be double checked aftterwards.
    – Ray J
    Commented Mar 30, 2021 at 11:01
  • 1
    To be correct I proposed to extend the limits of the extrapolation with knnidw(). This task is always performed but knnidw does not interpolate if there is no point less that 50 m away. Here I choose 75 m. This is a typical problem already encountered by other users in costal region
    – JRR
    Commented Mar 30, 2021 at 11:04
0

A slight modification to JRR's answer of removing infinites is to set NA values to NaN. In this way ArcGIS will read these as noData rather than a number. (I chose < -1 because it's below what I used in the filter for reading the las.)

myLas <- readLAS("284_4955_201901.laz", filter = "-drop_withheld -drop_overlap -drop_z_below -1")
x <- grid_terrain(myLas, 1, tin())

x[x < -1] <- NaN
NAvalue(x) <- NaN  # not sure if needed.
    
writeRaster(x, "dtm.tif", naflag = -9999)

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.