Extracting cell values and proportion of cell under grid object and then finding proportion of cells above a certain value in Terra (in R)

Question

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

Answer 1

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.