8

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.

1 Answer 1

10

You should look into the terra package which provides exactly the functionality you're looking for through virtual raster tiles (VRTs). We can use them to treat a collection of raster files on disk as a single raster file while taking advantage of the API to do a majority of the same tasks as you can do through the raster package.

First, let's create a sample of 4 rasters using the example straight from the ?terra::vrt() documentation.

library(terra)

r <- rast(ncols=100, nrows=100)
values(r) <- 1:ncell(r)
x <- rast(ncols=2, nrows=2)
filename <- paste0(tempfile(), "_.tif")
ff <- makeTiles(r, x, filename)
ff
#> [1] "/var/folders/b7/_6hwb39d43l71kpy59b_clhr0000gn/T//RtmpACJYNv/filedf6b65d4fca4_1.tif"
#> [2] "/var/folders/b7/_6hwb39d43l71kpy59b_clhr0000gn/T//RtmpACJYNv/filedf6b65d4fca4_2.tif"
#> [3] "/var/folders/b7/_6hwb39d43l71kpy59b_clhr0000gn/T//RtmpACJYNv/filedf6b65d4fca4_3.tif"
#> [4] "/var/folders/b7/_6hwb39d43l71kpy59b_clhr0000gn/T//RtmpACJYNv/filedf6b65d4fca4_4.tif"

Now, we'll read them in as a VRT, again, straight from the same example. This allows

vrtfile <- paste0(tempfile(), ".vrt")
v <- vrt(ff, vrtfile)
head(readLines(vrtfile))
#> [1] "<VRTDataset rasterXSize=\"100\" rasterYSize=\"100\">"                                                                                                                                                                                                                                                                                                             
#> [2] "  <SRS dataAxisToSRSAxisMapping=\"2,1\">GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]</SRS>"
#> [3] "  <GeoTransform> -1.8000000000000000e+02,  3.6000000000000001e+00,  0.0000000000000000e+00,  9.0000000000000000e+01,  0.0000000000000000e+00, -1.8000000000000000e+00</GeoTransform>"                                                                                                                                                                             
#> [4] "  <VRTRasterBand dataType=\"Float32\" band=\"1\">"                                                                                                                                                                                                                                                                                                                
#> [5] "    <NoDataValue>nan</NoDataValue>"                                                                                                                                                                                                                                                                                                                               
#> [6] "    <ColorInterp>Gray</ColorInterp>"
v
#> class       : SpatRaster 
#> dimensions  : 100, 100, 1  (nrow, ncol, nlyr)
#> resolution  : 3.6, 1.8  (x, y)
#> extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#> source      : filedf6b216a737.vrt 
#> name        : filedf6b216a737 
#> min value   :               1 
#> max value   :           10000

Now, we can just create a simple example polygon to restrict our region of interest from -180 to 180 to -90 to 90 longitude.

library(sf)

pl <- list(rbind(c(-90,-90), c(-90,90), c(90,90), c(90,-90), c(-90,-90)))
roi <- st_sfc(st_polygon(pl), crs = "EPSG:4326")

crop(v, roi)
#> class       : SpatRaster 
#> dimensions  : 100, 50, 1  (nrow, ncol, nlyr)
#> resolution  : 3.6, 1.8  (x, y)
#> extent      : -90, 90, -90, 90  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#> source      : memory 
#> name        : filedf6b216a737 
#> min value   :              26 
#> max value   :            9975
6
  • 4
    Just adding a comment to highlight the fact tha VRT is a gdal stuff not terra and is supported by raster, terra and stars
    – JRR
    Commented Nov 30, 2021 at 23:34
  • 1
    Very good point @JRR. While supported by others, to my understanding terra offers functionality to create a VRT while the other packages only support the use of it. Looking around, it seems gdalUtils also supports building a VRT.
    – caldwellst
    Commented Dec 1, 2021 at 8:17
  • 2
    stars, terra and sf can call gdalbuildvrt. I'm not sure for raster
    – JRR
    Commented Dec 1, 2021 at 12:34
  • Interesting @JRR. How do you call gdalbuildvrt through those packages? For instance, I know sf has a set of gdal functions but does not include gdalbuildvrt. How else would you access it?
    – caldwellst
    Commented Dec 1, 2021 at 14:22
  • 3
    sf::gdal_utils("buildvrt", ...), terra::vrt(...), stars::st_mosaic()
    – JRR
    Commented Dec 1, 2021 at 14:47

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.