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

Question

I have a raster stack representing Evapotranspiration (ET) with 396 layers (3 raster layers for a month for 11 years in total - 2009 to 2019). For each month the raster layer always represents 1st, 11th and 21st day of the month called dekads. Here is the sample dataset

library(raster)

#create a raster with random numbers
r <- raster(ncol=5, nrow=5, xmx=-80, xmn=-150, ymn=20, ymx=60)
values(r) <- runif(ncell(r))

#create a random raster stack for 3 raster a month for 11 years
n <- 396 #number of raster
s <- stack(replicate(n, r)) # convert to raster stack

#rename raster layers to reflect date
d =rep(c(1,11,21),132)
m =rep(1:12, 11, each =3)
y = rep (2009:2019, each =36)
df.date <- as.Date(paste(y, m, d,sep="-"), "%Y-%m-%d")
names(s) = df.date

I also have two other raster stacks with pixel values representing Season start (ss) 11 layers and season end (se)11 layers for years 2009 to 2019.

#create a raster stack representing season start (ss) and season end (se)
# The pixel value represents dekad number. Each raster layer covers exactly three calendar years with the target year in the middle.
# (1-36 for the first year, 37-72 for the target year, 73-108 for the next year). 
ss.1 = r # season start raster
values(ss.1)= as.integer(runif(ncell(ss.1), min=1, max=72))
se.1 = ss.1+10 # season end raster
yr = 11
ss <- stack(replicate(yr, ss.1)) # season start raster stack
se <- stack(replicate(yr, se.1)) #season end rasterstack

Now I need to estimate seasonal sum for each year from the "s" raster stack such that the time period for each pixels to sum should correspond to pixel values from "ss" and "se" by considering a 3 year moving window.

Here is an example of output I need for one time step (3yr window) with one season start (ss) raster and one season end (se) raster. But really struck at looping through three raster stacks (s - representing dataset, ss -representing season start date and se -representing season end date). Grateful for any help.

# Example to calculate pixel based sum for 1 time step
#subset first 3 years - equal to 108 dekads
s.sub = subset(s, 1:108)
# sum each grid cells of "s" raster stack using "ss.1" and "se.1" as an indicator for the three year subset.
for (i in 1:ncell(s.sub)) {
  x[i] <- sum(s[[ss.1[i]:se.1[i]]][i], na.rm = T)
}

Please provide a minimal, self-contained, reproducible example. In other words, do not refer to files you have but either use files that come with R or create some values in memory. Keep the example as small as possible; and show for at least one grid cell what the input and the expected output is. — Robert Hijmans
Sorry Robert Hijmans for being not including a working example. Now edited it with an working example. — karthikeyan matheswaran

Robert Hijmans Robert Hijmans · Accepted Answer · 2020-06-13T17:38:30

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)

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

1 Answers