CODEWITHSUNDEEP
rrasterr-raster
Lily Nature
Lily Nature

Reputation: 695

Subset R rasterstack based on difference in raster layers grid cell numbers

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

Upvotes: 0

Views: 430

Answers (1)

Robert Hijmans
Robert Hijmans

Reputation: 47146

Something like this may work for you.

Your example data

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)

Solution:

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)

Upvotes: 1

Related Questions