0

I'm trying to run an overlay function with rasters where I want to meet all of 3 different conditions at each cell (using == and > or < operators) and produce a single raster as output.

Running ifelse with the & operator seems to look at the conditions in a linear fashion from left to right - If the first two conditions are met then it will produce the if condition as output, regardless of the third condition. && can't be used here because the result is not vectorized.

You can see this with this example below where with the resulting raster it's clear that it is not evaluating all three arguments. First clue is that it produces results even if some values are NA in the third raster.

I think I can get the result I want by first checking the condition of x and y and then with this result separately checking the condition of z with a different function, but I was hoping to be able to do it all in one function (seems like this should be possible, at least).

Hoping someone can point me in the right direction.

library(raster)
fn <- system.file("external/test.grd", package="raster")
  s <- stack(fn, fn,fn)
  #Create grids
  s[[1]] <- round(runif(ncell(s), 1, 2))
  s[[2]] <- round(runif(ncell(s), 1, 2))
  s[[3]] <- round(runif(ncell(s), 1, 2))
  #convert some values in s[[3]] to NA
  s[[3]][s[[3]] == 1]<- NA

  #run overlay function
  result.rast <- overlay(s[[1]], s[[2]], s[[3]], fun =   
  function(x,y,z) { 
        ifelse( x == 2 & y == 1 & z ==2, 1, 0) 
  } )
DJ77
  • 3
  • 2

2 Answers2

0

Did you try stackApply?

You can also use each layer of the raster stack as vectors. Here is an example (it might be a better way to reference the cells in the rasterStack, though)

tt <- raster(ncol=4,nrow=5)
tt[] <- 1
tts <- stack(tt,tt,tt)
tts[[1]][4,2]<-NA
# now the condition
tt2 <- (tts[[1]] == 1 & tts[[2]] == 1 & tts[[3]] == 1)
plot(tt2)
Carlos Alberto
  • 598
  • 3
  • 9
  • Thanks. I tried your suggested code above but modified it to match the rasters and conditions I was trying to meet in my code above. `tt2 <- (s[[1]] == 2 & s[[2]] == 1 & s[[3]] == 2)` So it looks to me like it is not doing anything different. From reading the stackApply help it sounds like my function needs to be able to use `na.rm`, but my function as written does not. This may be my problem, if I could use `na.rm` as part of my condition my problem may be solved... – DJ77 Jan 15 '16 at 23:27
  • what do yo mean with "it is not doing anything different"? do you still get a non NA value in the position [4,2] of the final raster? – Carlos Alberto Jan 16 '16 at 03:37
0

I do not see evidence for the third condition not being used. NA values are a special case. See function f2 for some things you can do.

It is easier to see what's going on with a smaller raster

library(raster)
set.seed(0)
r <- raster(ncol=10, nrow=10, xmn=0, xmx=10, ymn=0, ymx=10)
r1 <- setValues(r, round(runif(ncell(r), 1, 2)))
r2 <- setValues(r, round(runif(ncell(r), 1, 2)))
r3 <- setValues(r, round(runif(ncell(r), 1, 2)))
r3[r3 == 1] <- NA
s <- stack(r1, r2, r3)


res1 <- overlay(s, fun =   
  function(x,y,z) { 
        ifelse( x == 2 & y == 1 & z ==2, 1, 0) 
  } )


#A more complex function, that keeps NAs

f2 <- function(x,y,z) {
    a <- rep(0, length(x))
    a[x == 2 & y == 1 & z ==2] <- 1
    a[is.na(x) | is.na(y) | is.na(z)] <- NA
    a 
} 

res2 <- overlay(s, fun = f2)
Robert Hijmans
  • 40,301
  • 4
  • 55
  • 63
  • I should have made more clear that NA values were an important consideration for me. What I wanted to happen was for the function to not be able to compute a value (i.e. return NA) if any cell it was considering was NA. RobertH's `f2`essentially does this after running through the initial set of conditions. Thanks so much everyone. – DJ77 Jan 17 '16 at 19:53