lasmergespatial fails for non-in-memory rasters

[barryrowlingson]

If you try and lasmergespatial with a raster that isn't in-memory, it fails silently returning NAs.

This line:

https://github.com/Jean-Romain/lidR/blob/c02abc6345138b10e6cd223c7dccae56e2fbc986/R/lasmergespatial.r#L218

is the culprit because raster@data@values[cells] does not contain the data for non-in-memory rasters:

> landcover@data@values[cells[1:10]]
 [1] NA NA NA NA NA NA NA NA NA NA

the entirety of that element being:

> landcover@data@values
logical(0)

if inMemory(landcover) is FALSE.

I think you can extract the data values by indexing directly by cell number:

> landcover[cells[1:10]]
 [1] 0 0 0 0 0 0 0 0 0 0

(yes this raster is all zeroes at those locations)

and I don't think that gives any efficiency loss to in-memory rasters. It might still be slow for large, non-in-memory rasters but that's the story with large rasters. You could test using raster::inMemory and print a warning...

SO Question that spotted this with (semi-) reproducible example data:

https://gis.stackexchange.com/questions/337359/classify-point-cloud-from-land-cover-raster-values-with-lidr-in-r/337388#337388

[Jean-Romain]

Thank you for reporting. I made a reproducible example and I fixed it as you suggested. This line of code was a very old one. I would not have made such error by now.

library(lidR)
LASfile <- system.file("extdata", "Megaplot.laz", package="lidR")
las   <- readLAS(LASfile)
shp <- system.file("extdata", "lake_polygons_UTM17.shp", package = "lidR")
lakes <- shapefile(shp)
ref = lidR:::rOverlay(las, 1)
f = tempfile(fileext = ".tif")

rlakes.ram = rasterize(lakes, ref)
inMemory(rlakes.ram)
#> [1] TRUE

writeRaster(rlakes.ram, f)

rlakes.disk = raster(f)
inMemory(rlakes.disk)
#> [1] FALSE

las = lasmergespatial(las, rlakes.ram, "rlake.ram")
las = lasmergespatial(las, rlakes.disk, "rlake.disk")

all(is.na(las$rlake.ram))
#> [1] FALSE
all(is.na(las$rlake.disk))
#> [1] TRUE
all.equal(las$rlake.ram, las$rlake.disk)
#> [1] FALSE

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

[spono]

Ciao JR,
using 2.1.5, I'm having the same issue while colorising a point cloud with a RasterStack, getting the following warning (plus NAs in data):

Warning messages:
1: In max(R, na.rm = TRUE) :
  no non-missing arguments to max; returning -Inf
2: In max(G, na.rm = TRUE) :
  no non-missing arguments to max; returning -Inf
3: In max(B, na.rm = TRUE) :
  no non-missing arguments to max; returning -Inf

Everything works when I force the load with readAll, as proposed in the mentioned post.

Just to understand: the desired behaviour from the your persective is "load in memory and do it"? because in that case a specific warning or note in the wiky could help the user in understanding why he gets NAs.

[Jean-Romain]

It is expected to work, please open an issue with a reproducible example.