[GIS] compare three rasters and three conditions with overlay in r

conditionaloverlayrraster

I'm trying to run an overlay function with rasters where I want to meet all of 3 different conditions at each cell 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. I think that to meet all three conditions I would need something equivalent to &&, but this can't be used here because it is not vectorized.

You can see this with this example below. Since the third raster has NA values and the result raster has fewer NA values, it is clear that it is not evaluating all three arguments.

Hoping someone can point me in the right direction.

library(raster)
fn <- system.file("external/test.grd", package="raster")
  s <- stack(fn, fn,fn)
  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

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

I copied some of the building code from this post Compare two rasters with raster overlay and replace cell value with variable

Best Answer

Perhaps you need to first look at how ifelse works. I get the same results when I use it "stand-alone" and within a call to raster::overlay.

a <- rep(2, 5)
b <- rep(1, 5)
d <- c(2, NA, 2, NA, 2)

library(raster)
r <- raster(nrow=1, ncol=5)
A <- setValues(r, a)
B <- setValues(r, b)
D <- setValues(r, d)
s <- stack(A,B,D)

ifelse(a==2 & b==1 & d==2, 1, 0)
#[1]  1 NA  1 NA  1
res <- overlay(s, fun = function(x,y,z) ifelse( x == 2 & y == 1 & z ==2, 1, 0) )
values(res)
#[1]  1 NA  1 NA  1

ifelse(a==1 & b==1 & d==2, 1, 0)
#[1] 0 0 0 0 0
res <- overlay(s, fun = function(x,y,z) ifelse( x == 1 & y == 1 & z ==2, 1, 0) )
values(res)
#[1] 0 0 0 0 0

This is why that happens:

FALSE & NA
#[1] FALSE
TRUE & NA
#[1] NA

So if you want NA values to trump FALSE you need to write another function. You can also do some post-processing after overlay:

msk <- calc(s, fun=function(x) any(is.na(x)))
resm <- mask(res, msk, maskvalue=TRUE)
Related Question