Search code examples
rxtszoor-raster

R How to calculate the seasonal average value for each year using stackApply?


I wanted to calculate the seasonal average value for each year and not for the entire period. I defined my seasons as follows: DJF (December-February), MAM (March-May), JJA (June-August) and SON (September-November).

Inspired by the solution of the question of Fredrick, I created a index "groups" which represents the seasons then I applied the command "stackApply" but this one calculates the average seasonal value for the all of period. I explain, the final layer obtained contains only 4 raster but for my case I wanted to calculate "the seasonal average value for each year, therefore it is necessary to have 4 rasters for each year and in total the rasterstack should have 136 rasters.

Below my code

Thanks for your help

library(raster)
set.seed(123)
r <- raster(ncol=10, nrow=10)
r_brick <- brick(sapply(1:408, function(i) setValues(r, rnorm(ncell(r), i, 3))))
dim(r_brick)

dates <- seq(as.Date("1982-01-01"), as.Date("2015-12-31"), by="month")  
months <- format(dates, "%Y-%m")

groups <- function(x) {
  d <- as.POSIXlt(x)

  ans <- character(length(x))
  ans[d$mon %in%  c(11,0:1)] <- "DJF"
  ans[d$mon %in%  2:4] <- "MAM"
  ans[d$mon %in%  5:7] <- "JJA"
  ans[d$mon %in% 8:10] <- "SON"
  ans
}

data.frame(dates, groups(dates))

r_brick.s <- stackApply(r_brick, indices=groups(dates), fun=mean,na.rm=TRUE) 
nlayers(r_brick.s)

Solution

  • Your example data

    library(raster)
    r <- raster(ncol=10, nrow=10)
    b <- brick(sapply(1:408, function(i) setValues(r, rnorm(ncell(r), i, 3))))    
    dates <- seq(as.Date("1982-01-01"), as.Date("2015-12-31"), by="month")  
    

    As Majid points out, if you want to group by years, you need to use these

    years <- as.integer(format(dates, "%Y"))
    months <- as.integer(format(dates, "%m"))
    

    Now the months need to be grouped. Note that as you start with December you must make sure to not combine January and December of the same year. Rather you want to combine December of year i with January and February of year i+1. Here is one way to accomplish that (make the year start in December!)

    n <- length(months)
    
    # move all months back one month   
    mnt <- c(months[-1], ifelse(months[n] < 12, months[n]+1, 1)) 
    
    # move the years along
    yrs <- c(years[-1],  ifelse(months[n] < 12, years[n], years[n]+1)) 
    
    # group by trimesters using integer division (or do: floor((mnt-1) / 3))
    trims <- (mnt-1) %/% 3  
    
    # get names instead of 0, 1, 2, 3
    trimnms <- c("DJF", "MAM", "JJA", "SON")[trims + 1]
    

    Combine years and names

    yt <- paste(yrs, trimnms, sep="_")
    

    Use the index

    s <- stackApply(b, indices=yt, fun=mean, na.rm=TRUE) 
    

    If the above business of moving the months backward is difficult to follow, try it out with just a few dates (the first 15 or so)