[R-sig-Geo] raster: stackApply problems..

[Mark R Payne]

Hi,

[reposted from r-help]

I am trying to use stackApply() to perform averages over subsets of a
brick. However, I am struggling with the indices argument, and how it
should be interpreted. Here is a simple working example illustrating my
problem:

r <- raster()
r[] <- 1

inp <- brick(r,r,r,r,r,r)*(1:6)
res <- stackApply(inp,c(2,2,3,3,1,1),mean)

Now if we look at the values of each object:

> inp
names       : layer.1, layer.2, layer.3, layer.4, layer.5, layer.6
min values  :       1,       2,       3,       4,       5,       6
max values  :       1,       2,       3,       4,       5,       6

> res
names       : layer.1, layer.2, layer.3
min values  :     3.5,     5.5,     1.5
max values  :     3.5,     5.5,     1.5

Now, the problem is that the names and order of the layers in "res" don't
line up with the indices that I provided. You can do the maths in your head
- e.g. the first two layers of "inp" have values of 1 and 2, so their mean
should be 1.5 - however, this is ending up as layer 3 in "res".

So how should the indices argument be interpreted in this context?

Suggestion: A more intuitive format for the "indices" argument in
stackApply might be as a factor - this way the order is implict and
stackApply ends up working similar to split() or tapply()...

Best wishes,

Mark

	[[alternative HTML version deleted]]