如何让terra::app()按顺序将函数列表应用于多波段栅格？

问题描述 #

使用R语言的terra和RStoolbox包处理Sentinel 2数据时，调用histMatch()函数对两幅影像做色彩平衡，出现异常结果：直方图匹配后影像的波段最大值有时会小于参考影像对应波段的最小值。已在RStoolbox的GitHub仓库提交该bug（issue #117）。

问题定位 #

通过分析源码，定位到问题根源：histMatch()会为参考栅格的每个波段生成对应函数，用于在累积分布函数中找到特定分位数对应的原始值，随后创建函数列表并尝试按顺序应用到待处理栅格的对应波段。

以下是histMatch()函数的相关代码片段：

totalFun <- function(xvals, f = layerFun) {
    if (is.vector(xvals)) 
        xvals <- as.matrix(xvals)
    app <- lapply(1:ncol(xvals), function(i) {
        f[[i]](xvals[, i])
    })
    do.call("cbind", app)
}
if (returnFunctions) {
    names(layerFun) <- names(x)
    return(layerFun)
}
.vMessage("Apply histogram match functions")
out <- terra::app(x, fun = totalFun, ...)

核心问题：terra::app()不会遵循传入函数的顺序，导致函数与波段不匹配，进而出现异常结果。

我正尝试编写修复方案提交到该bug报告中，为此制作了简化的复现示例，展示期望输出和目标实现方式：

####Install and/or load terra####
if(!require(terra)){
  install.packages(terra)
}


####Example Data Setup####
set.seed(23)#For reproducibility

##Creating a Fake multiband raster##
r1a<-rast(xmin=-123.1, xmax=-122.1, ymin=45, ymax=45.5, nrows=100, ncols=100, crs="epsg:4326")
r1b<-rast(xmin=-123.1, xmax=-122.1, ymin=45, ymax=45.5, nrows=100, ncols=100, crs="epsg:4326")
r1c<-rast(xmin=-123.1, xmax=-122.1, ymin=45, ymax=45.5, nrows=100, ncols=100, crs="epsg:4326")

values(r1a)<-rnorm(10000, 200, 35)
values(r1b)<-runif(10000, 500, 10000)
values(r1c)<-rnorm(10000, 600, 150)

r1<-c(r1a, r1b, r1c)

##Creating a list of functions
fxn1<-function(x){log(x+1)}
fxn2<-function(x){(x)^0.5}
fxn3<-function(x){x-25}

funtime<-list(fxn1, fxn2, fxn3)

####Applying each function to its corresponding band in the fake raster####
##Desired behavior##
out<-list()
for(i in 1:nlyr(r1)){
  out[[i]]<-app(r1[[i]], funtime[[i]])
}

do.call(c, out)

####Desired approach (throws error)####
totalfun<-function(x, f=funtime){
  lapply(1:length(x), function(i){
  out<-f[[i]](x[[i]])
})
  do.call(c, out)
}

result<-app(r1, fun=totalfun)

求助需求 #

希望能得到该问题的解决方法，以便完成修复方案并提交到GitHub的bug报告中，帮助其他进行多波段影像直方图匹配的用户。谢谢！

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