[GRASS-user] What's wrong with this r.mapcalc usage?

[Glynn Clements]

Tom Russo wrote:

> The other day I tried to perform an operation with r.mapcalc that didn't
> work the way I expected.
> 
> I have a map, SAO that is 1 in some places and null elsewhere.  It represents
> an area of interest.
> 
> I have two other maps, call 'em A and B, that are elevation models over
> two areas overlapped by SAO's non-null elements.  They are null outside the
> areas where they have valid elevation, and they represent adjoining areas.
> I wanted to select all points covered by non-null pixels of SAO that have
> elevation over 11,000 feet.  The first thing I tried was this input
> to r.mapcalc:
> 
>   if(!isnull(SAO),
>      if(!isnull(A)&&A>=11000,
>         1,
>         if(!isnull(B)&&B>=11000,
>            1,
>            null()
>           )
>         ),
>       null()
>      )
> 
> (Whitespace added here just so I can be sure to match parens as I type it
> in from memory)  Here, I'm expecting that anywhere that A is null, the third
> argument of the "if(!isnull(A)[...],1,[...] would be evaluated, and that
> it would result in the same thing as if I'd done an r.patch on A and B to
> get C, then just tested like 
> "if(!isnull(SAO),if(!isnull(C)&&C>=11000,1,null()),null())".    I thought I
> would be saving the time and space of doing an r.patch.
> 
> What I got, however, was as if I'd left out the if involving B --- I got
> all those pixels overlapping SAO and A where elevation was >=11000 feet, and
> nulls where SAO overlapped B.
> 
> (Yes, g.region was set so that all of the non-null extent of SAO was included,
> and yes, there were pixels that should have matched the criteria)
> 
> I tried a bunch of variations on this, never getting the result I wanted.  In 
> the end, I just did two separate runs of r.mapcalc:
>   m1=SAO*if(!isnull(A)&&A>=11000,1,null())
>   m2=SAO*if(!isnull(B)&&B>=11000,1,null())
> 
> and an r.patch, which saved me no effort at all, even counting the time wasted
> experimenting.
> 
> Can someone spot my blindingly obvious mistake or misconception on how
> r.mapcalc's 3-argument if function is supposed to work?  I can't.

It's not so much an issue with the way that if() works (although that
is a factor), but mainly with the way that && and || work.

If A is null, then !isnull(A) is false and A>=11000 is null (not
false). The result of "false && null" is null (not false[1]), and the
result of if(null,...) is always null.

[1] This is presumably the part that caught you out.

The behaviour of && and || is consistent with most other infix
operators operators, insofar as they return null if either argument is
null.

In newer versions of r.mapcalc (I don't recall exactly which
versions), you can use the &&& and ||| operators, which behave like &&
and || except that they follow the common boolean equivalences:

	x     &&& false = false
	false &&& x     = false

	x     ||| true  = true
	true  ||| x     = true

even when x is null.

Alternatively, you can use nested if() statements, e.g. replace:

	if (!isnull(A)&&A>=11000,<true-case>,<false-case>)
with:
	if (!isnull(A),if(A>=11000,<true-case>,<false-case>),<false-case>)

Finally, even with these changes, your r.mapcalc expression isn't
quite the same as the r.patch+r.mapcalc combination, as you're
treating A<11000 the same as null, whereas r.patch only cares about
null/non-null. So, if both A and B are non-null, A<11000 and B>=11000,
you'll get null from r.mapcalc but 1 from r.patch+r.mapcalc. OTOH, if
A and B are disjoint (I'm not sure from your description), then it
doesn't make any difference.

-- 
Glynn Clements <glynn at gclements.plus.com>