I'm starting from a Rasterbrick like that called r.brick:

class       : RasterBrick 
dimensions  : 17, 19, 323, 21915  (nrow, ncol, ncell, nlayers)
resolution  : 0.11, 0.11  (x, y)
extent      : 8.985, 11.075, 51.325, 53.195  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
data source : E:\Masterdaten\Rasterfiles\sfcwind\bc_raster_metropreg\windleistungsdichte\CNRM.CERFACS.CNRM.CM5_rcp45_r1i1p1_CLMcom.CCLM4.8.17_v1_day_WED_90m_ref_1976_2005_scen_2071_2100_bc.tif 
names       : CNRM.CERF//_2100_bc.1, CNRM.CERF//_2100_bc.2, CNRM.CERF//_2100_bc.3, CNRM.CERF//_2100_bc.4, CNRM.CERF//_2100_bc.5, CNRM.CERF//_2100_bc.6, CNRM.CERF//_2100_bc.7, CNRM.CERF//_2100_bc.8, CNRM.CERF//_2100_bc.9, CNRM.CERF//2100_bc.10, CNRM.CERF//2100_bc.11, CNRM.CERF//2100_bc.12, CNRM.CERF//2100_bc.13, CNRM.CERF//2100_bc.14, CNRM.CERF//2100_bc.15, ... 
min values  :          7.152665e+01,          1.498034e+02,          5.482229e+02,          2.722411e+02,          2.722821e+02,          1.871505e+02,          3.891109e+02,          6.919453e+02,          2.884659e+02,          1.538994e+02,          5.594246e+01,          2.979832e-01,          2.467846e+01,          2.492913e+01,          1.861346e+01, ... 
max values  :          4.258338e+02,          7.356049e+02,          1.613744e+03,          7.528055e+02,          7.839373e+02,          8.567026e+02,          9.895955e+02,          1.355797e+03,          1.178775e+03,          4.890604e+02,          2.062708e+02,          1.725067e+01,          1.128135e+02,          3.015463e+02,          1.362133e+02, ... 

Its my goal to calculate seasonal means from daily values for every cell. Now I kinda see, that it was unfortunate for me to combine the reference model data with the scenario model data, which it belongs to, in one raster file. Selecting the raster data for different time spans through giving the RasterBrick a time index with setZ creates 2 RasterStacks. It takes ages to process the function, which computes seasonal means for the different time spans for RasterStacks. Transforming the RasterStacks back to RasterBricks takes also too long. Is there any advice, how I can be more time efficient with the original RasterBrick I'm starting from? Is there the possibility to force subset not to create a RasterStack, so after using subset a RasterBrick still remains?

The grouping function was taken from here. It's grouping dates, which I've given to the raster file through setZ, into the meteorological season they belong to. Be aware, that using that function together with zApply places the layer names of the raster object (DJF, MAM, JJA, SON) alphabetically. This means, that MAM is the layer name of JJA for instance, which is wrong. Remember that.

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"

dates.reference <- seq(as.Date("1976-1-1"), as.Date("2005-12-31"), by="day")
dates.scenario <- seq(as.Date("2071-1-1"), as.Date("2100-12-31"), by="day")

dates.combined <- c(dates.reference, dates.scenario)

r.brick.reference_scenario <- setZ(r.brick, dates.combined)

#  Creates Rasterstacks, which zApply takes ages for
scenario.2071.2100.rs <- subset(r.brick.reference_scenario, which(getZ(r.brick.reference_scenario)>="2071-1-1" & getZ(r.brick.reference_scenario)<="2100-12-31")) 
reference.1976.2005.rs <- subset(r.brick.reference_scenario, which(getZ(r.brick.reference_scenario)>="1976-1-1" & getZ(r.brick.reference_scenario)<="2005-12-31")) 

scenario.2071.2100.rs <- setZ(scenario.2071.2100.rs, dates.scenario)
reference.1976.2005.rs <- setZ(reference.1976.2005.rs, dates.reference)

scenario.2071.2100.seas_mean.rs <- zApply(scenario.2071.2100.rs, by = groups(dates.scenario), fun = mean)
reference.1976.2005.seas_mean.rs <- zApply(reference.1976.2005.rs, by = groups(dates.reference), fun = mean)

# Creates RasterBricks from the RasterStacks, but the transformation takes too long, while the computing of the seasonal means from RasterBricks is time efficient
scenario.2071.2100.rs <- subset(r.brick.reference_scenario, which(getZ(r.brick.reference_scenario)>="2071-1-1" & getZ(r.brick.reference_scenario)<="2100-12-31")) 
reference.1976.2005.rs <- subset(r.brick.reference_scenario, which(getZ(r.brick.reference_scenario)>="1976-1-1" & getZ(r.brick.reference_scenario)<="2005-12-31")) 

scenario.2071.2100.rb <- brick(scenario.2071.2100.rs)
reference.1976.2005.rb <- brick(reference.1976.2005.rs)

scenario.2071.2100.rb <- setZ(scenario.2071.2100.rb, dates.scenario)
reference.1976.2005.rb <- setZ(reference.1976.2005.rb, dates.reference)

scenario.2071.2100.seas_mean.rb <- zApply(scenario.2071.2100.rb, by = groups(dates.scenario), fun = mean)
reference.1976.2005.seas_mean.rb <- zApply(reference.1976.2005.rb, by = groups(dates.reference), fun = mean)

The only option I see is to create files from the subsets and reading them in again with brick, which should be faster, but that's kinda silly.

2 Answers 2


First, you can just index a stack following something like:

Create a data string and then an index representing a query of month ranges. We then use grep to create an index to subset the raster stack to the query. Using which with %in% can be an alternative approach in some cases.

d <- seq(as.Date("2000-1-1"), as.Date("2001-12-31"), by="day")
( m <- grep(paste(c("May","June","July"), collapse ="|"), months(d)) )

Here is where you subset the raster stack, in a function or to a new stack.


You can also use this type of index directly in a function to reference specific rasters in a stack/brick.

raster::calc(r[[m]], mean)

You are really doing a lot of unnecessary gymnastic here. Since your reference data is not convolved with the scenario data, all you have to do is index the specific range of rasters to a new object.

ref <- r.brick[[1:10958]]

In some cases raster bricks can be more efficient but, you are loosing any gain in coercion. If you really want a brick object, why not just read the data as a brick. Just use the raster::brick function rather than raster::stack

Now that you have an idea on indexing a stack to create a subset I will aim you to the rts package which provides functions for analysis of time-series rasters. It is an extension of the raster package and allows you to have dates as your raster names. This facilitates functions that allow statistical summaries for day, months, years, quarterlies or defined time periods. See functions such as: apply.daily, apply.weekly, apply.yearly, period.apply, ...

The rts::rts function allows for coercion of a raster stack to a rst object. All you need for the coercion is a raster stack and the corresponding date string. The original raster stack object is just held in a slot as well as in any resulting objects. As such, the actual raster stack object can be accessed using x@raster which is necessary for writing out results.

  • Indexing like that is really convenient. I'm used to using far more code for indexing raster files through a time index. Thank you! I'm using brick and I've tried ref <- r.brick[[1:10958]] beforehand, but ref becomes a RasterStack and transforming it back into a RasterBrick through raster::brick(ref) took me 3 hours. When I have an indexed RasterBrick with daily data for a season, it takes me like 5 seconds to compute the 30 years mean for every cell in the raster object. Doing the same with a RasterStack takes hours for me.
    – Nucore
    Jun 20, 2018 at 16:51
  • Have you simply tried: ref <- brick(r.brick[[1:10958]]) BTW, you should not be seeing these type of processing time difference between stack and brick. I would contact the package maintainer over this. Jun 20, 2018 at 16:53
  • Yep, I've tried this too. Makes no difference in speed unfortunately.
    – Nucore
    Jun 20, 2018 at 17:02
  • Honestly, when I have these types of dimensions and ample memory, I coerce into a SpatialPixelsDataFrame and start using apply type functions. You can read directly to this class using rgdal::readGDAL By calling the data.frame directly from the @data slot you can start using tidyverse type functions at well. The ordering is not broken and can just be related directly back to the sp spatial object. Jun 20, 2018 at 17:55
  • Using stackApply(rs.JJA, indices = nlayers(rs.JJA), fun = "mean") to create the 30 years mean for every cell through daily summer data took 1 hour. Processing all of my data would take me 6 days of computing with this method. I need to look into rgdal, when I want to process raster data faster I guess. Basically what I want to do is the following: I want to compute the changes between a model reference time period and model scenario time period expressed as percentages. I want that for every cell, so I can create a pixel plot. Every pixel describes a change expressed as a percentage.
    – Nucore
    Jun 20, 2018 at 18:38

For time efficiency, I just transformed my brick into a matrix in the end, like that:

mat <- rasterToPoints(r.brick) 

For choosing specific seasonal data, like suggested from Jeffrey, I just created index vectors through a date vector.

Dates for the time span of interest:

dates.reference <- seq(as.Date("1976-1-1"), as.Date("2005-12-31"), by="day")

Index vector, which points at the values for the Spring season:

refidx.MAM <- grep(paste(c("March","April","May"), collapse ="|"), months(dates.reference))

reference.mat <- mat[, 3:10960]  % excludes the coordinates
reference.mat.MAM <- reference.mat[, refidx.MAM]

Do whatever you want with this seasonal data row wise (For instance I've calculated the 30yr means for every row).

I've needed my data as a raster again, so I did the following:

reference.mat.MAM <- cbind(mat[, 1:2], reference.mat.MAM)
reference.raster.MAM <- rasterFromXYZ(reference.mat.MAM, crs = crs(r.brick))

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.