R terra处理多类别栅格如何正确提取分类属性为独立图层且不丢失数据

安全阿岳

2026-10-2

问题描述

我在使用R的Terra包处理LANDFIRE发布的CONUS历史扰动数据集（下载地址：https://landfire.gov/version_download.php，数据集名称为HDist）时遇到问题。我的需求是将该数据集裁剪、投影到目标范围后，把单元格存储的不同属性值拆分为独立图层，比如扰动烈度层、扰动类型层等。该历史扰动数据集的所有属性都存储在同一个属性表中，在terra中这些属性表对应categories设置，我在裁剪和投影步骤都运行正常，但无法正确提取属性值并拆分为独立图层，现有代码如下：

library(terra)

setwd("your pathway to historical disturbance tif here")

h1 <- terra::rast("LC16_HDst_200.tif")  # 读入Hdist tif文件
h2 <- terra::project(h1, "EPSG:5070", method = "near") # 使用最邻近法投影
h3 <- crop(h2, ext(xmin,xmax,ymin,ymax)) # 按指定范围裁剪
h3

运行上述代码后得到了符合目标范围和投影的栅格对象，其分类属性如下：

categories  : Count, HDIST_ID, DISTCODE_V, DIST_TYPE, TYPE_CONFI, SEVERITY, SEV_CONFID, HDIST_CAT, FDIST, R, G, B

我了解这类数据集的取值存储在上述分类中，但使用plot(h3)绘图时仅显示Count分类的第一行值。我使用如下代码切换激活分类：

activeCat(h3) <- 4
h3

运行后输出如下，已将激活分类切换为第4项DIST_TYPE：

name       :   DIST_TYPE
min value  :   Clearcut
max value  : Wildland Fire Use

默认激活分类为count，切换逻辑无问题，但再次运行plot(h3)绘图时仅显示NoData，无其他有效内容。我尝试使用catalyze()函数将所有分类转换为数值图层：

h4 <- catalyze(h3)

运行后得到13个图层，对应13个分类，符合预期，但尝试绘制第4层（对应DIST_TYPE分类）时：

plot(h4, 4) # 绘制h4的第4层，对应DIST_TYPE分类

仅显示取值为8的内容，看起来只有NoData值，图层最小值和最大值均为8，其他分类图层也存在类似问题，似乎每个图层仅使用了属性表第一行的值，直接访问值还会导致程序崩溃。
综上，terra显然存储了arcgis中能正常查看的全部属性表内容，但无论绘图还是数据操作都仅能访问属性表的顶行，使用catalyze()转换后数据反而更加混乱。我清楚在arcgis pro中可以轻松实现需求，但为了文档连贯性希望全流程在R中完成，请问该如何解决该问题？我使用LANDFIRE evt数据时也遇到相同问题，dem、冠层覆盖等简单栅格无此问题，仅存在多分类（属性表多列）的栅格会出现该问题。
补充：当前即使尝试修复后绘图结果仍如下图所示：
修复后仍异常的绘图结果

解决方案

该问题是LANDFIRE类带属性表（RAT）的栅格在terra中处理的常见问题，根源是project、crop操作后栅格原始值与属性表的映射关系出错，可按以下方法解决：

优先在原始未投影的栅格上先完成属性拆分，再做投影、裁剪操作，能最大程度避免映射错误
若必须先做投影裁剪，可手动提取属性表做值匹配，规避activeCat和catalyze的兼容性bug，代码如下：

library(terra)

# 读入原始栅格
h1 <- rast("LC16_HDst_200.tif")

# 先提取完整属性表
rat <- levels(h1)[[1]]

# 执行原有投影、裁剪逻辑
h2 <- project(h1, "EPSG:5070", method = "near")
h3 <- crop(h2, ext(xmin,xmax,ymin,ymax))

# 提取栅格的原始编码值（不是分类标签）
h3_values <- values(h3, mat = FALSE)

# 新建空栅格堆栈存储拆分后的图层
output_stack <- rast()

# 遍历属性表的属性列（按需调整列范围）
for(col in colnames(rat)[-1]){
  # 新建单图层，将原始编码匹配为对应属性值
  temp_rast <- h3
  values(temp_rast) <- rat[[col]][match(h3_values, rat[[1]])]
  # 设置图层名
  names(temp_rast) <- col
  # 加入输出堆栈
  output_stack <- c(output_stack, temp_rast)
}

# 验证输出：绘制DIST_TYPE图层
plot(output_stack$DIST_TYPE)

若需要保留分类属性而非转为普通数值图层，可在匹配完成后重新给对应图层设置levels：

# 以DIST_TYPE为例
dist_type_levels <- unique(rat[,c("Value", "DIST_TYPE")])
levels(output_stack$DIST_TYPE) <- dist_type_levels

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