2

I would like to normalize a set of LAS files produced by Digital Aerial Photogrammetry (DAP) against a set of LAS files acquired by Lidar. The idea is to create the DEM raster from Lidar once and to reuse it when processing chunks for the DAP points.

The code creates two catalogues: one for the DAP, and one for the Lidar. Before entering into the function with catalog_apply(), I would like to create a DEM from the Lidar using the extent of the DAP. My code does reduce the LAS files referenced by the LAScatalog, but when applying the grid_terrain() function the now empty tiles are still processed but result in empty tiles (i.e. the catalogue recognizes the input option of the filter but the catalogue itself remains the original size). If I could reindex the catalogue it would only process LAS files that I need - is this possible? Or should I run grid_terrain() within the catalog_apply() function separately for each chunk.

library(lidR)
    
DAP_LASCatalogue <- readLAScatalog("haspeaks/LASDAP/", select = "xyzc")
        
# eventually I would like to input extents directly from the DAP Catalogue
DAPExtent <- extent(DAP_LASCatalogue)
min_x <- xmin(DAPExtent)
max_x <- xmax(DAPExtent)
min_y <- ymin(DAPExtent)
max_y <- ymax(DAPExtent)
print(c(min_x, min_y, max_x, max_y))
    
LIDAR_LASCatalogue <- readLAScatalog("J:/2018/BlockE_2017/LAS",
                                      select = "xyzc",
                                      filter = "-keep_xy 539052.1 5050814.6 556597.7 5059891.5",
                                      progress = TRUE)
plot(DAP_LASCatalogue)
plot(LIDAR_LASCatalogue)
    
LIDAR_DTMGrid <- grid_terrain(LIDAR_LASCatalogue, 1, knnidw(), use_class = c(2L, 9L), Wdegenerated = TRUE)

1 Answer 1

1

What I understand is that you want to select some tiles from the ALS point-cloud that encompass the DAP point cloud. I think catalog_intersect() can do the job

ctg1 = readLAScatalog("LiDAR/")
ctg2 = readLAScatalog("DAP/")
ctg3 = catalog_intersect(ctg1, ctg2)

plot(ctg1)
plot(ctg2, add = T)
plot(ctg3, add = T, col = rgb(1,0,0,0.3))

Created on 2020-12-15 by the reprex package (v0.3.0)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.