Hi, Navigating through the code, I learn that subset first makes a conversion from a RasterBrick to a RasterStack:
> getMethod('subset', 'Raster') Method Definition: function (x, ...) { .local <- function (x, subset, drop = TRUE, ...) { x <- stack(x) subset(x, subset = subset, drop = drop, ...) } .local(x, ...) } <environment: namespace:raster> Signatures: x target "Raster" defined "Raster" With a toy example I find that the slow part is due to this conversion: mat <- array(runif(12e4), dim=c(10, 10, 1200)) b <- brick(mat) idx <- seq(from=as.Date('1910-01-01'), to=as.Date('2009-12-31'), by='month') month <- function(x)as.numeric(format(x, '%m')) idxJuly <- which(month(idx)==7) ##conversion from RasterBrick to RasterStack system.time(s <- stack(b)) user system elapsed 35.570 0.000 36.603 ##subset with the index system.time(sJuly <- subset(s, idxJuly)) user system elapsed 0.080 0.000 0.098 Another approach is to build a list of RasterLayers with *only* the layers to keep: system.time({ l <- lapply(idxJuly, function(n)raster(b, n)) sJuly2 <- do.call(stack, l) } ) user system elapsed 1.684 0.000 1.755 Anyway, if your final objective is to calculate over subsets of the RasterBrick, then you may find useful the new zApply function: ##First a time index is added to the RasterBrick b <- setZ(b, idx) system.time(meanMonthly <- zApply(b, by=month, FUN=mean)) user system elapsed 0.696 0.000 0.718 Cheers, Oscar. -------------- Oscar Perpiñán Lamigueiro Dpto. de Ingeniería Eléctrica EUITI-UPM http://procomun.wordpress.com _______________________________________________ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo