如何基于fitQmapQUANT批量处理366层栅格栈的偏差校正与绘图？

前提确认 #

确保prp_r（RCM观测数据）和prp_g（GPM待校正数据）的栅格栈图层顺序完全对应（第i层为同一天的观测/模拟数据），且空间范围、分辨率一致，可通过compareRaster(prp_r, prp_g)验证。

方案1：For循环（直观易调试） #

适合新手理解，方便中途调试查看结果：

library(raster)
library(qmap)

# 初始化输出栅格栈，复用GPM的空间参数与图层数量
prp_corrected <- stack(prp_g)

# 遍历每个图层
for (i in 1:nlayers(prp_g)) {
  # 提取当日的RCM和GPM单图层
  r_layer <- prp_r[[i]]
  g_layer <- prp_g[[i]]
  
  # 转换为带坐标的数据框并合并（避免xy列重复）
  r_df <- as.data.frame(r_layer, xy = TRUE)
  g_df <- as.data.frame(g_layer, xy = TRUE)
  combined_df <- cbind(r_df, gpm_pr = g_df[, 3])
  
  # 拟合偏差校正模型
  qm_fit <- fitQmapQUANT(obs = combined_df[, 3], 
                         mod = combined_df$gpm_pr, 
                         qstep = 0.01, 
                         nboot = 1, 
                         wet.day = TRUE)
  
  # 执行校正并转回栅格
  corrected_pr <- doQmapQUANT(combined_df$gpm_pr, qm_fit, type = "tricub")
  corrected_layer <- raster(r_layer)
  values(corrected_layer) <- corrected_pr
  
  # 赋值到输出栅格栈
  prp_corrected[[i]] <- corrected_layer
  
  # 可选：打印进度
  if (i %% 10 == 0) cat("完成第", i, "个图层\n")
}

方案2：Lapply函数式编程（简洁高效） #

用函数式编程替代循环，代码更紧凑：

library(raster)
library(qmap)

# 定义单图层处理函数
process_single_layer <- function(i) {
  r_layer <- prp_r[[i]]
  g_layer <- prp_g[[i]]
  
  r_df <- as.data.frame(r_layer, xy = TRUE)
  g_df <- as.data.frame(g_layer, xy = TRUE)
  combined_df <- cbind(r_df, gpm_pr = g_df[, 3])
  
  qm_fit <- fitQmapQUANT(obs = combined_df[, 3], 
                         mod = combined_df$gpm_pr, 
                         qstep = 0.01, 
                         nboot = 1, 
                         wet.day = TRUE)
  
  corrected_pr <- doQmapQUANT(combined_df$gpm_pr, qm_fit, type = "tricub")
  corrected_layer <- raster(r_layer)
  values(corrected_layer) <- corrected_pr
  
  return(corrected_layer)
}

# 批量处理并转为栅格栈
corrected_layer_list <- lapply(1:nlayers(prp_g), process_single_layer)
prp_corrected <- stack(corrected_layer_list)

方案3：矩阵化处理（大数据场景更高效） #

减少栅格与数据框的反复转换，适合大尺寸栅格栈：

library(raster)
library(qmap)

# 将栅格栈转为矩阵：ncell行，nlayers列（每行对应一个单元格，每列对应一个日期）
rcm_mat <- as.matrix(prp_r)
gpm_mat <- as.matrix(prp_g)

# 初始化校正结果矩阵
corrected_mat <- matrix(NA, nrow = nrow(rcm_mat), ncol = ncol(rcm_mat))

# 遍历每一列（日期）
for (i in 1:ncol(rcm_mat)) {
  rcm_day <- rcm_mat[, i]
  gpm_day <- gpm_mat[, i]
  
  # 拟合模型并校正
  qm_fit <- fitQmapQUANT(obs = rcm_day, 
                         mod = gpm_day, 
                         qstep = 0.01, 
                         nboot = 1, 
                         wet.day = TRUE)
  corrected_day <- doQmapQUANT(gpm_day, qm_fit, type = "tricub")
  
  corrected_mat[, i] <- corrected_day
  
  if (i %% 10 == 0) cat("完成第", i, "个图层\n")
}

# 将矩阵转回栅格栈（注意转置：raster的values需要按图层顺序填充）
prp_corrected <- stack(prp_g)
values(prp_corrected) <- t(corrected_mat)

后续绘图示例 #

校正完成后，可批量绘制栅格图：

# 用rasterVis包绘制多图层面板图
library(rasterVis)
levelplot(prp_corrected[[1:3]], main = "校正后降水栅格（前3日）")

# 或用基础plot函数逐个绘制
plot(prp_corrected[[1]], main = "第1日校正后降水")

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