Negative values after CHM rasterization - lidR

Question

I don't understand why when rasterizing normalized point clouds (using lidR, R environment) with no negative value, I can get raster canopy height model with negative values?

An example based on the sample data from lidR package:

library(lidR)
LASfile <- system.file("extdata", "Megaplot.laz", package="lidR")
las <- readLAS(LASfile)
nlas <- normalize_height(las,tin())
summary(nlas$Z)
# > summary(nlas$Z) # NO Negative values
# Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
# 0.00    7.78   14.93   13.27   19.32   29.97 

If we look at the value of the rasterized CHM, we can find negative values. The phenomenon is less clear with this dataset but with my data, these can be of several meters!

chm <- grid_canopy(nlas, res = 1, pitfree(subcircle = 0.15))
# > chm
# class      : RasterLayer 
# dimensions : 236, 228, 53808  (nrow, ncol, ncell)
# resolution : 1, 1  (x, y)
# extent     : 684766, 684994, 5017772, 5018008  (xmin, xmax, ymin, ymax)
# crs        : +proj=utm +zone=17 +datum=NAD83 +units=m +no_defs 
# source     : memory
# names      : Z 
# values     : -0.0001215559, 28.97837  (min, max)

It also occur with dsmtin() algorithm, which is really similar than the one used for height normalization.

grid_canopy(nlas, res = 1, dsmtin())
# class      : RasterLayer 
# dimensions : 235, 228, 53580  (nrow, ncol, ncell)
# resolution : 1, 1  (x, y)
# extent     : 684766, 684994, 5017773, 5018008  (xmin, xmax, ymin, ymax)
# crs        : +proj=utm +zone=17 +datum=NAD83 +units=m +no_defs 
# source     : memory
# names      : Z 
# values     : -0.0001546422, 29.11114  (min, max)

Could someone explain me these negative values ?

Answer 1

For a given position p (x,y), the interpolation of a triangulation consists in finding in which triangle ABC the location p belongs to computes its z coordinates from the coordinates of the triangle.

For a given point, there is a mathematical solution to know if the point belongs in a triangle. However in computer science because of floating point accuracy when a point is very close to the edge the test is likely to fail. This is why the computation is made with a tolerance. It is like having a buffer around the triangles. Because of this tolerance a point may be found in the wrong triangle adjacent to the actual one. This is not a big deal, it leads to 10th of millimeter inaccuracies and you probably have many cases like that in your raster. But who cares? It is far below the actual accuracy of the sensor.

But when the expected value is 0, this inaccuracy becomes visible when negative. This is the reason of the -0.0001 you spotted. In v3.0.4 (released today release next week hopefully) the tolerance has been reduced + grid_canopy() rounds the pixel elevation to do not output too many decimal digits which are not relevant. The problem is gone.

Your problem with several meters of error might be somehow similar. A bug has been reported here some weeks ago. At the very edge a Delaunay triangulation is often very poor and generates irrelevant triangles. See below (left) where some almost vertical triangles were generated. Irrelevant triangles + computational inaccuracies may lead to weird results (middle). In version 3.0.4 this has been fixed by reducing the tolerance + by checking the normal of the triangles (right).

In a CHM, I guess it might happen as well in very steep triangles if any. Try v3.0.4 to see if it fixes your issue.