在R中使用镶嵌图来合并多个geotiff

Question

我有50个geotiff文件在同一个文件夹中。 它们都代表了世界不同地区的高程数据。 我想合并某些geotiff文件，我发现R中的mosaic可能会帮助我们。 我已经将这些geotiff移动到了同一个文件夹中，并编写了下面的R脚本：

setwd()
a<-list.files(pattern="*.tiff",file.name=TRUE)
combind<-merge(a,fun=mean)

然而，这个脚本返回了一个错误：error in as.data.frame(y)

请问我应该如何改进我的脚本？

- Bing-Hong Huang

---

你可能需要使用 do.call(merge, a)，而且 merge 函数没有 fun 参数。 - akrun

谢谢您的回复，但它返回“第二个参数必须是列表”。 - Bing-Hong Huang

您没有提供示例。它可能是Reduce(function(...) merge(...), a)。 - akrun

顺便说一下，“merge”来自“raster”包。 - Bing-Hong Huang

你需要提供一个可重现的例子。 - akrun

如果您查看raster中的?merge，则示例可以使用do.call(merge, x)运行，并且我也没有找到任何fun参数。 - akrun

Answer 1

编辑2023-06：使用terra::vrt()或terra::mosaic()

由于gdalUtils似乎已经有一段时间没有维护了，您可以使用与我原始答案类似的方法，只需使用更为更新的terra包，该包是raster包的继任者：

library(terra)

vrt(
  x = list.files(path = "folder/to/images", pattern = "*.tif$", full.names = TRUE), 
  filename = "dem.vrt"
)
# afterwards read it as if it was a normal raster:
dem <- rast("dem.vrt")

另一种可能性是使用terra::mosaic()。这样可以更好地控制重叠像素的聚合函数。首先，您需要使用sprc()从rast对象列表创建一个栅格集合，然后使用mosaic将它们合并成一个栅格。

library(terra)
# vector of file names
fls <- list.files("your/folder", ".tif$", full.names = TRUE)
# list of rast objects
r_lst <- lapply(fls, rast)
# create spatial raster collection
coll <- sprc(r_list)

# combine all rasters
mosaic(coll, function = "mean")

使用gdalUtils的原始答案

您可以利用强大的GDAL函数。根据我的经验，这些函数比纯R代码要快得多。

我的方法是使用library(gdalUtils)：

首先，构建一个虚拟栅格文件（vrt）：

library(gdalUtils)
setwd(...)
gdalbuildvrt(gdalfile = "*.tif", # uses all tiffs in the current folder
             output.vrt = "dem.vrt")

然后，将虚拟栅格复制到实际物理文件中：

gdal_translate(src_dataset = "dem.vrt", 
               dst_dataset = "dem.tif", 
               output_Raster = TRUE # returns the raster as Raster*Object
                                    # if TRUE, you should consider to assign 
                                    # the whole function to an object like dem <- gddal_tr..
               options = c("BIGTIFF=YES", "COMPRESSION=LZW"))

另一个纯粹（可能更慢）的光栅包解决方案是：

f <- list.files(path = "your/path", pattern = ".tif$", full.names = TRUE)
rl <- lapply(f, raster)

do.call(merge, c(rl, tolerance = 1))

你需要调整“容差”，因为光栅文件可能没有相同的起点。

- loki

---

感谢您的回复。gdalUtils最终返回警告信息""C:\Program Files (x86)\QGIS 2.18\bin\gdal_translate.exe" -of "GTiff" "dem.vrt" "elevation.tif"' had status 1.因此，我正在尝试第二个代码。 - Bing-Hong Huang

也许你需要设置 options = c("BIGTIFF=YES")。 - loki

我有时间检查了一下。在我的情况下，gdalbuildvrt（gdalfile =“.tif $”...）失败了。我编辑了我的答案。也许你可以检查一下。只需将“.tif $”替换为“.tif”即可。我还添加了一些选项。 - loki