Terra栅格含空值时条件赋值报错问题及解决问询

我有三个栅格图层：lc（土地覆被）、b（输入积雪覆盖）、c（修改后积雪覆盖）。需要实现：当lc取值为"anthropogenic"时，用b的数值替换c对应位置的数值。无空值时操作正常，但栅格包含空值（如研究区边缘）时，即使各栅格范围、分辨率完全匹配，仍会抛出错误：

[`[<-`] length of cells and values do not match

复现问题的测试代码 #

library(terra)
library(geodata)

# 导入世界国家数据集
w <- world(path=".", resolution = 1)
# 选择三个沿海国家
sel <- w[w$NAME_0 %in% c('Ghana', 'Togo', 'Benin'), ]
# 添加数值型虚拟变量
sel$num <- c(1, 2, 3)

# 导入土壤数据集
soil_full <- soil_af(var="pH", depth=15, path=".")
# 裁剪到三个国家范围
soil <- crop(soil_full, ext(sel))
soil_dummy <- soil
# 将国家矢量数据栅格化为与土壤栅格一致的范围和分辨率
values(soil_dummy) <- NA
sel_rast <- rasterize(x=sel, y=soil_dummy, field="num",
                     background=NA, update=TRUE, touches=TRUE)
names(sel_rast) <- "num"
par(mfrow=c(1,2))
plot(sel_rast)
plot(soil)

# 尝试赋值，触发错误
sel_rast[sel_rast==1] <- soil[sel_rast==1]
# Error: [`[<-`] length of cells and values do not match

# 查看栅格信息
sel_rast
# class       : SpatRaster 
# dimensions  : 921, 853, 1  (nrow, ncol, nlyr)
# resolution  : 0.008333333, 0.008333333  (x, y)
# extent      : -3.258333, 3.85, 4.741667, 12.41667  (xmin, xmax, ymin, ymax)
# coord. ref. : lon/lat WGS 84 (EPSG:4326) 
# source(s)   : memory
# name        : num 
# min value   :   1 
# max value   :   3

soil
# class       : SpatRaster 
# dimensions  : 921, 853, 1  (nrow, ncol, nlyr)
# resolution  : 0.008333333, 0.008333333  (x, y)
# extent      : -3.258333, 3.85, 4.741667, 12.41667  (xmin, xmax, ymin, ymax)
# coord. ref. : lon/lat WGS 84 (EPSG:4326) 
# source(s)   : memory
# name        : pH_5-15cm 
# min value   :      4.68 
# max value   :      7.85 

# 尝试用mask处理，仍失败
test <- mask(soil, sel_rast)
sel_rast[sel_rast==1] <- test[sel_rast==1] # 失败

# 裁剪到无空值的内陆区域，赋值成功
sel_rast.c <- crop(sel_rast, ext(sel_rast)/10)
soil.c <- crop(soil, ext(sel_rast)/10)
plot(soil.c)
plot(sel_rast.c)
sel_rast.c[sel_rast.c==3] <- soil.c[sel_rast.c==3] # 成功
plot(sel_rast.c)

问题原因 #

直接使用x[condition] <- y[condition]赋值时，若condition逻辑向量中包含NA（对应栅格空值区域的判断结果），会导致左右两边提取的向量长度或有效值数量不匹配——空值位置的逻辑判断结果为NA，这类位置在提取时会被特殊处理，最终引发长度不匹配错误。

解决方法 #

方法1：使用terra::ifel函数（推荐）

ifel是terra包专为栅格设计的条件赋值函数，可自动处理空值，语法简洁且高效：

# 对应需求：当lc取值为"anthropogenic"时，c替换为b的值，否则保留c原数值
c <- ifel(lc == "anthropogenic", b, c)

方法2：通过索引精准匹配

先提取符合条件的非空单元格索引，再进行赋值，确保左右两边长度完全一致：

# 获取lc等于目标值的非空单元格索引
idx <- which(lc[] == "anthropogenic")
# 执行赋值
c[idx] <- b[idx]

针对测试代码中的场景，修改为：

# 获取sel_rast等于1的非空单元格索引
idx <- which(sel_rast[] == 1)
# 赋值
sel_rast[idx] <- soil[idx]

内容的提问来源于stack exchange，提问作者kl-higgins