如何用R语言terra包将多波段土地覆盖栅格合并为单图层

问题描述 #

我用R语言合并多张土地覆盖栅格图，代码如下：

tif_files <- list.files(path = "LandCover", 
                        pattern = ".tif$", 
                        full.names = TRUE, 
                        ignore.case = TRUE)  
tif_rasts <- list()  
for (i in 1:length(tif_files)){  
  tif_rasts <- c(tif_rasts,  rast(tif_files[i]))  
}  

tif_mosaic <- do.call(mosaic, c(tif_rasts, fun="max"))  
plot(tif_mosaic)  

writeRaster(tif_mosaic,  
            filename = "land_cover.tif"),  
            overwrite = TRUE  
)  

但最终得到的土地覆盖数据包含3个波段，各波段值有相同也有不同，我需要把这些波段合并成单个图层，方便对应土地利用类型。运行unique(tif_mosaic)得到以下输出：

unique(tif_mosaic)  
   CHI_1_LULC_1 CHI_1_LULC_2 CHI_1_LULC_3
2             0            0            0
3             0            0          255
4            34          139           34
5           128            0            0
6           128          236          104
7           133          133          133
8           255            0            0
9           255            0          192
10          255          193           37
11          255          255          255

解决方案 #

问题根源 #

你得到多波段结果的原因是输入的原始tif文件本身是3波段（大概率是RGB格式的土地覆盖分类图），rast()读取后保留了多波段结构，后续的mosaic操作会按波段分别合并，最终输出3波段栅格。

方法1：将现有3波段栅格转成单波段唯一编码 #

把三个波段的数值组合成唯一字符串，再转成数值型ID，每个ID对应一组唯一的RGB值，后续可以映射到土地利用类型：

library(terra)

# 合并三个波段为唯一字符串编码
combined_str <- paste0(tif_mosaic[[1]], "_", tif_mosaic[[2]], "_", tif_mosaic[[3]])
# 转成因子型栅格，再转成数值型ID
single_band <- rast(as.factor(combined_str))

# 查看每个ID对应的原始RGB值
unique_mapping <- unique(cbind(values(tif_mosaic), land_id = values(single_band)))
print(unique_mapping)

方法2：根据RGB值直接重分类为土地利用类型 #

如果知道每个RGB组合对应的土地利用类型，可以创建映射表，用classify()直接生成带类型标签的单波段栅格：

library(terra)

# 示例映射表，根据你的实际土地分类规则修改
lut <- data.frame(
  R = c(0, 0, 34, 128, 128, 133, 255, 255, 255, 255),
  G = c(0, 0, 139, 0, 236, 133, 0, 0, 193, 255),
  B = c(0, 255, 34, 0, 104, 133, 0, 192, 37, 255),
  land_type = c("未分类", "水体", "林地", "耕地", "草地", "裸地", "建设用地", "湿地", "园地", "其他")
)

# 创建重分类矩阵：每行是R,G,B对应目标ID
reclass_mat <- cbind(lut[,1:3], 1:nrow(lut))
# 执行重分类
single_band <- classify(tif_mosaic, reclass_mat, right = TRUE)
# 给栅格添加土地类型属性
levels(single_band) <- data.frame(ID = 1:nrow(lut), land_type = lut$land_type)

# 查看并保存结果
plot(single_band)
writeRaster(single_band, filename = "single_band_landcover.tif", overwrite = TRUE)

优化原始合并流程 #

如果后续还要合并同类型的多波段tif，建议先将每个输入文件转成单波段编码再合并，避免生成多波段结果：

library(terra)

tif_files <- list.files(path = "LandCover", 
                        pattern = ".tif$", 
                        full.names = TRUE, 
                        ignore.case = TRUE)  

# 遍历文件，将每个多波段栅格转成单波段ID
tif_rasts <- lapply(tif_files, function(file) {
  r <- rast(file)
  combined_str <- paste0(r[[1]], "_", r[[2]], "_", r[[3]])
  rast(as.factor(combined_str))
})

# 合并单波段栅格
tif_mosaic <- do.call(mosaic, c(tif_rasts, fun = "max"))
plot(tif_mosaic)
writeRaster(tif_mosaic, filename = "land_cover_single.tif", overwrite = TRUE)

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