LAScatalog with user-defined grid

[antoniobelchiorferraz]

I think it would be interesting to provide a pre-defined grid cell (raster file) that would be constant independently of the LAS files. It would be useful, for instance, to resample different datasets into a common grid.

[Jean-Romain]

Do you mean providing a simple way to generate chunks that are always aligned the same independently of the collection? It is already the case. It is not easy to verify but below is a reproducible example that demonstrate the grid is common to both collections

library(lidR)
ctg1 = readLAScatalog("LiDAR/", chunk_size = 1000, progress = FALSE)

# Generate a fake LAScatalog of the same region but with misaligned tiles
ctg2 = ctg1
ctg2@data$geometry = ctg2@data$geometry + c(500, 500)
ctg2@data$Min.X = ctg2@data$Min.X + 500
ctg2@data$Max.X = ctg2@data$Max.X + 500
ctg2@data$Min.Y = ctg2@data$Min.Y + 500
ctg2@data$Max.Y = ctg2@data$Max.Y + 500

# Convert chunks into sfc to display
chunks_as_sf <- function(chunks) {
  u <- lapply(chunks, function(x) sf::st_as_sfc(st_bbox(x)))
  u <- do.call(c, u)
  u
}

# Check that the two collections are misaligned
plot(ctg1)
plot(ctg2, add = T)

# Generate the chunks
chunks1 <- engine_chunks(ctg1)
chunks2 <- engine_chunks(ctg2)

# Plot the chunks. They are aligned!
plot(chunks_as_sf(chunks1), col = rgb(0,0,1,0.3))
plot(chunks_as_sf(chunks2), col = rgb(1,0,0,0.3), add = TRUE)

`

[antoniobelchiorferraz]

Let's say I calculated metrics from RGB images according to a given raster with res=1m.

Then, I want to calculate metrics from a las point cloud using pixel_metrics, but the output raster needs to be defined in the same grid as the RGB-based metrics.

readLAScatalog("path/to/las/files/")
hmean <- pixel_metrics(ctg, ~mean(Z), res=1)

How can I ensure that both output rasters overlap? Does the output raster depend on the way I define the chunks?

[Jean-Romain]

hmean <- pixel_metrics(ctg, ~mean(Z), res=rgb)

with rgb being your raster template

Does the output raster depend on the way I define the chunks?

Absolutely not. lidR handles everything for you

[antoniobelchiorferraz]

[Jean-Romain]

Did you try? Yes it does. Check the doc.

res numeric. The resolution of the output. Can optionally be a RasterLayer or a stars or a SpatRaster. In that case the raster is used as the template

[antoniobelchiorferraz]

oh !! got it! I misunderstood your answer. I wasn't aware of that feature. Thank you!

[Jean-Romain]

You are not the first one to ask similar question. See also https://gis.stackexchange.com/questions/418332/change-extent-of-lidr-las-object-for-raster-alignment

In short: it just work the way you want it to work. You don't need to think too much about it.

[antoniobelchiorferraz]

@Jean-Romain I try it but the output raster (origin and extent) does not match the template (origin and extent).

work.dir<-".../lasFolder/"
ctg <- readLAScatalog(work.dir)
r <- raster(".../Large_ClassifiedPC_chm_20m.tif")
print(ncell(r))
#> [1] 50370
fg<-grid_metrics(ctg, ~mean(Z), res=r)
#> Processing [===========================================================================================================] 100% (4/4) eta:  0s
print(ncell(fg))
#> [1] 2500
writeRaster(fg,".../Large_ClassifiedPC_chm_20m_resampled4.tif")

[Jean-Romain]

I'm not excluding the existence of a bug but I'm more confident on the fact that it behaves as expected but not necessarily as you were expecting. The output fg is expected to be aligned with r and have the same resolution than r. The cells are expected to perfectly match. However the final extents may not match for different technical reasons. My guess here is that r is larger than ctg right? If r was smaller than ctg I think you would get the expected result.

Please show

plot(ctg)
plot(r, add = T)

[antoniobelchiorferraz]

Yes, r is much larger than ctg. I confirm that the cells of r and fg are aligned.

So, one solution might be to resample fg cells into r using a nearest neighbor approach.

[Jean-Romain]

Why do you want to resample fg? fg and r are matching. If you really want the exact same extent you can extend fg with raster or terra. It will automatically fill missing cells with NA