栅格气候数据偏差校正（qmap）报错：公式无法使用

我来帮你排查这个问题，先还原你的场景和报错：

你尝试用R的qmap包对栅格化的气候观测与模拟数据做偏差校正，但运行代码时遇到了如下错误：

Error in (function (x, fun, filename = "", recycle = TRUE, forcefun = FALSE, : cannot use this formula, probably because it is not vectorized

你的原始代码如下：

library(raster)
library(qmap)
#Create a rasterStack with observed and modelled data
r <- raster(ncol=20, nrow=20)
obs <- stack(lapply(1:100, function(x) setValues(r, runif(ncell(r))))) #observed data
mod <- stack(lapply(1:100, function(x) setValues(r, runif(ncell(r)))))*2 #modelled data (i want this unbiased)
#bias-correction function
f <- function(obs, mod, ...) {
 obs <- t(obs)
 qm.fit <- fitQmap(obs, t(mod), method="QUANT",qstep=0.01)
 t(doQmap(mod, qm.fit, type="linear") )
}
x <- overlay(obs, mod, fun=f)

错误原因分析 #

overlay函数处理栅格栈时，会把每个像元的所有层值整理成矩阵形式传递给自定义函数：比如你的obs和mod各有100层，那么传递给f的obs是ncell(r) × 100的矩阵（每行对应一个像元的100个时间步值），mod同理。

你的函数有两个核心问题：

1. 转置操作t(obs)打乱了overlay传递的矩阵维度，导致fitQmap接收的输入不符合「每行对应一个时间序列」的要求；
2. 函数没有适配overlay要求的向量化逻辑——即输入矩阵时，要返回行数一致的结果矩阵（每个像元对应一行输出）。

修正后的解决方案 #

我们需要调整函数逻辑，让它能处理overlay的矩阵输入，对每个像元的时间序列独立做偏差校正。这里提供两种高效的实现方式：

方式1：适配overlay的矩阵输入

library(raster)
library(qmap)

# 创建测试数据
r <- raster(ncol=20, nrow=20)
obs <- stack(lapply(1:100, function(x) setValues(r, runif(ncell(r))))) # 观测数据
mod <- stack(lapply(1:100, function(x) setValues(r, runif(ncell(r))))) * 2 # 模拟数据（存在偏差）

# 调整后的偏差校正函数
bias_correct_fun <- function(obs_mat, mod_mat, ...) {
  # obs_mat: ncell × nlayers 矩阵，每行是一个像元的观测时间序列
  # mod_mat: ncell × nlayers 矩阵，每行是一个像元的模拟时间序列
  
  # 逐行处理每个像元的时间序列
  corrected <- t(apply(cbind(obs_mat, mod_mat), 1, function(row) {
    obs_ts <- row[1:100] # 提取当前像元的观测序列
    mod_ts <- row[101:200] # 提取当前像元的模拟序列
    
    # 拟合偏差校正模型并应用
    qm_fit <- fitQmap(obs_ts, mod_ts, method="QUANT", qstep=0.01)
    doQmap(mod_ts, qm_fit, type="linear")
  }))
  
  return(corrected)
}

# 执行偏差校正
x <- overlay(obs, mod, fun=bias_correct_fun)

方式2：用calc合并处理（更直观）

把观测和模拟数据合并成一个栈，用calc逐像元处理，逻辑更清晰：

# 合并观测与模拟数据为一个栈
combined_stack <- stack(obs, mod)

# 用calc逐像元执行偏差校正
x <- calc(combined_stack, fun=function(vec) {
  obs_ts <- vec[1:100] # 前100个值是观测时间序列
  mod_ts <- vec[101:200] # 后100个值是模拟时间序列
  
  qm_fit <- fitQmap(obs_ts, mod_ts, method="QUANT", qstep=0.01)
  doQmap(mod_ts, qm_fit, type="linear")
})

改动说明 #

1. 适配了栅格函数的输入格式：每个像元的时间序列作为一维向量传入，符合fitQmap对输入的要求；
2. 逐像元独立拟合校正模型，这是气候栅格数据偏差校正的常规逻辑（每个栅格像元视为独立的格点）；
3. 避免了原函数中转置操作导致的维度不匹配问题，确保整个流程的向量/矩阵维度一致。

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