使用R中的terra包拼接不同范围的栅格图出现错误？

Question

使用R中的terra包拼接不同范围的栅格图出现错误？

4

我有一个文件夹，其中包含约191个GeoTIFF文件（每个文件都是更大区域的不同DEM（海拔）瓦片）。我想将所有瓦片合并成一个栅格文件。我使用terra包，并成功地能够加载每个栅格并将它们从2米分辨率聚合到30米分辨率。然而，当运行mosaic函数来合并它们时，遇到了错误（请参见下面的错误消息）。我已经能够在仅三个瓦片的较小子集上运行马赛克函数，但当我扩展到所有文件时，这就成为一个问题。

通过调用栅格总结（见下文），聚合确实会略微改变范围-这可能是问题吗？resample可能是一种选择，但每个单独的栅格具有不同的范围，我不确定如何实施这个修复方法。

不确定样本数据集是否有帮助，因为我知道函数有效。我在高性能集群上运行此代码，因此运行小批量代码效率不高。

library(terra)

files <- as.list(list.files("./DEM_tiles", full.names = TRUE))

raster.list <- lapply(files, rast)
 
for(i in 1:length(raster.list)){
 raster.list[[i]] <- aggregate(raster.list[[i]], fact = 15)
 }

raster.mosaic <- do.call(mosaic, raster.list)

> Error: [mosaic] internal error: extents do not match ()
> Execution halted

以下是两个瓷砖的示例：

### Before Aggregation
[[1]]
class       : SpatRaster 
dimensions  : 25000, 25000, 1  (nrow, ncol, nlyr)
resolution  : 2, 2  (x, y)
extent      : -1800000, -1750000, -6e+05, -550000  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
source      : 35_23_1_1_2m_v3.0_reg_dem.tif 
name        : 35_23_1_1_2m_v3.0_reg_dem 

[[2]]
class       : SpatRaster 
dimensions  : 25000, 25000, 1  (nrow, ncol, nlyr)
resolution  : 2, 2  (x, y)
extent      : -1800000, -1750000, -550000, -5e+05  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
source      : 35_23_1_2_2m_v3.0_reg_dem.tif 
name        : 35_23_1_2_2m_v3.0_reg_dem 


### After Aggregation
[[1]]
class       : SpatRaster 
dimensions  : 1667, 1667, 1  (nrow, ncol, nlyr)
resolution  : 30, 30  (x, y)
extent      : -1800000, -1749990, -600010, -550000  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
source      : memory 
name        : 35_23_1_1_2m_v3.0_reg_dem 
min value   :                 -15.62178 
max value   :                  233.6489 

[[2]]
class       : SpatRaster 
dimensions  : 1667, 1667, 1  (nrow, ncol, nlyr)
resolution  : 30, 30  (x, y)
extent      : -1800000, -1749990, -550010, -5e+05  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
source      : memory 
name        : 35_23_1_2_2m_v3.0_reg_dem 
min value   :                 -15.27713 
max value   :                  243.0772

- apple

3

一种方法是建立一个包含所有瓦片的vrt文件，然后将其写成tif格式。 - undefined

我是在寻找解决同样错误的方法时来到这里的。请确保 files <- as.list(list.files("./DEM_tiles", full.names = TRUE)) 只读取预期的文件。同一个文件夹中可能还有其他不应该合并的文件。 - undefined

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Robert Hijmans · Accepted Answer

错误是由于栅格不对齐和一个错误导致的。现在我得到了。

library(terra)
#terra version 1.2.1
crs <- "+proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +x_0=0 +y_0=0 +datum=WGS84 +units=m"
r1 <- rast(nrow=1667, ncol=1667, ext=c(-1800000, -1749990, -600010, -550000), crs=crs)
r2 <- rast(nrow=1667, ncol=1667, ext=c(-1800000, -1749990, -550010, -5e+05), crs=crs)
values(r1) <- 1:ncell(r1)
values(r2) <- 1:ncell(r2)
m <- mosaic(r1, r2)
#Warning message:
#[mosaic] rasters did not align and were resampled

m
#class       : SpatRaster 
#dimensions  : 3334, 1667, 1  (nrow, ncol, nlyr)
#resolution  : 30, 30  (x, y)
#extent      : -1800000, -1749990, -600010, -499990  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
#source      : memory 
#name        :   lyr.1 
#min value   :       1 
#max value   : 2778889

同时

@Elia的建议是一个不错的解决方案：

r1 <- writeRaster(r1, "test1.tif", overwrite=TRUE)
r2 <- writeRaster(r2, "test2.tif", overwrite=TRUE)
v <- vrt(c("test1.tif", "test2.tif"), "test.vrt", overwrite=TRUE)

v
#class       : SpatRaster 
#dimensions  : 3334, 1667, 1  (nrow, ncol, nlyr)
#resolution  : 30, 30  (x, y)
#extent      : -1800000, -1749990, -600020, -5e+05  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
#source      : test.vrt 
#name        :    test 
#min value   :       1 
#max value   : 2778889