-1

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))
}
10 Rep
  • 2,217
  • 7
  • 19
  • 33
  • I do not think your conclusion about why it's not working are correct. ?which.max says: "Missing and NaN values are discarded." Would you be able to provide a small subset of that object using dput()? – IRTFM Sep 17 '12 at 03:32
  • Aww. C'mon. That's a webpage, not an FTP server. You should provide the code you used to create this object! I have a solution in mind but cannot test it with this degree of ambiguity. – IRTFM Sep 17 '12 at 03:50
  • Please provide [reproducible example](http://stackoverflow.com/q/5963269/429846), one that we can run quickly to test out problems and solutions to your question. – Gavin Simpson Sep 17 '12 at 15:25
  • All, I will try and add some more information this evening. I'll pull together the different methods I have tried and there resulting errors. Thanks for the help so far. –  Sep 17 '12 at 19:31

2 Answers2

5

Here's a guess at a solution. It wasn't due to which.max not "supporting" the na.rm argument just that it already assumed it and only "had room" for a data argument. Tested on a small test case taken from the help pages but not on your data. You can use either one of these:

require(raster)
 which.max2 <- function(x, ...) which.max(x)           # helper function to absorb "na.rm"
 wsa <- stackApply(PRISM_stack, rep(1,12), fun=which.max2, na.rm=NULL)

Apparently this approach does not need the helper function to strip na.rm:

calc(PRISM_stack, which.max)

Wiht the new problem of all NA's in a cell this seems to succeed with either approach:

which.max2 <- function(x, ...) ifelse( length(x) ==sum(is.na(x) ), 0, which.max(x))

As does this:

which.max2 <- function(x, ...) ifelse( length(x) ==sum(is.na(x) ), NA, which.max(x))
IRTFM
  • 258,963
  • 21
  • 364
  • 487
  • DWin. Thanks for your work. I thought I provided enough information. A link to the data (simple enough to download) and all of the code used to load in the data. You would simple have to change the data search path. However, I will give your solution a shot this evening and share the results. I know the calc(PRISM_stack,which.max) was not working for me. I'll reproduce the error for you. –  Sep 17 '12 at 19:28
  • DWin. A side note. An issue that I was having was when all 12 layers of the RasterStack were "NA' this resulted in the character(0) return for which.max() which may be causing my issues. Again, i'll get more information together this evening. Thanks! –  Sep 17 '12 at 19:33
  • What I browsed to that page it had more than 48 links and a bunch of years and it was not clear to me what files had been used. So you were expecting us to repeat your "by hand" and not giving even explicit directions how to do it. I was expecting that you would build code like `sapply(filenames, function(x) x <- read.table(file=paste0(url, x)) )` – IRTFM Sep 17 '12 at 19:42
  • I have updated the question will all information needed to reproduce my errors. –  Sep 18 '12 at 03:21
  • Is there any benefit (efficiency?) to your which.max2() functions versus the which.max.na() function I submitted. I am still learning the R language and would love some insight here. –  Sep 18 '12 at 04:39
0

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