[GRASSLIST:3347] r.mapcalc problem

Emiliano_Giovanni_Va · May 4, 2004, 10:34pm

-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

Hi all,
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.
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.
Thank you in advance, I hope my post wasn't too long.

PS: I will paste here the r.report outputs of the three mentioned maps
and ./configure options I used to compile grass-5.3 on my Gentoo box.

*~~~~~~~~~~~~~~~~*

GRASS:~ > r.report mwaybuf units=cell
r.stats: 100%
+-----------------------------------------------------------------------------+
| RASTER MAP CATEGORY REPORT
~ |
|LOCATION: leics Fri Apr 30
22:07:57 2004|
|-----------------------------------------------------------------------------|
| north: 322000 east: 456000
~ |
|REGION south: 310000 west: 444000
~ |
| res: 50 res: 50
~ |
|-----------------------------------------------------------------------------|
|MASK:none
~ |
|-----------------------------------------------------------------------------|
|MAP: M1 Motorway Buffer Zone (mwaybuf in PERMANENT)
~ |
|-----------------------------------------------------------------------------|
| Category Information
| cell|
|#|description
|count|
|-----------------------------------------------------------------------------|
|1|1000m_around_m1. . . . . . . . . . . . . . . . . . . . . . . . . . .
|10237|
|*|no data. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
|47363|
|-----------------------------------------------------------------------------|
|TOTAL
|57600|
+-----------------------------------------------------------------------------+
GRASS:~ >

*~~~~~~~~~~~~~~~~*

GRASS:~ > r.report railbuf units=cell
r.stats: 100%
+-----------------------------------------------------------------------------+
| RASTER MAP CATEGORY REPORT
~ |
|LOCATION: leics Fri Apr 30
22:08:40 2004|
|-----------------------------------------------------------------------------|
| north: 322000 east: 456000
~ |
|REGION south: 310000 west: 444000
~ |
| res: 50 res: 50
~ |
|-----------------------------------------------------------------------------|
|MASK:none
~ |
|-----------------------------------------------------------------------------|
|MAP: Rail Buffer Zone (railbuf in PERMANENT)
~ |
|-----------------------------------------------------------------------------|
| Category Information
| cell|
|#|description
|count|
|-----------------------------------------------------------------------------|
|2|1000m_around_rail. . . . . . . . . . . . . . . . . . . . . . . . . .
|10347|
|*|no data. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
|47253|
|-----------------------------------------------------------------------------|
|TOTAL
|57600|
+-----------------------------------------------------------------------------+
GRASS:~ >

*~~~~~~~~~~~~~~~~*

GRASS:~ > r.report buffers units=cell
r.stats: 100%
+-----------------------------------------------------------------------------+
| RASTER MAP CATEGORY REPORT
~ |
|LOCATION: leics Fri Apr 30
22:13:17 2004|
|-----------------------------------------------------------------------------|
| north: 322000 east: 456000
~ |
|REGION south: 310000 west: 444000
~ |
| res: 50 res: 50
~ |
|-----------------------------------------------------------------------------|
|MASK:none
~ |
|-----------------------------------------------------------------------------|
|MAP: (untitled) (buffers in PERMANENT)
~ |
|-----------------------------------------------------------------------------|
| Category Information
| cell|
|#|description
|count|
|-----------------------------------------------------------------------------|
|*|no data. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
|57600|
|-----------------------------------------------------------------------------|
|TOTAL
|57600|
+-----------------------------------------------------------------------------+
GRASS:~ >

*~~~~~~~~~~~~~~~~*

./configure --host=i686-pc-linux-gnu --prefix=/opt --with-blas
- --with-lapack --without-postgres --with-tctk --with-gd --without-odbc
- --with-motif --with-motif-includes=/usr/X11R6/include --with-freetype
- --with-nls --with-x --without-gdal

- --

________________________________________________
Emiliano Giovanni Vavassori
Allievo Ordinario - Settore di Scienze Agrarie
Scuola Superiore di Studi Universitari e di Perfezionamento
"Sant'Anna"

Homepage: http://syntaxerrormmm.altervista.org/
Chiave pubblica GPG:
http://syntaxerrormmm.altervista.org/dwls/personal/testina.asc

GuIT --- Gruppo Utilizzatori Italiani TeX
~ --- Italian Official TUG ---
Home page: http://www.guit.sssup.it/

Per favore, non inviatemi allegati in formato Word o PowerPoint.
Si veda: http://www.fsf.org/philosophy/no-word-attachments.it.html

-----BEGIN PGP SIGNATURE-----
Version: GnuPG v1.2.4 (GNU/Linux)
Comment: Using GnuPG with Thunderbird - http://enigmail.mozdev.org

iD8DBQFAmBp094fVpJPtJ6oRApePAJ4gJgcV+MhfLhwZam+3qPYbbth/xACfWCZm
S4lX4R2pnjcT/DDsVmUKqPM=
=w7jX
-----END PGP SIGNATURE-----

H_B · May 5, 2004, 12:38am

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>

...

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.

...

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.

[here nan means NULL]

1 + nan = nan
1 - nan = nan
1 * nan = nan
1 / nan = nan
if(nan == nan) = nan

etc.

The exception is isnull(nan) = 1

You either need to use r.patch to merge the two maps, or convert the
NULL values to "0" with the 'r.null null=0' command if you want a
cumulative merge. (r.patch will take its value from the first map
listed)

See also r.series

Hamish

Paul_Kelly · May 5, 2004, 1:08am

Hello

On Wed, 5 May 2004, Emiliano Giovanni Vavassori wrote:

-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

Hi all,
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

Well the tutorial clearly states that it is designed for GRASS 4.2

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'.

If either of the maps have 'no data' cells then the corresponding output
cell will also be 'no data'. So you will end up with a map with no data.
GRASS 4.x did not support null cells and just used zero. So you could add
the two maps like in the tutorial and it worked out fine.

I think you should try joining the maps using r.patch instead. Or convert
the null cells to 0 with
r.null mwaybuf null=0
r.null railbuf null=0
for both the maps and then the tutorial will work as designed with GRASS 5.x

Paul K

Glynn_Clements · May 5, 2004, 1:20am

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@virgin.net>