如何在R语言中读取Landsat位值质量标志栅格？

问题背景 #

有一个8bit TIFF格式的Landsat QA栅格，标志按位编码，使用terra库提取位信息时出现错误。原代码及报错如下：

library(terra)

toBinary <- function(i){
  valuesbin = paste0(as.integer(rev(intToBits(i)[1:8])), collapse = "")
  valuesbin = unlist(strsplit(valuesbin, split = ""))
  as.numeric(valuesbin)
  
}

qa = rast('./qarast.tiff')
qa = app(qa,toBinary)

报错信息：

Error: [app] the number of values returned by 'fun' is not appropriate

解决方法 #

1. 修复terra代码实现位提取 #

报错根源是app()要求每个输入值返回单个结果，但你的函数返回了8个值（对应8位的每一位）。改用lapp()处理，可生成包含每一位的栅格栈：

library(terra)

# 定义提取指定位置位的函数（位位置从右往左数，1为最低位）
extract_single_bit <- function(x, bit_position) {
  # 位运算：将1左移对应位置，与原数值做按位与，判断是否非0
  as.integer(bitwAnd(x, bitwShiftL(1, bit_position - 1)) != 0)
}

# 读取QA栅格
qa_rast <- rast('./qarast.tiff')

# 提取8位中每一位，生成多图层栅格
qa_bit_stack <- lapply(1:8, function(pos) {
  lapp(qa_rast, extract_single_bit, bit_position = pos)
}) %>% rast()

# 给各图层命名，方便识别
names(qa_bit_stack) <- paste0("bit_", 1:8)

2. 使用专用库简化操作 #

推荐使用专门处理Landsat数据的库，无需手动编写位运算代码，更不易出错：

RStoolbox库

内置decodeQA()函数可直接解析Landsat QA标志：

library(RStoolbox)
library(terra)

qa_rast <- rast('./qarast.tiff')
# 解析QA栅格，返回包含各类质量标志的栅格栈
qa_flags <- decodeQA(qa_rast, type = "Landsat8")

# 查看生成的标志图层，比如云、云阴影、水体等
print(names(qa_flags))

Landsat8库

针对Landsat 8数据做了优化，可直接读取并解析QA文件：

library(Landsat8)

# 读取QA文件并解析
qa_parsed <- read_QA("qarast.tiff")

# 获取特定标志的栅格，例如云覆盖掩码
cloud_mask <- qa_parsed$cloud

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