[R-sig-Geo] Alternative to zonal for large images

[Oscar Perpiñán Lamigueiro]

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