3

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 ?

3
  • 1
    I'd say: floating point inaccuracies for the case you are reporting but this deserve further investigation to be sure. But if you encountered several meters of error there is another issue. Please first show your raster to see where are those errors then maybe open a bug report.
    – JRR
    Oct 9, 2020 at 19:13
  • Are you referring to the 0.1mm min value? This is greater accuracy than most ALS LiDAR systems produce. I agree with @JRR in first investigating the floating point inaccuracies or interpolation variability. Could you also please confirm that your point cloud does not include negative values?
    – Aaron
    Oct 9, 2020 at 19:48
  • After investigation, for the case reported in this question I confirm it is an accuracy artifact due to some internal tolerance. I'll give a more comprehensive answer later. But an error of several meters is something else that is not related.
    – JRR
    Oct 9, 2020 at 21:03

1 Answer 1

4

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.

enter image description here

2
  • Thanks a lot for these very clear in formations! It was indeed a edge artifact on my dataset
    – Adrien
    Oct 10, 2020 at 19:42
  • The best way is to always use a buffer around your region of interest. Many methods are prone to edge effect. Using a buffer you ensure that the region of interest will be free of artifacts.
    – JRR
    Oct 10, 2020 at 21:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.