2

I am trying to generate a normalized height values for a catalog of las tiles. I first generated a dtm using:

dtm = grid_terrain(ctg, knnidw(k=5))

Then, I normalize the values using:

normalized = lasnormalize(ctg, dtm)

The normalized data set shows some negative values, however they happen to be wildly large. I'm wondering what I might be doing that could cause this to happen. Here is some more info on the data I'm using:

> head(data.table(ctg@data)[, .(Min.Z)][order(Min.Z)])
    Min.Z
1: 396.91
2: 430.77
3: 456.20
4: 456.51
5: 458.07
6: 463.13
> dtm@data@max
[1] 753.7401
> head(data.table(normalized@data)[, .(Min.Z)][order(Min.Z)])
          Min.Z
1: -21474836.48
2: -21474836.48
3:      -183.38
4:      -171.77
5:       -19.05
6:       -17.92

Given the largest value in the dtm, I don't see how -21474836.48 could possibly be generated.

ctg is my catalog of las files.

dtm is the raster generated from grid_terrain.

normalized is the result of lasnormalize(ctg, dtm)


It turns out that none of the individual chunked/tiled .tif rasters have negative values. 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 grid_terrain:

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 too:

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

So a better set of questions:

Should I be using grid_terrain.vrt as my dtm for normalization?

If so, why could there be so many of these huge negatives in grid_terrain.vrt?

If not, how should I be using the output of grid_terrain() to normalize my las files?

  • Are you saying that each raster file looks correct but that the virtual mosaic is invalid? It looks like a bug that should be reported. – JRR Sep 27 '19 at 1:36
  • 1
    Issue submitted: github.com/Jean-Romain/lidR/issues/283 I have provided some further details there, but in summary the virtual mosaic has large negative values where as the .tifs have NAs . Both cases seem problematic. – Lucas Sep 27 '19 at 4:05
3

There are several questions here:

  1. Why do you have NAs in the DTM? NAs in the DTM are usually not a big deal. lidR interpolates within the convex hull of the point cloud to ensure to have a DTM in accordance with the point cloud especially with circular plots for example. A raster being rectangular you can have NAs in pixels with no points

  2. Why do you have -3.4e+38 in the GDAL virtual raster? NAs is not an existing value for a binary file. Usually NAs are replaced by an arbitrary value (typically -Inf) which is interpreted as NA at read time. Here the raster package did not interprets the NAs properly. Why? I don't know. Maybe a bug in raster, maybe in gdalUtils maybe lidR did not use those packages the best way.
    Edit: this one has been fixed in lidR v2.1.4 it was an issue with gdalUtils::gdalbuildvrt that badly handle the -Inf case.

  3. Is it a big deal? In theory no, we don't care. Those NA pixels do not contain any point. However in your case it seems that some NA pixels do contain points. It should have fail with an informative error but the -3.4e+38 instead of NAs make the code bypassing the conditional statement. It is a double bug. The second bug created a side effect and the first one has not been caught. That should be reported with a reproducible example.

Below a reproducible example but with no normalization issue:

library(lidR)

# Use an internal function to generate random values
las <- lidR:::dummy_las(500)
las@data[, Z := round(Z + 0.005*(X-30)^2 - 0.005*(Y - 50)^2 + 20,3)]

# Write it in a las file
f = tempfile(fileext = ".las")
writeLAS(las, f)

# Processing settings
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))
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)

# We have -3.4e+38 instead of NAs (fixed in v2.1.4)
range(dtm3[])
#> [1] -3.400000e+38  4.292965e+01

# But we can normalize the LAS object
lasn = lasnormalize(las, dtm3)
range(lasn$Z)
#> [1] -2.911 26.841
plot(lasn)

# And we can normalize the LAScatalog object
ctgn = lasnormalize(ctg, dtm3)
spplot(ctgn, "Min.Z")
| improve this answer | |

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.