terra包mosaic函数批量栅格镶嵌异常，输出结果异常求助

问题背景 #

需要合并一批CRS、分辨率一致但空间范围不同的栅格，重叠区域要求计算均值。使用terra::mosaic函数批量处理所有栅格时输出异常，但以下场景可正常运行：

• 仅处理部分栅格子集
• 分两步执行镶嵌操作
• 使用merge函数替代mosaic
• 将镶嵌函数从默认的mean改为first等其他函数

最小可复现代码 #

library(terra)

# 创建测试栅格
x1 <- rast(xmin=3323433, xmax=3668916, ymin=5774887, ymax=6119500, res=1000, vals=1, crs=terra::crs("EPSG:5677"))
x2 <- rast(xmin=3478201, xmax=3882596, ymin=5757906, ymax=6095527, res=1000, vals=2, crs=terra::crs("EPSG:5677"))
x3 <- rast(xmin=3278500, xmax=3625980, ymin=5410296, ymax=5828826, res=1000, vals=3, crs=terra::crs("EPSG:5677"))
x4 <- rast(xmin=3555086, xmax=3941507, ymin=5550139, ymax=5916728, res=1000, vals=4, crs=terra::crs("EPSG:5677"))
x5 <- rast(xmin=3459229, xmax=3875606, ymin=5228500, ymax=5661015, res=1000, vals=5, crs=terra::crs("EPSG:5677"))
x6 <- rast(xmin=3373358, xmax=3664922, ymin=5244482, ymax=5629050, res=1000, vals=6, crs=terra::crs("EPSG:5677"))

x  <- list(x1, x2, x3, x4, x5, x6)
x_ <- list(x1, x2, x3, x4)

1. 批量全量处理（异常） #

collection <- terra::sprc(x)
m1 <- mosaic(collection)
plot(m1)

2. 分步处理（正常） #

# 第一步合并前4个栅格
collection <- terra::sprc(x_)
m2 <- mosaic(collection)
plot(m2)

# 第二步合并剩余栅格
m3 <- mosaic(m2, x5, x6)
plot(m3)

3. 使用merge批量处理（正常） #

collection <- terra::sprc(x)
m4 <- merge(collection)
plot(m4)

运行环境 #

R version 4.1.2 (2021-11-01)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 22.04.3 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=de_DE.UTF-8        LC_COLLATE=en_US.UTF-8     LC_MONETARY=de_DE.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=de_DE.UTF-8       LC_NAME=C                  LC_ADDRESS=C               LC_TELEPHONE=C             LC_MEASUREMENT=de_DE.UTF-8 LC_IDENTIFICATION=C        

attached base packages:
[1] parallel  stats     graphics  grDevices utils     datasets  methods   base      

other attached packages:
[1] terra_1.7-71      dplyr_1.1.4       doParallel_1.0.17 iterators_1.0.14  foreach_1.5.2    

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.12        pillar_1.9.0       compiler_4.1.2     class_7.3-20       tools_4.1.2        xts_0.13.2         gstat_2.1-1        lifecycle_1.0.4    tibble_3.2.1      
[10] lattice_0.20-45    pkgconfig_2.0.3    rlang_1.1.3        DBI_1.2.2          cli_3.6.2          e1071_1.7-14       generics_0.1.3     vctrs_0.6.5        hms_1.1.3         
[19] classInt_0.4-10    grid_4.1.2         tidyselect_1.2.0   spacetime_1.3-1    glue_1.7.0         sf_1.0-15          R6_2.5.1           fansi_1.0.6        sp_2.1-3          
[28] tzdb_0.4.0         readr_2.1.5        magrittr_2.0.3     intervals_0.15.4   codetools_0.2-18   units_0.8-5        utf8_1.2.4         KernSmooth_2.23-20 proxy_0.4-27      
[37] FNN_1.1.4          zoo_1.8-12 

异常原因分析 #

该问题属于特定Terra版本（1.7-71）的已知bug：当mosaic函数处理多栅格（数量超过一定阈值）且使用mean作为合并函数时，重叠区域的均值累加计算逻辑存在边界错误。复杂的空间重叠层级会触发该bug，导致部分区域的均值计算异常；而分步处理时每次处理的栅格数量少，中间结果的累加逻辑可正常运行；merge函数的重叠区处理逻辑与mosaic不同，不受此bug影响；first等简单取数函数无需累加计算，因此不会触发错误。

解决办法 #

1. 分步镶嵌：将栅格列表拆分为若干批次，分多次执行mosaic操作，最后合并结果
2. 替换为merge函数：使用terra::merge替代mosaic，并指定fun=mean实现重叠区均值计算：

   m <- merge(collection, fun=mean)
   

3. 升级Terra版本：该bug已在Terra 1.7.x后续版本中修复，建议升级至最新稳定版：

   install.packages("terra")
   

4. 手动实现均值镶嵌：创建目标范围的空栅格，逐个累加栅格值并计数，最后计算均值：

   # 创建覆盖所有栅格范围的空栅格
   ext_all <- ext(do.call(c, lapply(x, ext)))
   r_template <- rast(ext_all, res=1000, crs=crs(x1))
   r_sum <- r_template
   r_count <- r_template
   
   # 累加值与计数
   for (ras in x) {
     r_sum <- cover(r_sum, ras) + ifelse(is.na(r_sum), 0, r_sum)
     r_count <- cover(r_count, 1) + ifelse(is.na(r_count), 0, r_count)
   }
   
   # 计算均值
   r_mean <- r_sum / r_count
   plot(r_mean)
   

内容的提问来源于stack exchange，提问作者JKupzig