grid_terrain() producing NAs and negative values

[lucas-johnson]

The grid_terrain.vrt file produced by grid_terrain() has several huge negative values (looks like it might be the most negative number possible). Here is a summary of the grid_terrain.vrt file:

vrt
class      : RasterLayer 
dimensions : 24000, 3000, 72000000  (nrow, ncol, ncell)
resolution : 1, 1  (x, y)
extent     : 535500, 538500, 4872000, 4896000  (xmin, xmax, ymin, ymax)
crs        : +proj=utm +zone=18 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
source     : memory
names      : grid_terrain 
values     : -339999999999999996123846586046231871488, 753.7401  (min, max)

It looks like there are plenty of these negative values too:

> length(vrt@data@values[vrt@data@values < 0])
[1] 21606964
> unique(vrt@data@values[vrt@data@values < 0])
[1] -339999999999999996123846586046231871488

Digging a little bit deeper I made a mosaic out of the .tif files generated by grid_terrain so that I could compare with the previous .vrt file:

dtm
class      : RasterLayer 
dimensions : 24000, 3000, 72000000  (nrow, ncol, ncell)
resolution : 1, 1  (x, y)
extent     : 535500, 538500, 4872000, 4896000  (xmin, xmax, ymin, ymax)
crs        : NA 
source     : memory
names      : layer 
values     : 470.7277, 753.7395  (min, max)

The mosaic of tifs is called dtm and at a first glance it looks better - no massive negative values. But really we just have NAs instead of the large negatives.

> length(dtm@data@values[dtm@data@values < 0])
[1] 21289879
> unique(dtm@data@values[dtm@data@values < 0])
[1] NA
> 

[Jean-Romain]

There are two things to figure out here.

• Why there are NAs in the DTM? Well, actually NAs are not a big deal as long as they belong in pixels that do not contains points
• Why you have -Inf instead of NAs? This is actually pretty obvious in theory. In practice I must investigate further to find a fix.

To answer these two questions I will need a reproducible example. Let try out of a dummy point cloud (see below). As you can see, I can reproduce but the -Inf values in the RasterLayer are actually the NAs that are not properly interpreted by the raster package. However it is consequent-less for the normalization because NAs are expected to be in pixels with no points.

In your case you seems to have NAs in pixels with point according to your gis question. So this is another issue that need a reproducible example I can't produce myself.

library(lidR)

las <- lidR:::dummy_las(500)
las@data[, Z := round(Z + 0.005*(X-30)^2 - 0.005*(Y - 50)^2 + 20,3)]
f = tempfile(fileext = ".las")
writeLAS(las, f)

ctg = readLAScatalog(f)
opt_chunk_size(ctg) = 50
opt_chunk_buffer(ctg) = 10
opt_progress(ctg) <- FALSE
plot(ctg, chunk = TRUE)

# From las to RasterLayer in memory
dtm1 = grid_terrain(las, 0.5, knnidw(5))
plot(dtm1)
dtm1
#> class      : RasterLayer 
#> dimensions : 200, 200, 40000  (nrow, ncol, ncell)
#> resolution : 0.5, 0.5  (x, y)
#> extent     : 0, 100, 0, 100  (xmin, xmax, ymin, ymax)
#> crs        : NA 
#> source     : memory
#> names      : Z 
#> values     : 8.605863, 42.92965  (min, max)

# From catalog to RasterLayer in memory
dtm2 = grid_terrain(ctg, 0.5, knnidw(5))
plot(dtm2)
dtm2
#> class      : RasterLayer 
#> dimensions : 200, 200, 40000  (nrow, ncol, ncell)
#> resolution : 0.5, 0.5  (x, y)
#> extent     : 0, 100, 0, 100  (xmin, xmax, ymin, ymax)
#> crs        : NA 
#> source     : memory
#> names      : Z 
#> values     : 8.605863, 42.92965  (min, max)

# From catalog to RasterLayers on disk
opt_output_files(ctg) <- paste0(tempfile(),"/{ID}")
dtm3 = grid_terrain(ctg, 0.5, knnidw(5))
#> Registered S3 method overwritten by 'R.oo':
#>   method        from       
#>   throw.default R.methodsS3
plot(dtm3)
dtm3
#> class      : RasterLayer 
#> dimensions : 200, 200, 40000  (nrow, ncol, ncell)
#> resolution : 0.5, 0.5  (x, y)
#> extent     : 0, 100, 0, 100  (xmin, xmax, ymin, ymax)
#> crs        : NA 
#> source     : /tmp/RtmpE5YSeL/file6f8642af2bf2/grid_terrain.vrt 
#> names      : grid_terrain 
#> values     : 8.605863, 42.92965  (min, max)

range(dtm3[])
#> [1] -3.400000e+38  4.292965e+01

lasn = lasnormalize(las, dtm3)
range(lasn$Z)
#> [1] -2.911 26.841
plot(lasn)

ctgn = lasnormalize(ctg, dtm3)
spplot(ctgn, "Min.Z")

^{Created on 2019-09-27 by the reprex package (v0.3.0)}

[Jean-Romain]

I've found why the NA are not properly interpreted. This has been fixed. As a temporary workaround you can force the NAflag in writeRaster

las <- lidR:::dummy_las(500)
las@data[, Z := round(Z + 0.005*(X-30)^2 - 0.005*(Y - 50)^2 + 20,3)]
f = tempfile(fileext = ".las")
writeLAS(las, f)

ctg = readLAScatalog(f)
opt_chunk_size(ctg) = 50
opt_chunk_buffer(ctg) = 10
opt_progress(ctg) <- FALSE
opt_output_files(ctg) <- paste0(tempfile(),"/{ID}")

# force NAflag not being -Inf because not interpreted properly by `gdalUtils`
ctg@output_options$drivers$Raster$param$NAflag = -9999999

dtm3 = grid_terrain(ctg, 0.5, knnidw(5))
plot(dtm3) # it works
range(dtm3[])

[lucas-johnson]

@Jean-Romain Thank you for looking into this, and providing a workaround.

I was able to normalize using technique outlined here: #184 (comment) This not ideal however, and is probably ignoring an underlying bug. I'll work on reproducing and submitting something for you.

[Jean-Romain]

Yes if you can make a reproducible example it will help me. I think that I know the issue (the same as in #184). I will change the method of selection of interpolation points for a better one. It actually relies on a old code that could be improved.

[Jean-Romain]

I close this issue, please open a new one if you reproduce. Here a summary of my guess on the issue:

• your point cloud is special in a way or an other and the interpolation did not covered all the points. You have a raster with NAs. This is not expected.
• when normalizing with such raster you should have an informative error but...
• because gdalUtils did not handled properly the NAs coded with -Inf your raster contains -3.8e+38 instead of NAs and the error that should have been thrown, has been bypassed leading to a dummy normalization.

It is a double bug. The first bug is handled internally an fails with an informative error but the second bug generated a side effect that enabled the code to work instead of failing.

[Jean-Romain]

I added an option full_raster in grid_terrain(). This will solve your issue (but I still need to find the bug).