I need to clip the national landcover database raster to a watershed boundary vector. This post works for some of the watersheds I am working with, but not all of them. For the ones that it does not work with, the issue appears to be that the mask step results in an entirely empty raster. For sake of brevity of this post, I am not including an example of two watersheds where this works and does not work. The example below is for a watershed that did not work in the workflow in the post mentioned above.
My solution for the watersheds that the mask step does not work for is to instead reproject the raster to the vector's CRS, then perform the clip and mask. To do so, I need to use terra::rast to get a SpatRaster so I can use terra::project. While this works to clip and mask the raster, a new issue has come up, where terra::rast is changing the cell values of the raster from numeric (i.e. the NLCD ID, which I want) to factor (i.e. the NLCD Class), which would be okay BUT the strings are different...e.g. the NLCD legend has ID 81 as "Pasture/Hay" but this gets changed to "Hay/Pasture" when using rast. Thus, this prevents me from using left_join correctly.
My end goal is to reclassify (i.e. simplify to 8 classes) the NLCD using a custom legend and then calculate the percent land for each reclassified land use.
Here is my reproducible example:
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))
What I am trying to understand is 1) why terra::rast actually changes the rasters cell values in the first place:
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
and 2) why are the string values changed from the NLCD legend (which comes from the function pal_nlcd() used above)? If the class names were the same after using terra::rast as the NLCD legend, the left_join would work by 'Class', but this is not the case since I am missing some classes in my final result, nlcd_df. Here you can see how Pasture/Hay in the NLCD legend is switched to Hay/Pasture when using terra::rast:
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"
Also if you are having trouble installing streamstats you may need to try these steps first:
# 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)
