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

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

提问人:madip 提问时间:10/13/2023 更新时间:10/13/2023 访问量:33

问:

我正在使用一个 30m 分辨率的栅格,表示草本覆盖百分比。我正在尝试确定每个 1km 网格像元中值大于 10(代表 10% 草本覆盖)的像元数。网格对象是一个 sf 数据框,其中每行表示一个 1km x 1km 的多边形,它有两列,cellid 和 geometry。

我有一个 SpatRasters(堆栈)列表,代表每年与 sf 网格对象处于同一投影中,我正在尝试在循环中运行此函数以提取每年的每个网格单元。最终,我想要的输出是每年的数据框,其中一列表示网格细胞,另一列表示 1km 网格内所有 30m 细胞中大于 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

循环 if-语句 提取 栅格 terra

评论

0赞 Robert Hijmans 10/13/2023
您能否编辑您的问题并提供一个独立的、最小的、可重复的例子。自包含:使用 R 附带的一些示例数据或使用代码创建它。越小越好。显示您期望的输出。最小:显示您想在循环外执行的操作。如果您的问题是关于编写循环,请删除所有空间数据复杂性。

答:

0赞 Robert Hijmans 10/13/2023 #1

您可以通过聚合栅格来解决此问题。与我在下面展示的内容类似的东西。

示例数据

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)

溶液

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

如果您必须使用多边形(最好不要使用)

p <- as.polygons(rast(a))
e <- extract(r > 10, p, fun=mean, na.rm=TRUE, bind=TRUE)
plot(e, 1)

如您所见,一个重要的简化是首先计算哪些单元格的值大于 10。这将返回一个具有 TRUE (1) / FALSE (0) 值的逻辑(布尔)栅格,以便您可以取像元值的平均值来获得比例。