[GRASSLIST:10635] Re: r.mapcalc problem

[Glynn Clements]

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.

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