9

I have what I thought would be a simple problem, but I haven't been able to find an appropriate answer. I have a multidimensional array v[x,y,z] and I would like to apply a function to the array along the z dimension using a grouping variable (group). Here is an example (in R):

v<-1:81
dim(v)<-c(3,3,9)
group<-c('a','a','a','b','b','b','c','c','c')

Given that the grouping variable has 3 levels (a, b and c), the result (out) I'm looking for is an array of dimension 3x3x3. I can obtain out using the following code for the above example:

out1<-apply(v[,,c(1:3)],c(1,2),sum)
out2<-apply(v[,,c(4:6)],c(1,2),sum)
out3<-apply(v[,,c(7:9)],c(1,2),sum)

library(abind)
out<-abind(out1, out2, out3, along=3) 

My question is if there is a a general means of obtaining the above result, which can be applied to large dimensional arrays and long grouping vectors.

2
  • 1
    Given your comment below @krlmlr 's answer it would be better if you could more accurately describe the data you are working with. It is frustrating to develop an answer to be told that actually, the data you have is very different from what you describe so will not work! Apr 21, 2013 at 21:55
  • 2
    Since you are dealing with remotely sensed data it makes sense to me for you to look at the raster package, and the stack and calc functions which are optimised for this kind of data. Apr 21, 2013 at 21:57

4 Answers 4

8

Easy:

out <- apply(v, c(1, 2), by, group, sum)

But to get the data in exactly the same order as you want:

out <- aperm(apply(v, c(1, 2), by, group, sum), c(2, 3, 1))
1
  • It is a function, look it up ?by.
    – flodel
    Apr 21, 2013 at 21:59
5

Using the package raster might be more appropriate for your needs. It has some code optimised for handling remotely sensed data, taking care of processing in chunks. Consider this example:

## Make 12 rasters, maybe one for each month of the year
for( i in seq(12) ){
    assign( paste0( "r" , i ) , raster( matrix(runif(1e3) , nrow = 1e2 ) ) )
}

## Create a raster stack from these
rS <- stack( mget( paste0("r",1:12) , envir = .GlobalEnv ) )

## Use calc to get mean, using by to group by a variable
## In this example I use the vector (1,1,1,2,2,2,3,3,3,4,4,4)
## meaning I get means for the first 3 rasters, then the next 3 etc
## So I get a mean for each quarter
rMean <- calc( rS , fun = function(x){ by(x , c( rep( 1:4 , each=3 ) ) , mean ) }  )

Which returns a raster brick with 4 layers (one mean for each quarter):

class       : RasterBrick 
dimensions  : 100, 10, 1000, 4  (nrow, ncol, ncell, nlayers)
resolution  : 0.1, 0.01  (x, y)
extent      : 0, 1, 0, 1  (xmin, xmax, ymin, ymax)
coord. ref. : NA 
data source : in memory
names       :         X1,         X2,         X3,         X4 
min values  : 0.02096586, 0.04015260, 0.04704145, 0.05884161 
max values  :  0.9727491,  0.9303025,  0.9804486,  0.9934670

I hope you can adapt this to your data.

4
  • Thanks a lot Simon, that looks interesting for use with raster data. I'm going to test it with a couple of data sets I have!
    – Arhopala
    Apr 21, 2013 at 23:24
  • @Arhopala how did the testing go? Did this solution work for you? Or do you need it to be more efficient/faster? Apr 23, 2013 at 13:07
  • Hi Simon, sorry for being a bit slow in replying, but was busy submitting a manuscript. I've had time to try your suggestion and it works fine. It's a bit slower than Flodel's suggestion however. I ran a system.time on both and user time for Flodel's was 14.476 and for yours 19.086 on 269 100x100 matrices. At present I'm converting satellite images to HDF5 or NetCDF, then importing into R for further analysis. Thanks for your help.
    – Arhopala
    Apr 24, 2013 at 0:56
  • The above times were in secs with a Intel® Core™ i7-2600K CPU @ 3.40GHz × 8 processor.
    – Arhopala
    Apr 24, 2013 at 11:41
2

This is much easier if your data is formatted as data frame:

library(plyr)
vd <- adply(v, 1:3)
head(vd)

  X1 X2 X3 V1
1  1  1  1  1
2  2  1  1  2
3  3  1  1  3
4  1  2  1  4
5  2  2  1  5
6  3  2  1  6

Then, you can simply attach your grouping...

vd$group <- rep(group, rep(3 * 3, length(group)))

...and split according to this grouping:

daply(vd, .(group), function(df) { ... } )

The anonymous function { ... } will be called once for each group, with df containing the sub-dataframe corresponding to that group. Here you could recombine and aggregate the data into a matrix using similar machinery. The function should return an array of dimensions 3x3x1, these will be concatenated by daply to form the desired result.

1
  • Thanks Krlmlr for your time. Unfortunately, that doesn't solve my problem. The example I gave was a using a small 'model'. The arrays I work with are very large remotely sensed data whereby each matrix in the array can represent a 1000 x 1000 point spatial matrix times the number of days (the z dimension), which consists of multiple years. I need to find the mean of each point in the 1000 x 1000 point spatial matrix for each month. In subsequent analyses I also need to maintain the array structure of my data. Thanks again for your time.
    – Arhopala
    Apr 21, 2013 at 21:27
0

split by the groups and go over the groups with lapply. Use the index to subset the array and use sum in apply. Simplify the list to an array with simplify2array.

x <- simplify2array( lapply(split(seq_along(group), group), \(i)
                    apply(v[,,i], 1:2, sum)) )
all.equal(x, out, check.attributes = FALSE)
#[1] TRUE

In this case rowSums could also be used.

x <- simplify2array( lapply(split(seq_along(group), group), \(i)
                    rowSums(v[,,i], dim=2)) )

Another way would be to use tapply inside apply where the order of the dimensions need to be reordered with aperm

x <- apply(v, 1:2, tapply, group, sum)
all.equal(aperm(x, c(2,3,1)), out, check.attributes = FALSE)
#[1] TRUE

Your Answer

Reminder: Answers generated by Artificial Intelligence tools are not allowed on Stack Overflow. Learn more

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.