[R-sig-Geo] Subsetting RasterBrick very slow

Oscar Perpiñan Lamigueiro oscar.perpinan at upm.es
Mon Jul 4 13:35:32 CEST 2011


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>

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'),

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:

  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 



Oscar Perpiñán Lamigueiro
Dpto. de Ingeniería Eléctrica


More information about the R-sig-Geo mailing list