在R中高效处理依赖型Raster Brick计算，优化循环性能

问题背景 #

我有一个Raster Brick（ncell=28536，nlayers=181），需要基于原始栅格砖运行数学函数，生成两个尺寸相同且相互依赖的栅格砖：

• outputBrick1的第i层基于outputBrick2的第i层计算
• outputBrick2的第i层（i>1）基于outputBrick1的第i-1层计算

当前代码在小数据集（24个像元、181层）上运行正常，但在28000个像元的数据集上速度极慢。尝试调整内存设置后，CPU和RAM使用率仍极低（CPU约6%、RAM约10%，R仅占用5%CPU和1GB内存），系统配置为64GB RAM、16GB GPU。

原始代码 #

初始化代码 #

library(raster)
b <- brick(ncols=5, nrows=5, nl=5)
inBrick <- setValues(b, runif(ncell(b) * nlayers(b)))
inBrick[c(1,2,3,22,23,24,25)] <- NA
outBrick1 <- inBrick
outBrick1[] <- NA
outBrick2 <- outBrick1

ini <- 0.3
p <- 0.15
p1 <- p/3
p2 <- p-(p/3)
fc <- 0.3
var1 <- which(!is.na(inBrick[[1]][]))

outBrick2[[1]][var1] <- ini
### now outBrick2 has initial values in 1st layer
weather <- c(0.1, 0, 0, 0, 0.3)

低效的循环计算代码 #

var3 <- 1:ncell(inBrick)
### outBrick1 Calculations
for (i in 1:nlayers(inBrick)) {
varr1 <- inBrick[[i]][]*(((outBrick2[[i]][]-p1)/(p2))^2)
for (j in 1:ncell(inBrick)) {
if(!is.na(outBrick2[[i]][j])){
  if(outBrick2[[i]][j]>p){
     outBrick1[[i]][j] <- inBrick[[i]][j]
  }else{
     outBrick1[[i]][j] <- varr1[j]
    }
  }
}
###outBrick2 Calculations
for (k in 2:nlayers(inBrick)) {
var2 <- outBrick2[[k-1]][] + (weather[k-1]-outBrick1[[k-1]][])/100
 for(l in 1:ncell(inBrick)){
  var3[l] <- min(fc, var2[l])
 }
 outBrick2[[k]][] <- var3
 }
}

内存设置尝试 #

rasterOptions(maxmemory = 5.17e+10)
rasterOptions(memfrac = 0.8)
rasterOptions(chunksize = 5.17e+10)

优化方案 #

核心思路是用向量化操作替代嵌套循环，R的原生循环（尤其是逐像元循环）效率极低，向量化函数基于底层编译语言实现，能大幅提升速度，同时充分利用内存。

优化后的代码 #

library(raster)

# 初始化部分（保留原逻辑）
b <- brick(ncols=5, nrows=5, nl=5)
inBrick <- setValues(b, runif(ncell(b) * nlayers(b)))
inBrick[c(1,2,3,22,23,24,25)] <- NA
ini <- 0.3
p <- 0.15
p1 <- p/3
p2 <- p-(p/3)
fc <- 0.3
weather <- c(0.1, 0, 0, 0, 0.3)

# 将栅格砖转换为数组，方便向量化计算
in_arr <- as.array(inBrick)
dim_in <- dim(in_arr)
out1_arr <- array(NA, dim = dim_in)
out2_arr <- array(NA, dim = dim_in)

# 初始化out2的第一层
out2_arr[,,1][!is.na(in_arr[,,1])] <- ini

# 按层顺序计算（层间有依赖，无法并行）
for (i in 1:dim_in[3]) {
  # 提取当前层的输入和out2数据
  current_in <- in_arr[,,i]
  current_out2 <- out2_arr[,,i]
  
  # 向量化计算out1的当前层
  varr1 <- current_in * ((current_out2 - p1)/p2)^2
  out1_arr[,,i] <- ifelse(!is.na(current_out2),
                          ifelse(current_out2 > p, current_in, varr1),
                          NA)
  
  # 计算out2的下一层（如果不是最后一层）
  if (i < dim_in[3]) {
    var2 <- out2_arr[,,i] + (weather[i] - out1_arr[,,i])/100
    # 向量化取最小值，保留NA
    out2_arr[,,i+1] <- ifelse(!is.na(var2), pmin(fc, var2), NA)
  }
}

# 将数组转换回Raster Brick
outBrick1 <- brick(out1_arr, template = inBrick)
outBrick2 <- brick(out2_arr, template = inBrick)

优化说明 #

1. 数组替代栅格逐元素访问：将Raster Brick转换为数组后，避免了raster包的逐像元访问开销，数组操作更直接高效。
2. 完全移除内层循环：用ifelse、pmin等向量化函数替代逐像元的for循环和条件判断，这些函数底层由编译语言实现，速度提升显著。
3. 仅保留必要的外层循环：由于层与层之间存在依赖关系（out2的第i层依赖out1的i-1层），必须按顺序遍历层数，但181层的循环开销可以忽略。
4. 充分利用内存：数组会直接加载到内存中，你的64GB RAM完全可以容纳当前数据量（约40MB），解决了之前内存利用率低的问题。

测试建议 #

• 先在小数据集上验证优化后代码的计算结果与原始代码一致，确保逻辑正确。
• 用system.time()函数对比优化前后的运行耗时，例如：

system.time({
  # 优化后的计算代码
})

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