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

Question

我想创建栅格堆栈的子集，并将它们写为新的堆栈，当上一层和下一层之间的差异在每个栅格层块之后都是 NA 时。如果没有团块，我将通过遵循罗伯特在中的回答（如下面的脚本）来实现这一点。但是，我也想运行通过考虑团块来做到这一点。每层可能有 1 或 2 个团块。因此，从下面示例数据堆栈中的 layer 1 开始，我想确定团块编号，并为每个团块创建一个栅格堆栈子集，直到上一层和下一层之间没有重叠像素（即，两层之间的差异都是NA）。所以我想要的是；从layer 1开始，对于每个clump，保留所有在上一层和下一层之间至少有一个公共像素的层，将它们写成1个stack，然后移动到下一层。在示例 r_stk 中，我想为 clump 1（顶部）保留层 1:8，将它们分配为 1 个堆栈，运行为 clump 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

如果我移除较小的团块并且每层只保留一个团块，效果很好我想运行为每一层上的每个团块做这个我想，我正在尝试运行在脚本之上增加一个 for 循环低于此考虑团块但无法成功运行

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]]

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

Subset R rasterstacks based on difference in grid cell number on each clumps

r

raster

r-raster