UPDATE 09/17/2012
Here is a piece of code using self-contained data that reproduces my issues:
Please keep in mind, the actual data dimensions I have are huge...
dimensions : 3105, 7025, 21812625, 12 (nrow, ncol, ncell, nlayers)
What I need is the index of the max for each row,col over the layers. All NA should return NA and multiple max copys should return the first max index (or something else, must be consistent)
# Create a test RasterStack
require(raster)
a <- raster(matrix(c(11,11,11,
NA,11,11,
11,11,13),nrow=3))
b <- raster(matrix(c(12,12,12,
NA,12,12,
40,12,13),nrow=3))
c <- raster(matrix(c(13,9,13,
NA,13,13,
13,NA,13),nrow=3))
d <- raster(matrix(c(10,10,10,
NA,10,10,
10,10,10),nrow=3))
corr_max <- raster(matrix(c(13,12,13,
NA,13,13,
40,12,13),nrow=3))
stack <- stack(a,b,c,d)
which.max2 <- function(x, ...)which.max(x)
# stackApply method
max_v_sApp <- stackApply(stack,rep(1,4),which.max2,na.rm=NULL)
# calc method
max_v_calc <- calc(stack,which.max)
Hopefully this provides enough information.
UPDATE:
This might work... testing now:
which.max2 <- function(x, ...){
max_idx <- which.max(x) # Get the max
ifelse(length(max_idx)==0,return(NA),return(max_idx))
}
So, here is the final issue I found and my solution.
The issue I was having with which.max() is how it handles a vector of all NA's
>which.max(c(NA,NA,NA))
integer(0)
>
When the stackApply() function tries to write this value to a new RasterLayer it fails. Where the function should return NA, it returns integer(0) which has length = 0.
The solution (my solution) was to write a wrapper for which.max()
which.max.na <- function(x, ...){
max_idx <- which.max(x)
ifelse(length(max_idx)==0,return(NA),return(max_idx))
}
This, implemented to my original RasterStack works fine. Thanks all for the help, and please suggest alternatives to this solution if you have them!
Thanks! Kyle