在R中处理大型光栅镶嵌图，而不将它们合并到一个文件中(如lidR目录)

Question

lidR包有一种处理大型(点云)数据集的简洁方法：catalog函数(此处的文档)避免将数据集加载到内存中，并且可以将分布在多个(不重叠)块上的马赛克数据集视为单个数据集。它以一种智能的方式在计算过程中加载所需的瓷砖.避免使用大型文件(多个GBs)，如果只处理数据集的小部分，则保持内存需求的低水平是很好的。

是否有一个类似的convenient/memory-efficient/"lidR-catalog-way“来处理R?或更多的大型光栅拼接：是否有一种方法可以处理R中的拼接光栅数据集而不首先合并它们？

我知道mosaic (文档)和merge函数，它们允许我将平铺的光栅马赛克合并成一个光栅数据集。我还发现，gdal将这样做，的速度更快，内存效率更高。以下是这方面的一个R片段：

  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)
  }

但是，两者都需要将整个数据集具体化为内存中的单个文件(或至少在磁盘上)。有什么办法可以避免这种情况吗？

My应用程序：我正在处理一个巨大的LAS点云(几个TB，所有东西都被平铺在100个MB文件中)和相应的栅格数据集(大约100 GB)，例如一个高分辨率地形模型。我通常处理小部分(通过空间多边形..shp/..kml分隔)或整个LAS瓷砖。这在lidR中非常节省内存，只需要加载所需的块：

# 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))

我也想对光栅数据集做同样的处理：

# 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)

我通常使用光栅包，但这似乎没有提供任何这样的功能。

Answer 1

您应该查看terra包，它提供了您通过虚拟光栅贴片(VRT)所要寻找的功能。我们可以使用它们将磁盘上的光栅文件集合作为单个光栅文件来处理，同时利用API执行与通过raster包执行的大部分任务相同的任务。

首先，让我们直接使用?terra::vrt()文档中的示例创建一个由4个光栅组成的示例。

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"

现在，我们将以VRT的形式阅读它们，再一次，直接从同一个示例中读取它们。这允许

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

现在，我们可以创建一个简单的多边形来限制我们感兴趣的区域从-180度到-90度到90度。

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