[GRASSLIST:3353] Re: r.mapcalc problem

[Glynn Clements]

Emiliano Giovanni Vavassori wrote:

> My name is Emiliano and I'm an Agricultural Sciences student from Italy.
> GIS will eventually be useful to my profession, so I am trying to learn
> how to use this tool for territorial analysis. I read the GRASS Seeds
> tutorial by ASSIST and the GRASS Tutorial by M.Lennert, but maybe
> there's something (that) I don't understand with r.mapcalc.
> I tested GRASS 5.0.2, 5.0.3 and 5.3 on my Gentoo Linux Box, compiled
> from sources and ebuild to try to understand what the problem could be.
> I have also tried on a Debian Linux box with GRASS 5.0.3 with no
> positive results.
> The problem I am experiencing with GRASS is related to simple algebra
> map operations such as '+' and '-'. On page B-13 of the GRASS Seeds
> tutorial, the subsequent step is required:
> 
> ~   r.mapcalc 'buffers = mwaybuf + railbuf' <enter>
> 
> (well, the tutorial requires interactive mode, but that doesn't
> matter---I still have the same problem). At that time, the 'mwaybuf'
> raster map contains only a category, with number 1, that represents a
> buffer zone around a motorway. 'railbuf' also represents a buffer with
> category 2, but around railways. In those cases, the rest of the maps
> are 'no data' cells. My guess is that with a map addition I should
> expect to have a 'buffers' raster map that contains all buffers, marked
> with two different colours (if displayed with d.rast) and as well, in
> r.report there should be 2 categories, but what I obtain is an empty
> 'buffers' map, where all cells are marked 'no data'.
> What also seems strange is that '-' operation on MASK (page B-16 of the
> same tutorial) results in an intersection (a logic AND) between two
> marked maps.

The "GRASS Seeds" tutorial was written for 4.3, where "no data" (null)
values were represented by zero. In 5.x, there is a distinguished null
value.

Most r.mapcalc functions, including the arithmetic operators, return
null if any of the arguments are null. So, in 5.x, the above
expression would result in a map where the only non-null cells are
those which are non-null in both the mwaybuf and railbuf maps. As the
two maps don't intersect, the entire result map is null.

You can convert null values using r.null, e.g.:

	r.null mwaybuf null=0
	r.null railbuf null=0

The maps should then resemble those which would have been created by
the version of r.buffer in 4.x, and the above r.mapcalc expression
should work.

Alternatively, you can use the following r.mapcalc expression:

	r.mapcalc 'buffers = if(isnull(mwaybuf),0,mwaybuf) + if(isnull(railbuf),0,railbuf)'

> I also read 'Performing Map Calculations on GRASS Data: r.mapcalc
> Program Tutorial' by M.Larson, M.Shapiro and S.Tweddale and 'r.mapcalc
> An Algebra for GIS and Image Processing' by M.Shapiro and J.Westervelt,
> but I did not learn anything that I did not know before.
> What I am trying to understand is whether I should submit a bug
> report or whether I'm missing something about r.mapcalc operations.

No, the problem is that the way in which "no data" cells are handled
has significantly changed since the tutorial was written. It used to
be the case that zero was interpreted as "no data" in some situations
(e.g. r.buffer would output zero for cells outside the buffer) but
interpreted as the value zero in other cases (e.g. r.mapcalc treats it
as zero).

In 5.x, r.mapcalc should consistently treat null values as "unknown",
in the sense that, for most operations, if one or more of the inputs
are unknown, the result will also be unknown. The main exceptions to
this rule are:

1. isnull() returns 1 if its argument is null and zero otherwise.

2. functions such as if() and eval() which ignore certain arguments
aren't influenced by whether the ignored arguments are null; e.g. 
if(1,a,b) will always return a, regardless of whether b is null.

OTOH, r.mapcalc doesn't check for "special cases"; e.g. if x is null,
then "x * y" will be null even if y is zero.

-- 
Glynn Clements <glynn.clements at virgin.net>