Calculating the area of each county that has a specific land cover value (R or QGIS)

Question

I need to calculate the area of each U.S. county that is covered by a specific type of land cover.

I am using a county shapefile (polygons) and the land cover data is a raster that comes from the NLCD CONUS 2011 (https://www.mrlc.gov/data/nlcd-2011-land-cover-conus-0).

I realize that there are similar questions to this on SE, however I haven't had any luck with them I think due to the size of the raster.

I have been at this for several days and tried various methods in both R (which I'm very proficient at but haven't used before to manipulate raster data) and QGIS (which I have only used a couple times).

I have been unable to convert the raster to a shapefile due to it's size. After hours of running it doesn't even get to 1% completion. I've tried the R packages stars, raster, and also GDAL through terminal as well as using polygonize from QGIS. I only need one of the land coverage types (value = 82, cultivated crops) and have tried eliminating the other values to reduce the number of polygons that would be created, but it has not helped. Here is the info on the raster:

class      : RasterLayer 
dimensions : 104424, 161190, 16832104560  (nrow, ncol, ncell)
resolution : 30, 30  (x, y)
extent     : -2493045, 2342655, 177285, 3310005  (xmin, xmax, ymin, ymax)
crs        : +proj=aea +lat_0=23 +lon_0=-96 +lat_1=29.5 +lat_2=45.5 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
source     : /[hidden...]/NLCD_2011_Land_Cover_L48_20190424/NLCD_2011_Land_Cover_L48_20190424.img 
names      : NLCD_2011_Land_Cover_L48_20190424 
values     : 0, 95  (min, max)
attributes :
        ID      COUNT Red Green Blue NLCD.Land.Cover.Class Opacity
 from:   0 7854240512   0     0    0          Unclassified       0
  to : 255          0   0     0    0                             0

I have also tried to convert the shapefile of counties into a raster object using both R and QGIS, but haven't had any luck getting a version that will combine with the landcover raster and allow me to look at the overlap that way.

I'd be grateful for any help in order to get to the final goal which is the area of each county covered by the specific land type. Either R or QGIS (or any other free tool) would be great!

---

[EDIT]

Here is some of the code I've tried so that people can actually help answer my question. I did not include the code because this is an issue of not understanding the process I should undertake using either QGIS, R, or both, not a code error or coding problem. Comments are mostly not mine, but instead from code I cribbed

library(tidyverse)
library(siverse)
library(conflicted)
conflict_prefer("select", "dplyr")
conflict_prefer("filter", "dplyr")
library(tigris)
library(sf)
library(albersusa)
library(raster)
library(rgdal)

nl <- raster("NLCD_2011_Land_Cover_L48_20190424.img")


nl@crs


nl %>% st_contour()

levels(nl)[[1]] %>% as_tibble %>% filter(NLCD.Land.Cover.Class == "Cultivated Crops")

n_values <- ncell(nl) * nlayers(nl)

keep82 <- function(x) {x[x!=82] <- NA; return(x) }

clusterR(nl, calc, args = list(fun=keep82))

library(raster)
rasterOptions(progress = "text", timer = T, chunksize = 1.28e+11, maxmemory = 1.64e+11, memfrac = .8)
aggnl <- aggregate(nl, fact = 33, fun = sum, na.rm = T)

m <- c(0, 81, 0, 81, 82, 1, 83, Inf, 0)
rclmat <- matrix(m, ncol=3, byrow=TRUE)
reclassify(nl, rclmat)

plot(nl)

# newsize <- raster(nrow = 3133, ncol = 4836)
asdf <- projectRaster(nl, res = c(1000,1000), crs = "+proj=aea +lat_0=23 +lon_0=-96 +lat_1=29.5 +lat_2=45.5 +x_0=0 +y_0=0 +datum=WGS84 +units=km +no_defs ")

library(raster)
renl <- resample(nl, asdf)

renl <- clusterR(nl, resample, args = list(y = asdf))



#New attempt 

#county_proj <- spTransform(county_shp, CRSobj = nl@crs)

county_shp <- county_shp %>% mutate(ID = row_number()) %>% 
  select(fips, ID)

ext <- extent(nl)
ext <- paste(ext[1], ext[3], ext[2], ext[4])

res <- paste(res(nl)[1], res(nl)[2])

write_sf(county_shp, "/Volumes/GoogleDrive/My Drive/Personal/Monsanto/data/NLCD_2011_Land_Cover_L48_20190424/county_shp.shp")

system(paste0("/usr/local/Cellar/gdal/3.2.0_2/bin/gdal_rasterize --config GDAL_CACHEMAX 128000 -a ID -te ",
              ext," -tr ", res,
              " -ot Int16 '/Volumes/GoogleDrive/My Drive/Personal/Monsanto/data/NLCD_2011_Land_Cover_L48_20190424/county_shp.shp' '/Volumes/GoogleDrive/My Drive/Personal/Monsanto/data/NLCD_2011_Land_Cover_L48_20190424/county_shp.img'"))










library(gdalUtils)
conflict_prefer("gdal_rasterize", "gdalUtils")
library(parallel)

nsc <- 23
cl <- makeCluster(nsc)

poly_county <- shapefile("/Volumes/GoogleDrive/My Drive/Personal/Monsanto/data/NLCD_2011_Land_Cover_L48_20190424/county_shp.shp")

poly_county %>% select(fips, ID)


r_template <- raster(county_shp, res = c(30, 30)) 

r_template <- raster(extent(nl), res = 30)

clusterExport(cl, c('r_template', 'county_shp'))
clusterEvalQ(cl, sapply(c('raster', 'gdalUtils'), require, char=TRUE))


s <- raster::stack(parLapply(cl, names(county_shp), function(nm) {
  f <- paste0('/Volumes/GoogleDrive/My Drive/Personal/Monsanto/data/NLCD_2011_Land_Cover_L48_20190424/cluster output/county_shp_', nm, '.tif')
  suppressWarnings(writeRaster(r_template, f))
  gdalUtils::gdal_rasterize(src_datasource='county_shp', dst_filename=f, a=nm,
                 output_Raster=TRUE)
}))



r <- raster(res=c(30,30), xmn=0, xmx=161190, ymn=0, ymx=104424, crs="+proj=longlat +datum=WGS84 +no_defs")    
for(n in 1:npar){
  writeRaster(r, filename=paste0(".\\gdal_p",n,".tif"), format="GTiff", overwrite=TRUE)
}                 
                 
                 
                 
cl <- makeCluster(npar)
clusterEvalQ(cl, sapply(c('raster', 'gdalUtils',"rgdal"), require, char=TRUE))
clusterExport(cl, list("r","npar"))
                 
# parallel apply the gdal_rasterize function against the vector parts that were written, 
# same number as processors, against the pre-prepared rasters
parLapply(cl = cl, X = 1:npar, fun = function(x) gdal_rasterize(src_datasource=paste0(".\\parts\\mydata_p",x,".shp"), dst_filename=paste0(".\\gdal_p",n,".tif"),b=1,a="code",verbose=F,output_Raster=T))

# There are now n rasters representing the n segments of the original vector file
# read in the rasters as a list, merge and write to a new tif. 
s <- lapply(X=1:npar, function(x) raster(paste0(".\\gdal_p",n,".tif")))
s$filename <- "myras_final.tif"
do.call(merge,s)
stopCluster(cl)


npar <- detectCores() - 1
features <- 1:nrow(county_shp[,])
parts <- split(features, cut(features, npar))

# write the vector parts out
for(n in 1:npar){
  writeOGR(county_shp[parts[[n]],], ".\\parts", paste0("mydata_p",n), driver="ESRI Shapefile")
}

# set up and write a blank raster for gdal_rasterize for EACH vector segment created above
r <- raster(res=c(30,30), xmn=0, xmx=161190, ymn=0, ymx=104424, crs="+proj=longlat +datum=WGS84 +no_defs")
for(n in 1:npar){
  writeRaster(r, filename=paste0(".\\gdal_p",n,".tif"), format="GTiff", overwrite=TRUE)
}

# set up cluster and pass required packages and objects to cluster
cl <- makeCluster(npar)
clusterEvalQ(cl, sapply(c('raster', 'gdalUtils',"rgdal"), require, char=TRUE))
clusterExport(cl, list("r","npar"))

# parallel apply the gdal_rasterize function against the vector parts that were written, 
# same number as processors, against the pre-prepared rasters
parLapply(cl = cl, X = 1:npar, fun = function(x) gdalUtils::gdal_rasterize(src_datasource=paste0(".\\parts\\mydata_p",x,".shp"),
                                                                dst_filename=paste0(".\\gdal_p",n,".tif"),b=1,a="code",verbose=F,output_Raster=T))

# There are now n rasters representing the n segments of the original vector file
# read in the rasters as a list, merge and write to a new tif. 
s <- lapply(X=1:npar, function(x) raster(paste0(".\\gdal_p",n,".tif")))
s$filename <- "myras_final.tif"
do.call(merge,s)
stopCluster(cl)