在 R 中处理大型栅格马赛克，无需将它们合并到单个文件（如lidR目录）

Question

litR 包有一种巧妙的方法来处理巨大（点云）数据集：该catalog函数（此处为doc）避免将数据集加载到内存中，并且可以将马赛克（分布在多个（不重叠）图块上的数据集）视为单个数据集。它在计算过程中以智能方式动态加载所需的图块。如果只处理数据集的一小部分，最好避免使用大文件（多个 GB）并保持内存需求精简。

R 中是否有类似的方便/内存高效/“lidR-catalog-way”来处理大型光栅马赛克？ 或者用更通用的方式来说： 是否有一种方法可以在 R 中使用镶嵌栅格数据集而无需先合并它们？

我知道mosaic（doc）和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)
  }

然而，两者都需要将整个数据集具体化为内存中（或至少在磁盘上）的单个文件。有办法避免这种情况吗？

我的应用程序： 我正在使用一个巨大的 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)

我通常使用raster package，但这似乎没有提供任何此类功能。

Answer 1

您应该研究一下terra通过虚拟光栅图块 (VRT) 提供您正在寻找的功能的软件包。我们可以使用它们将磁盘上的光栅文件集合视为单个光栅文件，同时利用 API 执行与通过包执行的大部分相同任务raster。

首先，让我们使用直接来自文档的示例创建 4 个栅格的样本?terra::vrt()。

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 到 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