Skip to content

raster package - zonal - standard deviations

5 messages · Rafael Wüest, Christian Levers, Robert J. Hijmans

#
Hello,

I am using the raster package to calculate zonal statistics for my stack 
of ca. 50 rasters (each raster is about 33 mio. pixels) and I am 
wondering if there is a simple way to calculate standard deviations for 
the zones since "sd" is not implemented in the presets for the "stat" 
argument. If not, I guess defining a function which calculates std.dev 
is the easiest solution.

Thanks in advance and a nice weekend,
Christian
2 days later
#
Hi Christian

What version of raster do you use. Using my installation of raster (version 2.0-12), 'sd' as a function in zonal works fine (see example below).

r <- raster(ncols=10, nrows=10)
r[] <- runif(ncell(r)) * 1:ncell(r)
z <- r
z[] <- rep(1:5, each=20)
zonal(r, z, 'sd')
     zone      "sd"
[1,]    1  4.894058
[2,]    2  9.216353
[3,]    3 14.144341
[4,]    4 19.973913
[5,]    5 27.934773

I guess updating the raster package will solve the problem.

Best,
Rafael
On 02.11.2012, at 18:35, Christian Levers wrote:

            
--
Rafael W?est
Swiss Federal Research Institute WSL
Z?rcherstrasse 111
8903 Birmensdorf
Switzerland

+41 44 7392126
rafael.wueest at wsl.ch
http://www.wsl.ch/info/mitarbeitende/wueest/index_EN
#
Hi Rafael,

thanks for the quick reply. I'm using version 2.0-12, on a Win7 platform 
though. Since your little code works also fine for me, I guess its a 
matter of my input raster(s). Here is the printout:

 > relnrg_mskd <- zonal(slope_final_mskd.tif, zone.raster, 'sd')
Fehler in .local(x, z, ...) : stat can be 'sum', 'mean', 'min', or 'max'

f?r:
 > slope_final_mskd.tif
class       : RasterLayer
dimensions  : 5320, 6270, 33356400  (nrow, ncol, ncell)
resolution  : 1000, 1000  (x, y)
extent      : 1500000, 7770000, 740000, 6060000 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 
+ellps=GRS80 +units=m +no_defs
names       : slope_final_mskd
values      : 0, 46  (min, max)

 > zone.raster
class       : RasterLayer
dimensions  : 5320, 6270, 33356400  (nrow, ncol, ncell)
resolution  : 1000, 1000  (x, y)
extent      : 1500000, 7770000, 740000, 6060000 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 
+ellps=GRS80 +units=m +no_defs
names       : harvestmap_efi
values      : 1, 460  (min, max)

Cheers,
Christian

Am 05.11.2012 16:35, schrieb Rafael W?est:

  
    
#
Thanks a lot Robert!

Indeed it's (obviously) a matter of memory (although our server has 
approx 200GB RAM). While simulating RAM limitations, I obtain the same 
message: error in .local(x, z, ...) : stat can be 'sum', 'mean', 'min', 
or 'max'.

Changing maxmemory doesn't improve things, even with 1E+20. I was 
thinking on looping through the zones (my layer has just 460) and using 
cellStats, too.

Thanks again and all the best from Berlin,
Christian


Am 05.11.2012 19:11, schrieb Robert J. Hijmans: