[R-sig-Geo] Subsetting RasterBrick very slow

[Oscar Perpiñan Lamigueiro]

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