提取网格对象下的单元格值和单元格比例，然后在Terra（R中）中查找高于某个值的单元格比例

Question

我正在使用代表草本覆盖百分比的 30m 分辨率栅格。我试图确定每个 1km 网格单元内值大于 10（代表 10% 草本覆盖）的单元数量。网格对象是一个 sf 数据框，其中每行代表一个 1km x 1km 多边形，它有两列，cellid 和 Geometry。

我有一个代表每年的 SpatRaster（堆栈）列表，它们与 sf 网格对象处于同一投影中，并且我尝试在循环中运行此函数以提取每年的每个网格单元。最终，我想要的输出是每年的数据框，其中一列代表网格 cellid，另一列代表 1 公里网格内所有 3000 万个细胞中草本覆盖率大于 10% 的比例。

这是我目前正在使用的循环；然而，它 a) 看起来非常非常慢，b) 当我尝试在 tortgrid_1km_strata 对象的整个长度（而不仅仅是 1:100）运行它时，它似乎正在中断。它抛出一个错误，显示“data.frame(layer = i,aggregate(w,list(x),sum,na.rm = FALSE))中的错误： 参数意味着不同的行数：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

您可能可以通过聚合栅格来解决这个问题。类似于我下面显示的内容。

示例数据

library(terra)
r <- rast(res=30, xmin=0, xmax=3000, ymin=0, ymax=3000)
set.seed(1)
values(r) <- runif(ncell(r), 1, 100)

解决方案

a <- aggregate(r > 10, 33, mean, na.rm=TRUE)