I am working with a 30m resolution raster representing percentage herbaceous cover. I am trying to determine the number of cells within each 1km grid cell that has a value greater than 10 (representing 10% herbaceous cover). The grid object is a sf data frame where each row represents a 1km x 1km polygon and it has two columns, cellid and geometry.
I have a list of SpatRasters (stack) representing each year that are in the same projection as the sf grid object, and I am trying to run this function within a loop to extract for each gridcell in each year. Ultimately, the output I want is a data frame for each year with a column representing grid cellid and another representing the proportion of all 30m cells within the 1km grid that are greater than 10% herb cover.
This is the loop I currently am working with; however, it is a) seemingly very very slow and b) seems to be breaking when I try to run it for the entire length of the tortgrid_1km_strata object (rather than just 1:100). It throws an error that says "Error in data.frame(layer = i, aggregate(w, list(x), sum, na.rm = FALSE)) : arguments imply differing number of rows: 1, 0".
herb_10_LIST <- vector("list",length(stack))
for (i in 1:nlyr(stack)){ ## Just running for first year raster
herb_10 <- c() ## Making an empty vector to store the proportion values with >10% herb cover for each grid cell
for (j in 1:nrow(tortgrid_1km_strata)){
cell <- terra::vect(tortgrid_1km_strata[j,]) # 1 km grid cell
ext <- terra::extract(x = stack[[i]], y = cell, fun=table, weights=TRUE, exact=FALSE) # Extract 30m cells under 1km grid and find vlaues and weights
vals <- colnames(ext) # For some reason the way the table comes out, the raster values are the column names...
## If there are no raster cells under the 1km grid, paste NA
if (length(vals) <=2){ ## The first two columns are not actual raster values!
herb_10[j] <- NA }
## Otherwise,
else {
vals <- as.numeric(vals[3:length(vals)])
perc <- ext[1,] ## The first row of the output table is the cell weights
perc <- perc[,c(3:ncol(perc))] %>% as.numeric()
tab <- as.data.frame(cbind(vals,perc)) # Making into a df of values and weights for each 30 cell
length <- as.numeric(nrow(tab)) # Number of cells under the grid
perc_10 <- filter(tab, vals >= 10) # Filtering for cells > 10% grass cover
## If there is at least one cell with > 10% grass cover
if (nrow(perc_10) >= 1) {
prop_10 <- as.numeric(unname(colSums(perc_10)))
prop_10 <- prop_10[2]/ length } # Finding prop of all cells that had > 10% cover
## Otherwise, assign prop as 0
else { prop_10 <- 0 }
## Put the proportion into the vector
herb_10[j] <- prop_10
}
} # End of inner loop
herb_10_LIST[[i]] <- as.data.frame(cbind(cellids,herb_10)) ## Binding to cellids (object I created that represents just the gridcell ids
} # End of outer loop
You could probably approach this problem by aggregating the raster. Something along the lines of what I show below.
Example data
library(terra)
r <- rast(res=30, xmin=0, xmax=3000, ymin=0, ymax=3000, nlyr=2)
set.seed(1)
values(r) <- runif(size(r), 1, 100)
Solution
a <- aggregate(r > 10, 33, mean, na.rm=TRUE)
If you have to use polygons (it is much better not to)
p <- as.polygons(rast(a))
e <- extract(r > 10, p, fun=mean, na.rm=TRUE, bind=TRUE)
plot(e, 1)
As you can see, one important simplification is to first compute which cells have values larger than 10. That returns a logical (Boolean) raster with TRUE (1) / FALSE (0) values such that you can take the mean of the cell values to get a proportion.