基于栅格层网格单元数差异的子集 R rasterstack [英] Subset R rasterstack based on difference in raster layers grid cell numbers

问题描述

我想创建栅格堆栈的子集，并在上一层和下一层之间的差异全部为NA时将它们写为新的堆栈.即，从第1层开始，我想创建一个光栅堆栈的子集，直到前一层和下一层之间没有重叠像素(即两层之间的差异都是NA)所以我要的是；从第 1 层开始，保留上一层和下一层之间至少有 1 个公共像素的所有层，将它们写为 1 堆栈，然后移动到下一层.下面是一个示例数据和不成功的 for 循环.在这个例子中，我想保留层 1:8，命名并编写它们，然后从层 9 重新开始，依此类推.

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. I.e., starting from layer 1, I want to 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, 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. Below are a sample data and unsuccessful for-loop. In this example, I want to retain layers 1:8, name and write them and start again from layer 9 and so on.

r <- raster(ncol=5, nrow=5)
set.seed(0)
#create raster layers with some values 
s <- stack(lapply(1:8, function(i) setValues(r, runif(ncell(r)))))
s1<-extend(s,c(-500,100,-400,100))

#to recreate the condition I am looking for, create 2 layers with `NA` vlaues
s2 <- stack(lapply(1:2, function(i) setValues(r, runif(ncell(r)))))
s1e<-extend(s2,c(-500,100,-400,100))
s1e[]<-NA

#Stack the layers
r_stk<-stack(s1,s1e)
plot(r_stk)

#here is the sample code showing what i am expecting here but could not get

required_rst_lst<-list() # sample list of raster layers with overlapping pixels I am hoping to create
for ( i in 1: nlayers(r_stk))
  # i<-1
  lr1<-subset(r_stk,i)
lr1
lr2<-subset(r_stk,i+1)
lr2
diff_lr<-lr1-lr2  
plot(diff_lr)

if ((sum(!is.na(getValues(diff_lr)))) ==0)) #??

required_rst_lst[[i]] #?? I want layers 1: 8 in this list 
#because the difference in these layers in not NA

推荐答案

这样的事情可能对你有用.

Something like this may work for you.

您的示例数据

library(raster)
r <- raster(ncol=5, nrow=5)
set.seed(0)
s <- stack(lapply(1:8, function(i) setValues(r, runif(ncell(r)))))
s1 <- extend(s,c(-500,100,-400,100))

s2 <- stack(lapply(1:2, function(i) setValues(r, runif(ncell(r)))))
s1e <- extend(s2,c(-500,100,-400,100))
values(s1e) <- NA 
r_stk <- stack(s1,s1e)

解决方案:

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

#[[1]]
#class      : RasterStack 
#dimensions : 14, 9, 126, 8  (nrow, ncol, ncell, nlayers)
#resolution : 72, 36  (x, y)
#extent     : -468, 180, -414, 90  (xmin, xmax, ymin, ymax)
#crs        : +proj=longlat +datum=WGS84 +no_defs 
#names      :  layer.1.1,  layer.2.1,    layer.3,    layer.4,    layer.5,    layer.6,    layer.7,    layer.8 
#min values : 0.06178627, 0.01339033, 0.07067905, 0.05893438, 0.01307758, 0.03554058, 0.06380848, 0.10087313 
#max values :  0.9919061,  0.8696908,  0.9128759,  0.9606180,  0.9926841,  0.9850952,  0.8950941,  0.9437248 
#
#[[2]]
#class      : RasterLayer 
#dimensions : 14, 9, 126  (nrow, ncol, ncell)
#resolution : 72, 36  (x, y)
#extent     : -468, 180, -414, 90  (xmin, xmax, ymin, ymax)
#crs        : +proj=longlat +datum=WGS84 +no_defs 
#source     : memory
#names      : layer.1.2 
#values     : NA, NA  (min, max)
#
#[[3]]
#class      : RasterLayer 
#dimensions : 14, 9, 126  (nrow, ncol, ncell)
#resolution : 72, 36  (x, y)
#extent     : -468, 180, -414, 90  (xmin, xmax, ymin, ymax)
#crs        : +proj=longlat +datum=WGS84 +no_defs 
#source     : memory
#names      : layer.2.2 
#values     : NA, NA  (min, max)

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