提问人:madip 提问时间:10/13/2023 更新时间:10/13/2023 访问量:33
提取网格对象下的单元格值和单元格比例,然后在 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)
问:
我正在使用一个 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
答:
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) 值的逻辑(布尔)栅格,以便您可以取像元值的平均值来获得比例。
评论