栅格属性的快速重命名与合并方法探究

问题背景 #

使用R语言terra包处理栅格数据时，通过属性表将栅格ID映射为不同层级的类别（如category1仅重命名、category2/category3合并部分ID），但用terra::subst()处理大栅格（如400MB）时耗时极长，希望找到无需逐像素修改的快速方案，同时保留绘图与后处理能力。

现有可行的重命名方案：通过levels()关联属性表，绘图时直接显示类别标签而非ID：

levels(example_raster) <- example_attribute_table 
plot(example_raster)

快速解决方案 #

1. 基于levels的多类别可视化（无需修改栅格） #

利用terra的栅格属性表（RAT）存储所有层级类别，绘图时手动指定映射规则，实现不同类别层级的可视化，完全不修改原始栅格数据：

# 将完整属性表关联到栅格的levels
levels(example_raster) <- example_attribute_table

# 按category2合并类别绘图
# 获取category2的唯一值与对应ID映射
cat2_unique <- unique(example_attribute_table[, c("id", "category2")])
# 绘图时指定颜色、图例标签与断点
plot(example_raster,
     col = terra::hcl.colors(length(unique(cat2_unique$category2))),
     legend.values = cat2_unique$category2,
     breaks = cat2_unique$id)

# 按category3合并类别同理
cat3_unique <- unique(example_attribute_table[, c("id", "category3")])
plot(example_raster,
     col = terra::hcl.colors(length(unique(cat3_unique$category3))),
     legend.values = cat3_unique$category3,
     breaks = cat3_unique$id)

2. 转换为分类栅格（Categorical Raster） #

将栅格转为因子类型后，可直接基于levels中的多列属性进行层级分析与可视化，操作更简洁：

# 转换为分类栅格（因子类型）
example_raster <- as.factor(example_raster)
levels(example_raster) <- example_attribute_table

# 直接调用属性列进行绘图
# 提取category3的唯一标签
cat3_labels <- unique(levels(example_raster)[[1]]$category3)
plot(example_raster,
     col = terra::hcl.colors(length(cat3_labels)),
     legend.values = cat3_labels,
     breaks = match(cat3_labels, levels(example_raster)[[1]]$category3))

3. 用Geopackage存储分类栅格与属性 #

Geopackage原生支持存储栅格的完整属性表（RAT），可一次性将栅格与多类别属性存入文件，后续读取直接调用，避免重复关联属性的操作：

# 将分类栅格写入Geopackage，保留所有属性列
writeRaster(example_raster, "classified_soil.gpkg", overwrite = TRUE, datatype = "INT1U")

# 读取时自动加载属性表
gpkg_raster <- terra::rast("classified_soil.gpkg")

# 直接按category2绘图
plot(gpkg_raster,
     legend.values = levels(gpkg_raster)[[1]]$category2)

4. 向量级快速生成合并类别的栅格（若需实体栅格） #

如果后续分析必须使用合并后的栅格数据，可通过向量匹配替代subst()，大幅提升效率：

# 基于属性表生成category2的映射向量
cat2_mapping <- example_attribute_table$category2[match(values(example_raster), example_attribute_table$id)]

# 创建新栅格并赋值（用因子类型节省存储空间）
raster_cat2_fast <- example_raster
values(raster_cat2_fast) <- factor(cat2_mapping, levels = unique(cat2_mapping))

# 查看结果
plot(raster_cat2_fast)

方案优势总结 #

• 前三种方案均不修改原始栅格像素值，仅通过属性映射实现类别可视化，处理大栅格几乎瞬时完成；
• 第四种方案通过向量级操作替代逐像素替换，效率远高于subst()或reclassify()；
• Geopackage存储方案适合长期复用，避免重复处理属性关联。

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