为何terra::rast转换NLCD栅格时将数值型转为因子且类名变更？

问题描述 #

我需要将国家土地覆盖数据库（NLCD）栅格裁剪至流域边界矢量。部分流域使用mask步骤会生成空栅格，因此改用terra::rast将栅格转为SpatRaster，投影至矢量CRS后执行裁剪与掩膜操作，但遇到两个问题：

• terra::rast将原数值型的NLCD ID转为因子型的NLCD Class；
• 因子类名与pal_nlcd()生成的图例不一致（如原“Pasture/Hay”变为“Hay/Pasture”），导致left_join无法正常匹配，影响后续自定义重分类的土地利用占比计算。

解决方案 #

1. 核心思路：基于数值ID匹配，避免因子转换问题 #

terra::rast读取带属性表的栅格时默认会转成因子，且类名与pal_nlcd()的标准命名存在差异。最可靠的解决方式是保留原始数值ID，直接用ID做匹配，彻底规避类名字符串不一致的问题。

修改后的完整代码 #

library(streamstats) # devtools::install_github("markwh/streamstats")
library(FedData)
library(terra)
library(sf)
library(tidyverse)

# 创建自定义重分类图例（基于NLCD数值ID匹配）
legend <- pal_nlcd() %>%
  mutate(ID2 = c(1,2,3,3,3,3,2,4,4,4,4,4,5,5,2,2,6,7,8,8))
legend_2 <- data.frame(ID2 = unique(legend$ID2)) %>% 
  mutate(Class2 = c("Open Water", "Other", "Developed", "Forest", "Grassland", "Pasture/Hay", "Cultivated Crops", "Wetlands"))
legend <- left_join(legend, legend_2, by = 'ID2')

# 下载流域边界并转为sf格式
setTimeout(400)
ws1 <- delineateWatershed(xlocation = -73.65167, ylocation = 42.93556, crs = 4326, includeparameters = "false")
ws1_sp <- toSp(ws1, what = 'boundary')
ws1_sf <- st_as_sf(ws1_sp)

# 下载NLCD栅格数据
nlcd <- get_nlcd(template = ws1_sf, label = 'HUDSON RIVER AT STILLWATER NY', year = 2016)

# 转为SpatRaster时强制保留数值ID，不转因子
nlcd_rast <- rast(nlcd, factors = FALSE)
# 投影、裁剪、掩膜（使用terra原生函数替代raster包函数）
nlcd_rast <- project(nlcd_rast, crs(ws1_sf))
nlcd_rast <- crop(nlcd_rast, ext(ws1_sf))
nlcd_rast <- mask(nlcd_rast, ws1_sf)

# 计算重分类后的土地利用占比（基于数值ID匹配）
nlcd_df <- as.data.frame(nlcd_rast, xy = FALSE) %>%
  drop_na() %>%
  rename(ID = NLCD.Land.Cover.Class) %>%
  # 用数值ID匹配，彻底解决类名不一致问题
  left_join(., legend[,c(1,2,5,6)], by = c('ID'='ID')) %>%
  group_by(Class2) %>%
  summarise(n = n()) %>%
  mutate(pland = round(n / sum(n), 4)) %>%
  arrange(Class2) %>%
  dplyr::select(c(1,3)) %>%
  complete(., Class2 = unique(legend$Class2), fill = list(pland = NA))

关键修改点 #

• rast(nlcd, factors=FALSE)：强制保留原始数值型ID，避免自动转为因子；
• 匹配逻辑从类名字符串改为数值ID：彻底规避类名表述差异的问题；
• 替换raster::extent为terra::ext：保持terra包生态的一致性，减少跨包调用风险。

验证修改效果 #

# 验证ID仍为数值型
class(as.data.frame(nlcd_rast, xy = FALSE)$NLCD.Land.Cover.Class) # 输出应为numeric
# 查看匹配后的占比结果
head(nlcd_df)

streamstats包安装问题解决 #

若安装streamstats遇到报错，可按以下步骤操作：

# 步骤1 - 重启R
# Session->Restart R

# 步骤2 - 将Rtools添加到系统PATH
old_path <- Sys.getenv("PATH")
Sys.setenv(PATH = paste(old_path, "C:\\rtools40\\usr\\bin", sep = ";"))

# 重新安装包
devtools::install_github("markwh/streamstats")
library(streamstats)

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