3

I want to use watershed with lidR package on a catalog. I paste this command:

rep <- "/home/lasfiles"
new_ctg <- catalog(rep)
lidR:::catalog_laxindex(new_ctg)
opt_output_files(new_ctg) <- paste0(rep, "/outputs/mnh/{ID}_cano")
opt_chunk_size(new_ctg) <- 500
opt_chunk_buffer(new_ctg) <- 30
opt_cores(new_ctg) <- 7L
algo <- pitfree(thresholds = c(0,10,20,30,40,50), subcircle = 0.2)
chm <- grid_canopy(new_ctg, 0.5, algo)

The result CHM is ok,

summary(chm)

        grid_canopy
 Min.       1305.953
 1st Qu.    1406.373
 Median     1468.121
 3rd Qu.    1539.043
 Max.       1716.370
 NA's        180.000
 Warning message:
 In .local(object, ...) :
     summary is an estimate based on a sample of 1e+05 cells (5% of 
     all cells)

but when I pass this command:

algo <- watershed(chm, th = 4)

Error in assert_is_all_of(chm, "RasterLayer") : chm is not a RasterLayer it is RasterStack.

Why this error?

  • 1
    Seems like a bug. It should be a RasterLayer. Please make a reproducible example. Your code does not enable us to reproduce. This is regularly tested internally in the package and grid_canopy returns always a RasterLayer. You either made a mistake somewhere else or you found a limit case with a bug. We need more context. – JRR Feb 14 at 12:28
  • 1
    Edit your question and if you can't demonstrate this on data we can use you could at least show the output of summary(chm) for context. – Spacedman Feb 14 at 14:38
1

You chose to write the rasters into files with the option opt_output_files. By choosing this option, the output of grid_canopy is no longer a RasterLayer but a Virtual Raster Mosaic (see gdalUtils::gdalbuildvrt).

In lidR VRT are loaded back as RasterStack and not a RasterLayer as expected. This is why you encountered such error. This is an improvement that should be made.

You can extract only the first (and actually unique) layer as RasterLayer:

chm = chm[[1]]

That being said, lidR is made to process point-clouds and potentially wide point-clouds as far as possible but it is not made for processing wide rasters. Thus the watershed function has not been designed to efficiently process your rasters. I'm not even sure it can process a VRT. When created, it was more designed to work with small size CHM. But it should work anyway with a RasterLayer. Just be careful it will be all loaded in memory.

What you found may be considered as a bug. I opened an issue to report an enhancement.

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.