问题描述
我想创建栅格堆栈的子集,并在上一层和下一层之间的差异全部为 NA
时将它们作为新堆栈写入。即,从第 1 层开始,我想创建一个栅格堆栈的子集,直到前一层和下一层之间没有重叠像素(即两层之间的差异都是 NA
) 所以我想要的是;从第 1 层开始,保留上一层和下一层之间至少有 1 个公共像素的所有层,将它们写为 1 堆栈,然后移动到下一层。以下是示例数据和不成功的 for 循环。在这个例子中,我想保留层 1:8,命名并编写它们,然后从层 9 重新开始,依此类推。
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,runif(ncell(r)))))
s1e<-extend(s2,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
解决方法
这样的事情可能对你有用。
您的示例数据
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,runif(ncell(r)))))
s1e <- extend(s2,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,126 (nrow,ncell)
#resolution : 72,ymax)
#crs : +proj=longlat +datum=WGS84 +no_defs
#source : memory
#names : layer.1.2
#values : NA,NA (min,max)
#
#[[3]]
#class : RasterLayer
#dimensions : 14,ymax)
#crs : +proj=longlat +datum=WGS84 +no_defs
#source : memory
#names : layer.2.2
#values : NA,max)