基于每个团块上网格单元数差异的子集R栅格堆栈

Question

我想创建栅格堆栈的子集，并在前一层和下一层之间的差值为每个栅格层簇之后的所有NA时将它们写入新堆栈。如果没有团块，我将按照Robert对这个问题的回答来实现这个目标(如下面的脚本)。然而，我也想通过考虑团块来运行它。每层可能有1或2个团块。因此，从下面示例数据堆栈中的layer 1开始，我想要识别团块数量，并为每个团块创建栅格堆栈的子集，直到前一层和下一层之间没有重叠像素(即，两层之间的差异都是NA)。所以我想要的是;从layer 1开始，对于每个团块，保留所有在前一层和下一层之间至少有1个公共像素的层，将它们写为1堆栈，并移动到下一层。在示例r_stk中，我希望将层1(顶部)保留为1:8，将它们分配为一个堆栈，运行为簇2(底部)，并再次保留层1:5，将它们分配为一个新堆栈，依此类推。下面是示例数据和代码，如果没有团块，它们将按照这个答案工作得很好。

library(raster)
library(tidyverse)
#Create null raster, fill values and get stack with clumps 
r<-raster(extent(-180,-140,34,83),
crs="+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0",
resolution=10, vals=NULL)
r
#Make series of raters with clumps and stack
r1<-r
values(r1)<-c(70,69,NA,NA,73,70,NA,NA,NA,NA,NA,NA,100,99,NaN,NA,101,99,76,NA)
r2<-r
values(r2)<-c(89,81,72,NA,87,77,69,NA,NA,NA,NA,NA,89,99,NaN,NA,89,100,84,NA)
r3<-r
values(r3)<-c(112,103,86,76,90,82,78,NaN,NA,NA,NA,NA,79,93,NaN,NA,78,93,88,NA)
r4<-r
values(r4)<-c(125,115,98,88,84,81,82,NaN,NA,NA,NA,NA,69,81,NaN,NA,69,80,83,NA)
r5<-r
values(r5)<-c(132,125,110,100,77,76,82,NaN,NA,NA,NA,NA,NaN,71,NaN,NA,70,71,74,NA)
r6<-r
values(r6)<-c(118,114,103,93,72,75,77,NaN,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA)
r7<-r
values(r7)<-c(98,92,76,69,70,70,76,NaN,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA)
r8<-r
values(r8)<-c(76,73,68,NA,76,73,NaN,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA)
r9<-r
values(r9)<-c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA)
r_stk<-stack(r1,r2,r3,r4,r5,r6,r7,r8,r9)
plot(r_stk) #raster stack with clumps

如果我移除较小的团块并且每层只保留一个团块，则可以正常工作我想对每一层上的每一个团块运行这个我想，我正试图运行一个额外的循环上的脚本

singleclump_lst<-list()
for (i in 1: nlayers(r_stk)){
rasi<-subset(r_stk,i)
#Classify based on clumps
clumps<-clump(rasi,directions=8)
clumpFreq2 <- as.data.frame(freq(clumps))
clumpFreq_na2<-clumpFreq2%>%
drop_na()
clumpFreq_na2
excludeID_i <-clumpFreq_na2$value[which(clumpFreq_na2$count == max(clumpFreq_na2$count))]
excludeID_i
subNA_i <- function(a, b) {
a[!b %in% excludeID_i] <- NA
return(a)}
rasclmp_i<-overlay(rasi,clumps,fun=subNA_i)

singleclump_lst[[i]]<-rasclmp_i
}
rr_stk<-stack(singleclump_lst)  
rr_stk
plot(rr_stk)
out <- lst <- list()
nc <- ncell(rr_stk)   
for (i in 1:nlayers(rr_stk)) {
if (i==1) {
j <- 1
s <- rr_stk[[i]]
} else {
s <- s + rr_stk[[i]]
}
if (freq(s, value=NA) == nc) {
ii <- max(j, i-1)   
out <- c(out, rr_stk[[j:ii]])
s <- rr_stk[[i]]
j <- i
}
}
out <- c(out, rr_stk[[j:i]])
out

Answer 1

您的示例数据

library(raster)
b <- brick(extent(-180,-140,34,83), nrow=5, ncol=4,
crs="+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")
values(b) <- cbind(
c(70,69,NA,NA,73,70,NA,NA,NA,NA,NA,NA,100,99,NaN,NA,101,99,76,NA),
c(89,81,72,NA,87,77,69,NA,NA,NA,NA,NA,89,99,NaN,NA,89,100,84,NA),
c(112,103,86,76,90,82,78,NaN,NA,NA,NA,NA,79,93,NaN,NA,78,93,88,NA),
c(125,115,98,88,84,81,82,NaN,NA,NA,NA,NA,69,81,NaN,NA,69,80,83,NA),
c(132,125,110,100,77,76,82,NaN,NA,NA,NA,NA,NaN,71,NaN,NA,70,71,74,NA),
c(118,114,103,93,72,75,77,NaN,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA),
c(98,92,76,69,70,70,76,NaN,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA),
c(76,73,68,NA,76,73,NaN,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA),
c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA))

没有团块(即单个团块)的解决方案，就像前面的答案一样，但是包装成一个函数

one_clump <- function(r_stk) {
out <- lst <- list()
nc <- ncell(r_stk)   
for (i in 1:nlayers(r_stk)) {
if (i==1) {
j <- 1
s <- r_stk[[i]]
} else {
s <- s + r_stk[[i]]
}
if (freq(s, value=NA) == nc) {
ii <- max(j, i-1)   
out <- c(out, r_stk[[j:ii]])
s <- r_stk[[i]]
j <- i
}
}
out <- c(out, r_stk[[j:i]])
out
}

获取团块及其唯一id

clm <- clump(b[[1]])
u <- unique(clm)

屏蔽单个簇数据的函数

f <- function(i) {
rr <- clm == i
bb <- mask(b, rr, maskvalue=0)
one_clump(bb)
}

为每个ID调用f

x <- lapply(u, f)

x是一个列表。每个元素是一个簇

的结果。

length(x) 
#2

集群#1的列表

r1 <- x[[1]]