使用Terra包镶嵌多波段Sentinel-2栅格时所有波段值异常一致

问题重现 #

处理相邻的Sentinel-2分块影像（均为7波段SpatRaster）时，使用terra::mosaic()得到的结果所有波段极值完全相同，可视化效果一致；但使用terra::merge()可得到符合预期的结果（保留各波段原始极值差异）。

原始栅格信息 #

> sent1
class       : SpatRaster 
dimensions  : 10980, 10980, 7  (nrow, ncol, nlyr)
resolution  : 10, 10  (x, y)
extent      : 399960, 509760, 6290220, 6400020  (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / UTM zone 9N (EPSG:32609) 
source      : T09VVD_2022_mean.tif 
names       :   B02,  B03,   B04,        B05,        B06,       B07, ... 
min values  :     0,    0,    29,   379.8125,   523.3125,   385.875, ... 
max values  : 10456, 9464, 15768, 13930.6875, 15912.1250, 16714.938, ... 

> sent2
class       : SpatRaster 
dimensions  : 10980, 10980, 7  (nrow, ncol, nlyr)
resolution  : 10, 10  (x, y)
extent      : 399960, 509760, 6390240, 6500040  (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / UTM zone 9N (EPSG:32609) 
source      : T09VVE_2022_mean.tif 
names       :  B02,  B03,   B04,       B05,        B06,      B07, ... 
min values  :    0,  211,     1,   379.125,   331.3125,   342.75, ... 
max values  : 9104, 8448, 10480, 12210.500, 16052.9375, 16814.38, ... 

mosaic()异常结果 #

sent3 <- mosaic(sent1, sent2)

> show(sent3)
class       : SpatRaster 
dimensions  : 20982, 10980, 7  (nrow, ncol, nlyr)
resolution  : 10, 10  (x, y)
extent      : 399960, 509760, 6290220, 6500040  (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / UTM zone 9N (EPSG:32609) 
source      : spat_fubTXtzTOH0dho5_70406.tif 
names       :        B02,        B03,        B04,        B05,        B06,        B07, ... 
min values  :   355.6161,   355.6161,   355.6161,   355.6161,   355.6161,   355.6161, ... 
max values  : 12130.0361, 12130.0361, 12130.0361, 12130.0361, 12130.0361, 12130.0361, ... 

即使将栅格加入SpatRasterCollection后镶嵌，结果仍一致：

sent.sprc <- sprc(sent1, sent2)
sent4 <- mosaic(sent.sprc)

> show(sent4)
class       : SpatRaster 
dimensions  : 20982, 10980, 7  (nrow, ncol, nlyr)
resolution  : 10, 10  (x, y)
extent      : 399960, 509760, 6290220, 6500040  (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / UTM zone 9N (EPSG:32609) 
source      : spat_hDQlRCrwBUOApYB_49105.tif 
names       :        B02,        B03,        B04,        B05,        B06,        B07, ... 
min values  :   355.6161,   355.6161,   355.6161,   355.6161,   355.6161,   355.6161, ... 
max values  : 12130.0361, 12130.0361, 12130.0361, 12130.0361, 12130.0361, 12130.0361, ... 

merge()正常结果 #

sent5 <- merge(sent1, sent2)

> show(sent5)
class       : SpatRaster 
dimensions  : 20982, 10980, 7  (nrow, ncol, nlyr)
resolution  : 10, 10  (x, y)
extent      : 399960, 509760, 6290220, 6500040  (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / UTM zone 9N (EPSG:32609) 
source      : spat_fqpzxaTDk0yKNGb_50974.tif 
names       :   B02,  B03,   B04,       B05,        B06,      B07, ... 
min values  :     0,    0,     1,   379.125,   331.3125,   342.75, ... 
max values  : 10456, 9464, 15768, 13930.688, 16052.9375, 16814.38, ... 

环境信息 #

• 运行环境：HPC及本地机器
• Terra版本：开发版
• 加载包：仅terra和here

问题原因分析 #

1. Terra开发版潜在bug：mosaic()错误地跨波段统一计算数值（如全局均值/极值），而非按波段独立处理，导致所有波段特征完全一致。这属于多波段镶嵌逻辑的异常。
2. 默认函数参数的异常执行：mosaic()默认fun=mean，正常应按波段独立计算均值，但异常结果显示函数执行时混淆了波段维度的处理逻辑。

解决方法 #

方法1：用merge()替代mosaic() #

对于相邻且重叠范围小的Sentinel-2分块影像，merge()更适配此类拼接场景，能正确保留各波段原始数值特征，如测试结果所示。

方法2：显式指定mosaic()的按波段处理函数 #

强制函数按像素取第一个非NA值（适用于重叠区域无需均值的场景），确保按波段独立处理：

sent3 <- mosaic(sent1, sent2, fun = first)

或自定义函数：

sent3 <- mosaic(sent1, sent2, fun = function(x) x[!is.na(x)][1])

方法3：切换到Terra稳定版 #

开发版可能存在未修复的逻辑问题，安装并使用最新稳定版尝试解决：

install.packages("terra")
library(terra)

方法4：提交bug报告 #

若上述方法无效，可将测试文件（或最小复现示例）提交到Terra的官方代码仓库，帮助开发者定位修复问题。

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