基于每个团块上网格单元数差异的子集 R 栅格堆栈 [英] Subset R rasterstacks based on difference in grid cell number on each clumps

问题描述

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

I want to create subsets of raster stacks and write them as new stacks when the difference between the previous layer and the next layer is all NA following each clump of the raster layers. Without clumps, I would achieve this by following Robert's answer in this question ( as below in script). However, I want to run this by considering the clumps too. There may be 1 or 2 clumps in each layer. So starting from layer 1 in the example data stack below, I want to identify the clumps numbers and for each clump, create a subset of raster stacks until there are no-overlapping pixels between the previous and next layers (i.e., the difference between the two layers is all NA). So I want is; starting from layer 1, for each clump, retain all the layers that have at least 1 common pixel between the previous and next layer, write them as a 1 stack, and move to the next. In the sample r_stk, I want to retain layers 1:8 for clump 1 (top) assign them as 1 stack, run for clump 2 (bottom), and again retain layers 1:5 assign them as a new stack, and so on. Below are the sample data and code that would work fine following this answer if there would be no clumps.

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 循环下面考虑团块但无法使其成功运行

Works fine if i remove smaller clump and keep only one clump per layer however i want to run this for each clump on each layer I guess, i am trying to run one additional for loop on top of the script below that consider clumps but could not make it successfully run

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,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))

没有团块(即单个团块)的解决方案与上一个答案相同，但包装成函数

The solution without clumps (that is, for a single clump) as in the previous answer, but wrapped into a function

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

Get the clumps and their unique IDs

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

屏蔽单个丛的数据的函数

A function that masks out the data for a single clump

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

为每个 ID 调用 f

x <- lapply(u, f)

x 是一个列表.每个元素都是一团的结果

x is a list. Each element is the result for one clump

length(x) 
#2

团块#1 的列表

r1 <- x[[1]]

这篇关于基于每个团块上网格单元数差异的子集 R 栅格堆栈的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！