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

Tom_Russo · January 30, 2008, 6:07pm

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.

--
Tom Russo KM5VY SAR502 DM64ux http://www.swcp.com/~russo/
Tijeras, NM QRPL#1592 K2#398 SOC#236 AHTB#1 http://kevan.org/brain.cgi?DDTNM
"And, isn't sanity really just a one-trick pony anyway? I mean all you get is
one trick, rational thinking, but when you're good and crazy, oooh, oooh,
oooh, the sky is the limit!" --- The Tick

Glynn_Clements1 · January 31, 2008, 9:11am

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@gclements.plus.com>

neteler · January 31, 2008, 8:31am

On Jan 31, 2008 10:11 AM, Glynn Clements <glynn@gclements.plus.com> wrote:
...

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.

For the record: You kindly implemented that in May 2006:
http://lists.osgeo.org/pipermail/grass-dev/2006-May/022848.html

Markus

Tom_Russo · January 31, 2008, 2:07pm

On Thu, Jan 31, 2008 at 09:11:39AM +0000, we recorded a bogon-computron collision of the <glynn@gclements.plus.com> flavor, containing:

Tom Russo wrote:

> The other day I tried to perform an operation with r.mapcalc that didn't
> work the way I expected.
>

[...]

>
> if(!isnull(SAO),
> if(!isnull(A)&&A>=11000,
> 1,
> if(!isnull(B)&&B>=11000,
> 1,
> null()
> )
> ),
> null()
> )

[...]

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.

Yep. That's it. My C/C++ background caught me, coz in those languages
the right-hand operand of && isn't evaluated if the left-hand operand is
false.

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:

Yes, I saw that &&& and ||| existed in my version of r.mapcalc, but didn't
RTFM carefully enough to realize they were what I wanted.

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.

Yes, A and B are disjoint DEMs of adjoining USGS quadrangles.

--
Tom Russo KM5VY SAR502 DM64ux http://www.swcp.com/~russo/
Tijeras, NM QRPL#1592 K2#398 SOC#236 AHTB#1 http://kevan.org/brain.cgi?DDTNM
"And, isn't sanity really just a one-trick pony anyway? I mean all you get is
one trick, rational thinking, but when you're good and crazy, oooh, oooh,
oooh, the sky is the limit!" --- The Tick