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

编程问答 2022-05-24

问题描述

我想创建光栅堆栈的子集，并在前一层和下一层之间的差异在每个光栅层簇之后都是 NA 时将它们写为新堆栈。如果没有团块，我会按照罗伯特在 this question 中的回答来实现这一点（如下脚本中所示）。但是，我也想通过考虑团块来运行它。每层可能有 1 或 2 个团块。因此，从下面示例数据堆栈中的 layer 1 开始，我想确定团块编号，并为每个团块创建一个栅格堆栈子集，直到前一层和下一层之间没有重叠像素（即两层之间的差异都是NA）。所以我想要的是；从 layer 1 开始，对于每一块，保留上一层和下一层之间至少有 1 个公共像素的所有层，将它们写为 1 堆栈，然后移动到下一层。 在示例 r_stk 中，我想保留层 1:8 为丛 1（顶部）将它们分配为 1 个堆栈，为丛 2（底部）运行，并再次保留层 1:5 将它们分配为新堆栈， 等等。 下面是示例数据和代码，如果没有团块，它们可以在 this answer 之后正常工作。

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",resolution=10,vals=NULL)
r
#Make series of raters with clumps and stack
r1<-r
values(r1)<-c(70,69,NA,73,70,100,99,NaN,101,76,NA)
r2<-r
values(r2)<-c(89,81,72,87,77,89,84,NA)
r3<-r
values(r3)<-c(112,103,86,90,82,78,79,93,88,NA)
r4<-r
values(r4)<-c(125,115,98,80,83,NA)
r5<-r
values(r5)<-c(132,125,110,71,74,NA)
r6<-r
values(r6)<-c(118,114,75,NA)
r7<-r
values(r7)<-c(98,92,NA)
r8<-r
values(r8)<-c(76,68,NA)
r9<-r
values(r9)<-c(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

解决方法

您的示例数据

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