Loop through time values in a raster brick in R

Question

I have raster bricks of temperature and rainfall consisting of one layer for each month of the year. The raster brick looks like this:

> tmx
class      : RasterBrick
dimensions : 3875, 8017, 31065875, 240  (nrow, ncol, ncell, nlayers)
resolution : 0.04490319, 0.04490319  (x, y)
extent     : -180.0393, 179.9495, -90.04088, 83.95898  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
source     : /scratch/bntjoa002/chelsa_cmip5_ts/tasmax/CHELSAcmip5ts_tasmax_ACCESS1-3_rcp8                                      5_2070-2089.nc
names      : X2070.01.15, X2070.02.15, X2070.03.15, X2070.04.15, X2070.05.15, X2070.06.15,                                       X2070.07.15, X2070.08.15, X2070.09.15, X2070.10.15, X2070.11.15, X2070.12.15, X2071.01.15                                      , X2071.02.15, X2071.03.15, ...
Date       : 2070-01-15, 2089-12-15 (min, max)
varname    : air_temperature

I have been using the following code to get the average monthly values for a specified period (e.g. resulting in 12 averaged layers, one for each month, over a X year period):

new_tmx <- subset(tmx, which(getZ(tmx) >= as.Date("2070-01-15") & getZ(tmx) <= as.Date("2089-12-15"))) 
indices <- format(as.Date(names(new_tmx), format = "X%Y.%m.%d"), format = "%m")
indices <- as.numeric(indices)
Monthlytmx<- stackApply(new_tmx, indices, fun = mean)

I do this separately for the maximum temperature, minimum temperature, and precipitation bricks.

What I need now is to do the same process, but instead to output one raster layer for each year consisting of the averaged values for one-year periods (i.e. for 2070, 2071, 2072...). I've only done very simple loops in R, and I'm not sure how to deal with the time/date function in this example?

Answer 1

I think what you need to do is create a new index from the year in the date - you can use the "format" function on a date to get the year out:

> indices = format(getZ(tmx),"%Y")
> indices
  [1] "2070" "2070" "2070" "2070" "2070" "2070" "2070" "2070" "2070" "2070"
 [11] "2070" "2070" "2071" "2071" "2071" "2071" "2071" "2071" "2071" "2071"
 [21] "2071" "2071" "2071" "2071" "2072" "2072" "2072" "2072" "2072" "2072"
...

then stackApply with that:

> Yearlytmx<- stackApply(new_tmx, indices, fun = mean)
> Yearlytmx
class      : RasterBrick 
dimensions : 3, 4, 12, 20  (nrow, ncol, ncell, nlayers)
resolution : 0.25, 0.3333333  (x, y)
extent     : 0, 1, 0, 1  (xmin, xmax, ymin, ymax)
crs        : NA 
source     : memory
names      : index_2070, index_2071, index_2072, index_2073, index_2074, index_2075, index_2076, index_2077, index_2078, index_2079, index_2080, index_2081, index_2082, index_2083, index_2084, ... 
min values :  0.2871428,  0.3901444,  0.3483226,  0.2906245,  0.2947116,  0.3599598,  0.3503530,  0.3295298,  0.4077774,  0.3867792,  0.3916475,  0.3714683,  0.3746881,  0.3817867,  0.3694289, ... 
max values :  0.6998091,  0.6197922,  0.6100803,  0.5943091,  0.6323317,  0.6812757,  0.7028588,  0.6113572,  0.5682294,  0.6368843,  0.6010818,  0.6163263,  0.6534441,  0.5922601,  0.5980651, ... 

That brick now has 20 layers, which should be one for each year.

You can then add the year as a Z to the data:

> Yearlytmx = setZ(Yearlytmx,unique(indices),"Year")

Or use the "names", which have an "index_" attached to them:

> names(Yearlytmx)
 [1] "index_2070" "index_2071" "index_2072" "index_2073" "index_2074"
 [6] "index_2075" "index_2076" "index_2077" "index_2078" "index_2079"
[11] "index_2080" "index_2081" "index_2082" "index_2083" "index_2084"
[16] "index_2085" "index_2086" "index_2087" "index_2088" "index_2089"