在 R 中处理大型光栅马赛克，无需将它们合并到单个文件（如 lidR 目录）
Work with large raster mosaics in R without merging them to a single file (like lidR catalog)
lidR 包有一种处理大型（点云）数据集的巧妙方法：catalog 函数（此处为 doc）避免将数据集加载到内存中，并且可以处理马赛克 [数据集是分布在多个（非重叠）图块]作为单个数据集。它以智能方式在计算过程中即时加载所需的图块。如果只处理数据集的一小部分，最好避免处理大文件（多个 GB）并保持内存需求精简。
在 R 中是否有类似的 convenient/memory-efficient/"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 包，但它似乎没有提供任何此类功能。
您应该查看 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 到 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
lidR 包有一种处理大型（点云）数据集的巧妙方法：catalog 函数（此处为 doc）避免将数据集加载到内存中，并且可以处理马赛克 [数据集是分布在多个（非重叠）图块]作为单个数据集。它以智能方式在计算过程中即时加载所需的图块。如果只处理数据集的一小部分，最好避免处理大文件（多个 GB）并保持内存需求精简。
在 R 中是否有类似的 convenient/memory-efficient/"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 包，但它似乎没有提供任何此类功能。
您应该查看 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 到 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