error in reasoning: logical operators in raster calculator

Question

I'm processing a formula in QGIS raster calculator and have an error in reasoning.

Having the layers a and b I want to create output layer c under the following conditions:

(1) if a = 1 or b = 1 --> c = 1
(2) if a = 5 or b = 5 --> c = 5, whereas condition (1) should overrule condition (2)
(3) if a + b = 0 --> c = 0
(4) else a + b / (number of layers which are != 0; !=1; !=5

Below is the formula with comments and hopefully readable.

It actually works, but there is a NoData Output for a+b=0.

• a=1 OR b=1 --> c=1:

  (
      ("a@1"=1) OR ("b@1"=1)
  )+ 
  

• a=5, b not equal 1; OR b=5, a not equal 1 --> c=5:

  (
      (
          (("a@1"=5) AND ("b@1"!=1)) OR 
          (("b@1"=5) AND ("a@1"!=1))
      )*5
  )+
  

• a=0 and b=0, -->c=0:

  (
      (("a@1"=0) AND ("b@1"=0)) *0
  )+
  

• a and b not equal 1; 5 --> c= average of a and b, disregarding 0 values:

  (
      (
          ("a@1"!=1)  AND 
          ("b@1"!=1)  AND 
          ("a@1"!=5)  AND 
          ("b@1"!=5)  AND
          (("a@1"!=0) AND ("b@1"!=0))
      )
   ) * 
      (("a@1"+"b@1")/
          (
              (("a@1">=2) AND ("a@1"<=4)) + 
              (("b@1">=2) AND ("b@1"<=4))
          )
      )
  )
  

What am I missing?

I have a similar question pending for SAGA Raster Calculator, but would prefer using QGIS Raster Calculator.

Answer 1

I changed the code a little. Now it works, I think.

(
    ("a@1"=1) OR ("b@1"=1)
)+ 

(
    (
        (("a@1"=5) AND ("b@1"!=1)) OR (("b@1"=5) AND ("a@1"!=1))
    )*5
)+


(
    (
        ("a@1"=0) AND ("b@1"=0))
    *0
)+


(   
    (
        (("a@1">1) AND ("a@1"<5)) AND 
         ("b@1" = 0)
     )*"a@1"
) +


(
    (
        (("b@1">1) AND ("b@1"<5)) AND 
         ("a@1" = 0)
    )*"b@1"
)+


(
    (
        (
            (
                ("a@1">1) AND 
                ("a@1"<5)) AND 
                (
                    ("b@1">1) AND 
                    ("b@1"<5)
                )
            )*
        ("a@1"+"b@1")
    )/2
)

Answer 2

With RQGIS, the code is easiest to understand, check this:

##a=raster
##b=raster
##c=output raster
library(raster)
c = stack(a)[[1]]
c[] = -9999 # Temporal NoData value
c[a == 1 | b ==1] = 1 # Rule 1
c[(a == 5 | b == 5) & c !=1] = 5 # Rule 2
c[a + b == 0] = 0 # Rule 3
r = (a + b)/((!a %in% c(0,1,5)) + (!b %in% c(0,1,5))) # Rule 4
c[c == -9999] = r[c == -9999]
c[c == -9999] = NA # Back to NoData

Applied to an example dataset:

Evaluations:

Note: rasters need to have the same extension and resolution, if not. extend() or projectRaster() functions can be applied