[GRASSLIST:10637] Re: r.mapcalc problem

[Dylan Beaudette]

On Tuesday 28 February 2006 03:24 pm, Glynn Clements wrote:
> Silvia Franceschi wrote:
> > I am trying to do a mapcalc between two maps, the first one is a dsm
> > (digital surface model) and the second is a dtm (digital terrain model).
> > The dsm is defined only where there are features on the surface and the
> > dtm is defined all over the region.
> > I would like to fill the holes in the dsm map with the correspondets
> > values of the dtm using a mapcalc:
> >
> > r.mpacalc "new_map=if(dsm,dsm,dtm)"
> >
> > Using this command I obtain a new map equal to the dsm map.
> >
> > If I use these instead:
> >
> > r.mpacalc "new_map=if(dsm!=null(),dsm,dtm)"
> > r.mpacalc "new_map=if(dsm==null(),dtm,dsm)"
> >
> > I obtain two maps with all values nan.
> >
> > Do someone knows the right way to do this mapcalc?
>
> As Ivan suggested, you need to use isnull(). I'll elaborate on this
> slightly.
>
> The behaviour of null values in r.mapcalc is similar to the way that
> floating-point arithmetic handles NaN. Almost any expression which
> involves a null value will itself evaluate to null.
>
> The easiest way to understand how null values are handled is to treat
> null as representing an "unknown" value.
>
> E.g. in the above, "dsm == null()" is interpreted as "dsm == x", where
> x is unknown. As the value of x is unknown, whether or not it is equal
> to dsm is also unknown. Hence, the expression "dsm == null()"
> evaluates to null.
>
> Similarly, "null() == null()" is interpreted as "x == y", where both x
> and y are unknown, so again the result is null.
>
> The only cases where expressions involving null values don't evaluate
> to null are isnull() and cases where the operands which are null never
> have any effect upon the result, i.e.:
>
> 1. if() expressions where the "unused" branch is null:
>
> 	if(0,null,b) == b
> 	if(1,a,null) == a
>
> 2. eval() expressions, where arguments other than the last are null:
>
> 	eval(null,null,null,x) == x
>
> In every other case, if any operand is null, the result is null.
>
> However, there is one slight flaw in the "null means unknown"
> interpretation. r.mapcalc doesn't do any arithmetic optimisations, so
> you can get a null result even in specific cases where the result
> isn't actually unknown, e.g. the following all evaluate to null:
>
> 	null * 0
> 	eval(x = null(), x - x)
> 	eval(x = null(), sin(x) ^ 2 + cos(x) ^ 2 - 1)
> 	eval(x = int(null()), y = int(null()), z = int(null()), n = int(null()), \
> 		x > 0 && y > 0 && z > 0 && n > 2 && x ^ n + y ^ n == z ^ n)
>
> Theoretically, all of these should evaluate to zero according to the
> "null means unknown" interpretation (well, the last one probably
> does). In practice, it would be trivial to make the first expression
> evaluate to zero, harder for the second, much harder for the third,
> and essentially impossible for the last.

Wow. I think that some of these little snippets of wisdom should find their 
way into some of the relevant documentation... just a matter of finding the 
time...

Cheers,

-- 
Dylan Beaudette
Soils and Biogeochemistry Graduate Group
University of California at Davis
530.754.7341