我需要将国家土地覆盖数据库栅格裁剪为流域边界矢量。This post适用于我正在使用的一些流域，但不是所有流域。对于它不起作用的那些，问题似乎是mask步骤导致完全空的光栅。为了这篇文章的简洁起见，我没有包括两个流域的例子，这两个流域是有效的，也是无效的。下面的示例是针对在上述帖子中的工作流中不起作用的分水岭。
对于mask步骤不适用的流域，我的解决方案是将栅格重新投影到矢量的CRS，然后执行裁剪和遮罩。为此，我需要使用terra::rast来获得SpatRaster，以便可以使用terra::project。虽然这可以裁剪和遮罩栅格，但出现了一个新问题，其中terra::rast将栅格的像元值从数字（即我想要的NLCD ID）进行分解（即NLCD类），这是可以的，但是字符串是不同的。NLCD图例的ID 81为“Pasture/Hay”，但在使用rast时会更改为“Hay/Pasture”。因此，这阻止了我正确使用left_join。
我的最终目标是重新分类（即）。使用自定义图例将NLCD简化为8类），然后计算每个重新分类的土地使用的土地百分比。
以下是我的可重复示例：

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

# create the custom reclassified legend (needed below):

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')

# download the watershed boundary and convert to sf:

setTimeout(400)
ws1 <- delineateWatershed(xlocation = -73.65167, ylocation = 42.93556, crs = 4326,includeparameters = "false") # takes a minute to run
ws1_sp<-toSp(ws1, what = 'boundary') # also takes a minute to run
ws1_sf<-st_as_sf(ws1_sp)

# download the NLCD raster data:

nlcd<-get_nlcd(template = ws1_sf,label = 'HUDSON RIVER AT STILLWATER NY',year = 2016)

# convert to SpatRaster, reproject to vector's crs, then clip and mask:

nlcd_rast<-rast(nlcd)
nlcd_rast<-project(nlcd_rast, crs(ws1_sf))
nlcd_rast<-crop(nlcd_rast, raster::extent(ws1_sf))
nlcd_rast<-mask(nlcd_rast, ws1_sf)

# calculate the percent land of the reclassified land uses:

nlcd_df<-as.data.frame(nlcd_rast, xy = FALSE)%>%
  drop_na()%>%
  rename(ID = NLCD.Land.Cover.Class)%>%
  left_join(.,legend[,c(1,2,5,6)], by = c('ID'='Class'))%>%
  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))

我试图理解的是1）为什么terra::rast实际上首先改变了栅格单元值：

class(as.data.frame(nlcd, xy = FALSE)$NLCD.Land.Cover.Class) # numeric
class(as.data.frame(nlcd_rast, xy = FALSE)$NLCD.Land.Cover.Class) # factor

和2）为什么字符串值从NLCD图例（来自上面使用的函数pal_nlcd()）更改？如果在使用terra::rast作为NLCD图例后，类名相同，则left_join将按“Class”工作，但情况并非如此，因为我在最终结果nlcd_df中缺少一些类。在这里，您可以看到使用terra::rast时，NLCD图例中的Pasture/Hay如何切换为Hay/Pasture：

unique(as.character(as.data.frame(nlcd_rast, xy = FALSE)$NLCD.Land.Cover.Class))
 [1] "Evergreen Forest"             "Herbaceous"                  
 [3] "Open Water"                   "Woody Wetlands"              
 [5] "Emergent Herbaceous Wetlands" "Shrub/Scrub"                 
 [7] "Deciduous Forest"             "Barren Land"                 
 [9] "Mixed Forest"                 "Developed, Open Space"       
[11] "Developed, Low Intensity"     "Developed, Medium Intensity" 
[13] "Developed, High Intensity"    "Hay/Pasture"                 
[15] "Cultivated Crops"            
> legend$Class
 [1] "Open Water"                   "Perennial Ice/Snow"          
 [3] "Developed, Open Space"        "Developed, Low Intensity"    
 [5] "Developed, Medium Intensity"  "Developed High Intensity"    
 [7] "Barren Land (Rock/Sand/Clay)" "Deciduous Forest"            
 [9] "Evergreen Forest"             "Mixed Forest"                
[11] "Dwarf Scrub"                  "Shrub/Scrub"                 
[13] "Grassland/Herbaceous"         "Sedge/Herbaceous"            
[15] "Lichens"                      "Moss"                        
[17] "Pasture/Hay"                  "Cultivated Crops"            
[19] "Woody Wetlands"               "Emergent Herbaceous Wetlands"

如果您在安装streamstats时遇到问题，您可能需要先尝试以下步骤：

# step 1 - restart R

# Session->Restart R

# step 2 - add Rtools to PATH

old_path <- Sys.getenv("PATH")

Sys.setenv(PATH = paste(old_path, "C:\\rtools40\\usr\\bin", sep = ";"))

# now try:

devtools::install_github("markwh/streamstats")

library(streamstats)