使用terra包为分类栅格添加属性数据出现异常结果的求助

问题描述 #

将栅格数据分为3类并为每类添加对应数值，使用terra包执行以下代码后，生成的第二个栅格图像显示异常（如图所示），需修正代码：

library(terra)

f <- system.file("ex/elev.tif", package="terra") 
r <- rast(f)

# classify the values into three groups
m <- c(min(global(r, "min", na.rm=TRUE)), 275, 1,
       275, 410, 2,
       410, max(global(r, "max", na.rm=TRUE)), 3)

#Reclassify the raster stack
rclmat <- matrix(m, ncol=3, byrow=TRUE)
rc <- classify(r, rclmat, include.lowest=TRUE)
plot(rc)

#Add some values to raster attribute table
rc$code <- c(0.002, 0.0056, 0.0124)

#Generate raster using the new column
y <- as.numeric(rc, 2)
plot(y)

问题原因 #

as.numeric(rc, 2)的用法不符合terra包的设计逻辑：该函数的第二个参数是指定提取的波段索引，但rc是单波段栅格，因此返回全NA的结果，导致图像异常。它并不能实现将属性表中的自定义列映射为栅格值的需求。

修正方案 #

使用terra::subs()函数，将重分类后的栅格类别值（1、2、3）替换为属性表中对应的code值，这是terra包中实现类别值映射的标准方法。

修正后的完整代码 #

library(terra)

f <- system.file("ex/elev.tif", package="terra") 
r <- rast(f)

# 分类为三组
m <- c(min(global(r, "min", na.rm=TRUE)), 275, 1,
       275, 410, 2,
       410, max(global(r, "max", na.rm=TRUE)), 3)
rclmat <- matrix(m, ncol=3, byrow=TRUE)
rc <- classify(r, rclmat, include.lowest=TRUE)
plot(rc)

# 添加自定义属性列
rc$code <- c(0.002, 0.0056, 0.0124)

# 利用属性表映射生成新栅格
y <- subs(rc, rc, which="code")
plot(y)

代码说明 #

• subs(rc, rc, which="code")：第一个rc是待替换的源栅格，第二个rc提供属性表作为查找依据，which="code"指定用属性表中的code列值替换原栅格的类别值。
• 也可以手动构建查找表来实现替换，适合更复杂的映射场景：

  lookup <- data.frame(ID = c(1,2,3), code = rc$code)
  y <- subs(rc, lookup, by="ID", which="code")
  

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