Oscar, I find your idea ingenious, however I don't see how this method
could still work with out of memory raster (since you load the whole raster
using getValues()). data.table storage is more efficient, so maybe it's the
fact that you can shrink the space taken by the two layers that allow to
process rasters in memory ? Or is canProcessInMemory() too conservative ? I
think it's using a 10% buffer.

Etienne


2013/2/12 Oscar Perpiñán Lamigueiro <oscar.perpinan@gmail.com>

> Hello,
>
> I think that the combination of the raster and data.table packages can
> be really useful. For example, here it is a very simplified version of
> zonal using the data.table methods. It provides the same results about
> 10 times faster.
>
> Best,
>
> Oscar.
>
> --
> Oscar Perpiñán Lamigueiro
> Dpto. Ingeniería Eléctrica
> EUITI-UPM
> URL: http://procomun.wordpress.com
> Twitter: @oscarperpinan
>
> ----------------------------------------------------------------
> library(data.table)
>
> myZonal <- function (x, z, stat = "mean", digits = 0, na.rm = TRUE,
>                      ...) {
>   fun <- match.fun(stat)
>   vals <- getValues(x)
>   zones <- round(getValues(z), digits = digits)
>   rDT <- data.table(vals, z=zones)
>   setkey(rDT, z)
>   rDT[, lapply(.SD, fun), by=z]
> }
>
> ## A RasterLayer with 1e6 cells
> r <- raster(ncols=1000, nrows=1000)
> r[] <- runif(ncell(r)) * 1:ncell(r)
> z <- r
> z[] <- rep(1:10, each=ncell(r)/10)
>
> > system.time(print(zonal(r, z)))
>       zone    "mean"
>  [1,]    1  24936.45
>  [2,]    2  75138.38
>  [3,]    3 124880.63
>  [4,]    4 175018.38
>  [5,]    5 224632.18
>  [6,]    6 274356.22
>  [7,]    7 325141.16
>  [8,]    8 375625.27
>  [9,]    9 424195.16
> [10,]   10 474595.73
>    user  system elapsed
>   3.988   0.220   4.215
> > system.time(print(myZonal(r, z)))
>      z      vals
>  1:  1  24936.45
>  2:  2  75138.38
>  3:  3 124880.63
>  4:  4 175018.38
>  5:  5 224632.18
>  6:  6 274356.22
>  7:  7 325141.16
>  8:  8 375625.27
>  9:  9 424195.16
> 10: 10 474595.73
>    user  system elapsed
>   0.464   0.024   0.489
>
> ## Now with a RasterBrick
> s <- brick(r, nl=4)
> s[] <- runif(ncell(s)*4)
>
> > system.time(print(zonal(s, z)))
>       zone   layer.1   layer.2   layer.3   layer.4
>  [1,]    1 0.4986786 0.5007106 0.5001397 0.4994188
>  [2,]    2 0.4993135 0.5009089 0.5006163 0.5001419
>  [3,]    3 0.4989671 0.5013064 0.4983836 0.5004493
>  [4,]    4 0.5005992 0.5009225 0.4971442 0.4999080
>  [5,]    5 0.4998065 0.4989140 0.5000887 0.5014669
>  [6,]    6 0.5010002 0.5008529 0.4987527 0.5012223
>  [7,]    7 0.5002876 0.5006790 0.4982280 0.4999350
>  [8,]    8 0.5006809 0.4997827 0.5000287 0.4995565
>  [9,]    9 0.4993470 0.5010943 0.5005641 0.5003999
> [10,]   10 0.4994278 0.4991750 0.4978472 0.5003966
>    user  system elapsed
>  10.009   0.560  10.584
> > system.time(print(myZonal(s, z)))
>      z   layer.1   layer.2   layer.3   layer.4
>  1:  1 0.4986786 0.5007106 0.5001397 0.4994188
>  2:  2 0.4993135 0.5009089 0.5006163 0.5001419
>  3:  3 0.4989671 0.5013064 0.4983836 0.5004493
>  4:  4 0.5005992 0.5009225 0.4971442 0.4999080
>  5:  5 0.4998065 0.4989140 0.5000887 0.5014669
>  6:  6 0.5010002 0.5008529 0.4987527 0.5012223
>  7:  7 0.5002876 0.5006790 0.4982280 0.4999350
>  8:  8 0.5006809 0.4997827 0.5000287 0.4995565
>  9:  9 0.4993470 0.5010943 0.5005641 0.5003999
> 10: 10 0.4994278 0.4991750 0.4978472 0.5003966
>    user  system elapsed
>   0.652   0.284   0.939
>
>
>

	[[alternative HTML version deleted]]

