3

I'm working with the R package lidR to process point clouds. I have a raster with 4 coarse land cover classes and the corresponding tile with the LiDAR data in las format (both in the same coordinate system). Now I want to attach the land cover classes to the point cloud with the lasmergespatial tool from the lidR package.

#import las
mylas<-readLAS(files="xyz.las", select = "xyzia", filter = "")


#import raster as rasster (not as stack)
landcover<-raster("LC_raster.tif")

#classify las with raster
class_las<-lasmergespatial(mylas, landcover, attribute = "class")

So basically I don't get errors but the attribute class in the new created las file is empty with NA. What might be a problem?

  • You say they are in the same coordinate system but you should show us. What is the extent and size and projection of those R objects? If you can't share the data, at least show the output when printing the objects so we can see for ourselves. – Spacedman Oct 1 at 8:28
  • I prepared the few code lines and the two datasets in this folder. feel free to download it. Thanks a lot for any helpful comments. polybox.ethz.ch/index.php/s/9pPMM3i6691BW83 – Reto Oct 1 at 9:53
2

I'd say this was a bug in lasmergespatial that occurs when the raster is not stored in memory.

The code eventually ends up here:

lidR:::lasmergeRasterLayer = function (las, raster) 
{
    cells <- raster::cellFromXY(raster, coordinates(las))
    return(raster@data@values[cells])
}

but extracting the cell values like that only works for in-memory rasters. For your large raster:

> landcover@data@values
logical(0)
> inMemory(landcover)
[1] FALSE

By default the function raster() does not load in memory. You can force the data to be loaded in memory with readAll() if your raster is not too big:

> landcover = readAll(landcover)

readAll() will pull the values into memory and it will work. But it might also kill your pc...

You could also try this. First get the cell indices of each of the LAS points by doing what the function does:

cells <- raster::cellFromXY(landcover, lidR:::coordinates(mylas))

that gets you a vector of about 29,000,000, one value for each LAS point, and then:

v = landcover[cells]

will (eventually) pull out the 29,000,000 values of the raster at each of the raster locations. But it may be slow (1000000 took a few seconds).... I've given up waiting after a minute. Test with a small raster first.

Have reported this to the lidR package: https://github.com/Jean-Romain/lidR/issues/285

  • Nice catch. It has been fixed in v2.1.4. You can indeed force the data loading with readAll. raster("filename") does not read in memory by default even for small raster (at least on linux) – JRR Oct 1 at 12:09
  • Thanks! It worked with windows as well. – Reto Oct 1 at 12:45

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.