基于均值与标准差对RasterBrick/Stack各图层重分类的技术求助

问题背景 #

我有一个名为wnt的RasterBrick/Stack（含8个图层），需要对每个图层按自身的均值（mn）和标准差（sd）重分类，规则如下：

• x > mn + sd → 赋值1
• mn + 0.5sd < x < mn + sd → 赋值2
• mn - 0.5sd < x < mn + 0.5sd → 赋值3
• mn - sd < x < mn - 0.5sd → 赋值4
• x < mn - sd → 赋值5

原尝试代码未得到预期结果：

# create an empty raster brick to store the reclassified data
wnt_reclass <- brick(nrows=nlayers(wnt), ncols=ncell(wnt), 
                     xmn=extent(wnt)[1], xmx=extent(wnt)[2], 
                     ymn=extent(wnt)[3], ymx=extent(wnt)[4], 
                     crs=projection(wnt))

# loop through each layer of the raster brick, calculate the mean and standard deviation 
# for the current layer, and reclassify the data based on the current layer's mean and sd
for (i in 1:nlayers(wnt)) {
  # extract the current layer
  layer <- wnt[[i]]
  
  # calculate the mean and sd for the current layer
  mn <- mean(layer[], na.rm = TRUE)
  sd <- sd(layer[], na.rm = TRUE)
  
  # reclassify the data based on the current layer's mean and sd
  x <- layer[]
  x[x > (mn + sd)] <- 1
  x[(mn + 0.5 * sd) < x & x < (mn + sd)] <- 2
  x[(mn - 0.5 * sd) < x & x < (mn + 0.5 * sd)] <- 3
  x[(mn - sd) < x & x < (mn - 0.5 * sd)] <- 4
  x[x < (mn + sd)] <- 5
  
  # set the reclassified data for the current layer in the output raster brick
  wnt_reclass[[i]] <- setValues(layer, x)
}

# Set the names of the reclassified layers
names(wnt_reclass) <- paste0("wnt_", 1:nlayers(wnt_reclass))

# Set the extent and projection of the reclassified rasterBrick to match the original 'wnt' rasterBrick
extent(wnt_reclass) <- extent(wnt)
projection(wnt_reclass) <- projection(wnt)

错误分析 #

1. 重分类逻辑覆盖问题：最后一行x[x < (mn + sd)] <- 5会覆盖前面所有赋值——除了x > mn + sd的像元，其余所有像元都会被设为5，完全破坏分类规则。
2. 空栅格创建参数错误：nrows=nlayers(wnt)和ncols=ncell(wnt)参数颠倒，导致输出栅格的行列数完全错误，应该使用原栅格的nrow()和ncol()。
3. 变量名冲突：使用sd作为变量名，会与R内置的sd()函数重名，存在潜在风险。

修正后的代码 #

library(raster)

# 快速创建与原栅格结构完全匹配的空输出Stack
wnt_reclass <- stack(wnt)

# 遍历每个图层执行重分类
for (i in 1:nlayers(wnt)) {
  current_layer <- wnt[[i]]
  
  # 计算当前图层的均值和标准差（忽略NA值）
  mn <- mean(current_layer[], na.rm = TRUE)
  sd_val <- sd(current_layer[], na.rm = TRUE)
  
  # 初始化重分类结果为NA，确保未匹配像元保留缺失值
  reclass_vals <- rep(NA, ncell(current_layer))
  
  # 按从极端到中间的顺序赋值，避免区间覆盖问题
  reclass_vals[current_layer[] > mn + sd_val] <- 1
  reclass_vals[current_layer[] > mn + 0.5*sd_val & current_layer[] <= mn + sd_val] <- 2
  reclass_vals[current_layer[] > mn - 0.5*sd_val & current_layer[] <= mn + 0.5*sd_val] <- 3
  reclass_vals[current_layer[] > mn - sd_val & current_layer[] <= mn - 0.5*sd_val] <- 4
  reclass_vals[current_layer[] <= mn - sd_val] <- 5
  
  # 将重分类结果赋值给输出图层
  wnt_reclass[[i]] <- setValues(current_layer, reclass_vals)
}

# 设置输出图层名称
names(wnt_reclass) <- paste0("wnt_", 1:nlayers(wnt_reclass))

关键优化点 #

• 用stack(wnt)直接创建结构匹配的空栈，避免手动输入参数出错
• 调整赋值顺序为从极端区间到中间区间，同时用<=明确边界值归属（离散栅格像元用闭区间更严谨）
• 变量名sd改为sd_val，避免与内置函数冲突
• 初始化结果为NA，确保异常值或未匹配像元保留缺失状态

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