如何用Terra分类栅格掩膜？大栅格高效掩膜方法咨询

问题背景 #

我有一个超大栅格rasterA，值范围为1:928863，需要将另一栅格rasterB中对应rasterA值>300的区域设为NA。栅格尺寸模拟代码如下：

library(terra)
set.seed(0)
rasterA <- rast(nrows=42893, ncols=52031)
rasterB <- rasterA
values(rasterA) <- sample(1:928563, ncell(rasterA), replace=TRUE)
values(rasterB) <- 1

当前采用的方法运行数小时仍未完成：

rasterA <- subst(rasterA, 300:928863, NA)
rasterB  <- mask(rasterB, rasterA)

尝试将rasterA转为分类栅格后掩膜，但未达到预期效果：

leveller <- data.frame(id = 1:928863, dist= c(rep("good", 300), rep(NA, 928563)))
levels(rasterA) <- leveller
rasterB <- mask(rasterB, rasterA)

以下小示例可复现问题：希望在分类栅格rf对应NA标签的区域掩膜栅格r，但运行后r仍为全值1：

library(terra)
set.seed(0)
rf <- rast(nrows=10, ncols=10)
values(rf) <- sample(3, ncell(rf), replace=TRUE)
cls <- data.frame(id=1:3, cover=c("forest", "water", NA))
levels(rf) <- cls

r <- rast(nrows=10, ncols=10)
values(r) <- 1

r <- mask(r, rf)

需求：能否通过分类栅格的因子水平进行掩膜？或有无更高效的方法实现最初的掩膜需求？

解决方案 #

1. 最优高效方法：直接使用逻辑运算生成掩膜 #

核心需求是保留rasterA值≤300的区域，其余设为NA，无需遍历大量值进行替换，直接通过向量化逻辑运算生成掩膜，效率远超subst：

针对超大栅格的代码：

# 方法1：用乘法直接赋值
mask_layer <- rasterA <= 300
rasterB <- rasterB * mask_layer

# 方法2：结合mask函数更直观
rasterB <- mask(rasterB, rasterA <= 300)

这种方法无需修改原始rasterA，逻辑运算为底层优化的向量操作，处理超大栅格速度会大幅提升。

2. 分类栅格掩膜失效原因及修正 #

terra::mask默认基于栅格单元格的实际NA值进行掩膜，而你设置的NA是分类属性表中的标签，并非单元格的实际值。也就是说，rf的单元格值仍为1、2、3，属性表的NA标签不会改变单元格值，因此mask不生效。

若一定要用分类栅格实现，需先将对应标签的单元格转为实际NA：

# 小示例修正代码
library(terra)
set.seed(0)
rf <- rast(nrows=10, ncols=10)
values(rf) <- sample(3, ncell(rf), replace=TRUE)
cls <- data.frame(id=1:3, cover=c("forest", "water", NA))
levels(rf) <- cls

# 将属性表中对应NA标签的单元格转为实际NA
rf[rf == which(is.na(cls$cover))] <- NA

r <- rast(nrows=10, ncols=10)
values(r) <- 1

r <- mask(r, rf)

但该方法对超大栅格的效率不如直接逻辑运算，不推荐。

3. 超大栅格内存优化建议 #

若栅格为磁盘存储的大文件，可启用分块处理避免内存溢出：

options(terra.rasterize.chunksize = 1e7) # 设置分块大小，根据内存调整
rasterB <- mask(rasterB, rasterA <= 300)

内容的提问来源于stack exchange，提问作者we need a Mat. Stat.