Sum Nlayers of a Rasterstack in R

Sum nlayers of a rasterStack in R

You can use stackApply to do this. Using your example data, it looks like the name of each raster layer is the date. You can use that to build the indices that you need to pass to stackApply.

The list of indices needs to have 31 1s for the January days etc.

You could do:

    #get the date from the names of the layers and extract the month
    indices <- format(as.Date(names(data), format = "X%Y.%m.%d"), format = "%m")
    indices <- as.numeric(indices)

    #sum the layers
    datasum<- stackApply(data, indices, fun = sum)

The result will be a raster stack of 12 layers.

To subset raster layers from the stack, you can do data[[c(1,2]]

Sum conditioned through raster stack in R

You need a simple function that computes (in)stability. For example, function f

f <- function(x) sum(diff(x)==0)

computes the number of steps that are not transitions from one state to another.

Try it

f(c(0,0,0,0,0,0))
#[1] 5
f(c(0,0,0,0,1,1))
#[1] 4
f(c(0,1,0,1,0,1))
#[1] 0

Now with Raster data

library(raster)
s <- stack(system.file("external/rlogo.grd", package="raster"))
s <- stack(s, s)
x <- calc(s, f)
plot(x)

Raster series sum

In general, it is recommended not to use loops in R for things that can be easily vectorised. Rather than trying to fix the (several) problems with your loop, I show instead a better way. You can perform the whole calculation in a single vectorised line:

sum(cellStats(myras1==150, stat="sum")) * 100/m

Breaking this down: cellStats performed on a raster stack will return a vector of values, one for each layer. sum then adds these together. Then we divide by number of cells in the whole stack (all layers combined) and multiply by 100 to convert to a percnetage.

Testing this on some reproducible dummy test data:

set.seed(123)
myras1 = list(
  raster(nrows = 100, ncols = 100, vals = sample(140:150,10000,T)),
  raster(nrows = 100, ncols = 100, vals = sample(140:150,10000,T)),
  raster(nrows = 100, ncols = 100, vals = sample(140:150,10000,T)),
  raster(nrows = 100, ncols = 100, vals = sample(140:150,10000,T))
)
myras1 = stack(myras1)
m = ncol(myras1) * nrow(myras1) * nlayers(myras1)

sum(cellStats(myras1==150, stat="sum")) * 100/m
# [1] 8.815

Sum pixel values in a Raster stack based on another raster stack

The last part of your example did not work and I changed it to this (for one year)

x <- rep(NA, ncell(s))
for (i in 1:ncell(s)) {
  x[i] <- sum(s[i][ss.1[i]:se.1[i]], na.rm = T)
}
x <- setValues(ss.1, x)
x
#class      : RasterLayer 
#dimensions : 5, 5, 25  (nrow, ncol, ncell)
#resolution : 14, 8  (x, y)
#extent     : -150, -80, 20, 60  (xmin, xmax, ymin, ymax)
#crs        : +proj=longlat +datum=WGS84 +no_defs 
#source     : memory
#names      : layer 
#values     : 0.6505058, 10.69957  (min, max)

You can get that result like this

idx <- stack(ss.1, se.1)
thefun <- function(x, y){
    apply(cbind(y, x), 1, function(i) sum(i[(i[1]:i[2])+2], na.rm = T))
}  
z <- overlay(s, idx, fun=thefun)
    

There are more examples here for a similar question.

Given that this is a general problem, I have added a function rapp (range-apply) for it in terra (the replacement for raster) --- available here; this should be on CRAN in early July.

library(terra)
r <- rast(ncols=5, nrows=5, xmin=-150, xmax=-80, ymin=20, ymax=60)
values(r) <- 1:ncell(r)
s <- rast(replicate(36, r))

ss.1 <- r 
values(ss.1) <- as.integer(runif(ncell(ss.1), min=1, max=72))
se.1 <- ss.1+10 

x <- rapp(s, ss.1, se.1, sum)
    

Function to sum each grid cells of raster stack using other rasters as an indicator

Your approach

x <- s$layer.1

system.time(
for (i in 1:ncell(x)) {
     x[i] <- sum(s[[r_start[i]:r_end[i]]][i], na.rm = T)
   }
)

   user  system elapsed 
  0.708   0.000   0.710

My proposal

You can add the rasters used as indices at the end of your stack and then use calc to highly speed up the process (~30-50x).

s2 <- stack(s, r_start, r_end)
sum_time <- function(x) {sum(x[x[6]:x[7]], na.rm = T)}

system.time(
   output <- calc(s2, fun = sum_time)
)

   user  system elapsed 
  0.016   0.000   0.015

all.equal(x, output)

[1] TRUE

Sample Data

library(raster)

# Generate rasters of random values
r1 <- r2 <- r3 <- r4 <- r5 <- r_start <- r_end <- raster(ncol=10, nrow=10)

r1[] <- rnorm(ncell(r1), 1, 0.2)
r2[] <- rnorm(ncell(r2), 1, 0.2)
r3[] <- rnorm(ncell(r3), 1, 0.2)
r4[] <- rnorm(ncell(r4), 1, 0.2)
r5[] <- rnorm(ncell(r5), 1, 0.2)
s <- stack(r1,r2,r3,r4,r5)

r_start[] <- sample(1:2, ncell(r_start),replace = T)
r_end[] <- sample(3:5, ncell(r_end),replace = T)

---

Related Topics