The lidR package has a neat way to work with huge (pointcloud) datasets: The catalog
function (doc here) avoids loading the dataset to memory and can treat mosaics [datasets that are spread across multiple (non-overlapping) tiles] as single dataset. It loads required tiles on-the-fly during computations in an intelligent way. It is great to avoid working with huge files (multiple GBs) and keep memory requirements lean if only working with small parts of the dataset.
Is there a similar convenient/memory-efficient/"lidR-catalog-way" to process large raster mosaics in R? Or more put in a more general way: Is there a way to work with mosaic raster datasets in R without merging them first?
I am aware of the mosaic
(doc) and merge
functions, which allow me to merge my tiled raster mosaic into a single raster dataset. I also found that gdal
will do so a lot faster and memory efficient. Here is a R-snippet for this:
mosaicgdal <- function(files, out) {
in_files = do.call(paste, c(as.list(files), sep = " "))
cmd = paste("gdal_merge.py -a_nodata -32767 -ps 25 25 -o", out, in_files)
system(cmd)
}
However both require materializing the whole dataset as one single file in memory (or at least on disk). Is there a way to avoid this?
My application: I am working with an enormous LAS pointcloud (several TB, everything tiled in 100 MB files) and corresponding raster datasets (around 100 GB), e.g. a high-res Terrain Model. I am usually processing small portions (delimited via spatial polygons [.shp/.kml]) or whole LAS tiles. This works beautifully memory-efficient in lidR, only loading necessary tiles:
# load several TB las tiles as catalog (only file-paths and metadata)
las_data = catalog("D:/FOLDER/WITH/11TB/OF/LAS/TILES")
# load polygon of region of interest
region_of_interest = readOGR("D:/EG/SHP/FILE/OF/ROI")
# cut out the portion of las_data overlapping the polygon
# this will load only the required tiles in memory
las_data_roi = clip_roi(las_data, extent(region_of_interest ))
# ... and create a DTM for example
dtm_roi = grid_terrain(las_data_roi , algorithm = kriging(k = 10L))
I'd like to do the same with a raster dataset:
# load raster dataset only as pointers to files (and metadata such as extent)
raster_data = raster("D:/FOLDER/WITH/LOTS/OF/RASTER/TILES")
# e.g. crop from mosaic without having to call mosaic/merge first
# avoiding having a huge single file/reading the whole dataset
raster_roi = crop(raster_data, roi)
I am usually using the raster package, but that doesn't seem to provide any such functionality.