1

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.

  • Order of precedence is sufficiently arcane, and varies with language and plarform, that you ought to use parenthesis and indentation and line continuation and comments to make complex expressions as readable as possible. – Vince Feb 15 '18 at 11:41
  • hopefully it is better readable now! – Kevkev770 Feb 15 '18 at 12:18
0

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:

enter image description here

enter image description here

Evaluations:

enter image description here enter image description here enter image description here enter image description here

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

  • this looks very straight forward! What do know need Raster reclass for? – Kevkev770 Feb 15 '18 at 16:03
  • Raster reclass is the custom name of this script, you can name it as you want – aldo_tapia Feb 15 '18 at 16:04
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
)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.