[GRASS-user] Problems creating new map with r.mapcalc

GRASS GRASS user
1

Hi!
I wonder if anyone could explain to me how to create a new map (raster) from
two pre-existing maps, in mapcalc. Because I have problems with understanding
how to do this in r.mapcalc

I have two raster datasets that I would like to combine like so; if the cell
is in map a, assign the value 7; if it is in map B, assign the correct value
from B, (i.e keep the value the cell already have in B).

I have tried to read "Performing Map Calculations on GRASS Data: r.mapcalc
Program Tutorial", by Marji Larson, Michael Shapiro and Scott Tweddale, U.S.
Army Construction Engineering Research Laboratory (December 1991) , and they
give this example:

MAPCALC> aspect.str = if(streams,26,aspect)

This should then give a new map that contains the values from "aspect", but
with the value 26 assigned to cells with "streams".

Example from my experiments:
R_test9 = if(R_VannA@larsf == 1,7, R_paramNew@larsf )

This makes a new map; but it contains only 7, were there was data in R_vannA,
and else, NULL (I guess), see otput below of r.info.

I also tried "R_test9 = if(if(R_VannA@larsf == 1),7, R_paramNew@larsf ) ; and
this gives the same (non) result. And "R_test10 = if(R_VannA@larsf, 7,
R_paramNew@larsf)"; also the same as; I.E not what I wanted. The problem is
that it never assigns a value from R_paramNew@larsf!!!

What am I missing? Or doing wrong?

Regards
Larsf

Running; SUSE 10.3, grass63 dated 12.09 from
http://download.opensuse.org/repositories/Application:/Geo/openSUSE_10.3/i586/

---------------------------------------------r.info output (example)

(3) r.info map=R_test9@larsf

+--------------------------------------------------------------------------
--+

| Layer: R_test9@larsf Date: Sun Sep 21 15:28:36 2008
| | Mapset: larsf Login of Creator: larsf
| | Location: Frosta
| | DataBase: /home/larsf/grassgis
| | Title: ( R_test9 )
| | Timestamp: none
| |
| ------------------------------------------------------------------------
|----|
|
| Type of Map: raster Number of Categories: 7
| | Data Type: CELL
| | Rows: 3360
| | Columns: 3720
| | Total Cells: 12499200
| | Projection: UTM (zone 32)
| | N: 7062800 S: 7046000 Res: 5
| | E: 596600 W: 578000 Res: 5 |
| Range of data: min = 7 max = 7
| |
|
| Data Description:
| | generated by r.mapcalc
| |
|
| Comments:
| | if(R_VannA@larsf == 1, 7, R_paramNew@larsf)
| |

+--------------------------------------------------------------------------
--+

2

Argh!
The NULL's were biting me :wink:

Did this; R_VannB@larsf = if(isnull(R_VannA@larsf), 0, 1)

And lo and behold this did work; R_test14 = if(R_VannB@larsf, 7,
R_paramNew@larsf)

And produced the wanted map!

But could someone enlighten me as to what exactly is the reason why this
worked?

regards
larsf

Søndag 21. september 2008 skreiv Lars Forseth:

Hi!
I wonder if anyone could explain to me how to create a new map (raster)
from two pre-existing maps, in mapcalc. Because I have problems with
understanding how to do this in r.mapcalc

I have two raster datasets that I would like to combine like so; if the
cell is in map a, assign the value 7; if it is in map B, assign the correct
value from B, (i.e keep the value the cell already have in B).

I have tried to read "Performing Map Calculations on GRASS Data: r.mapcalc
Program Tutorial", by Marji Larson, Michael Shapiro and Scott Tweddale,
U.S. Army Construction Engineering Research Laboratory (December 1991) ,
and they give this example:

MAPCALC> aspect.str = if(streams,26,aspect)

This should then give a new map that contains the values from "aspect", but
with the value 26 assigned to cells with "streams".

Example from my experiments:
R_test9 = if(R_VannA@larsf == 1,7, R_paramNew@larsf )

This makes a new map; but it contains only 7, were there was data in
R_vannA, and else, NULL (I guess), see otput below of r.info.

I also tried "R_test9 = if(if(R_VannA@larsf == 1),7, R_paramNew@larsf ) ;
and this gives the same (non) result. And "R_test10 = if(R_VannA@larsf, 7,
R_paramNew@larsf)"; also the same as; I.E not what I wanted. The problem is
that it never assigns a value from R_paramNew@larsf!!!

What am I missing? Or doing wrong?

Regards
Larsf

Running; SUSE 10.3, grass63 dated 12.09 from
http://download.opensuse.org/repositories/Application:/Geo/openSUSE_10.3/i5
86/

---------------------------------------------r.info output (example)

> (3) r.info map=R_test9@larsf
>
> +------------------------------------------------------------------------
>-- --+
>
> | Layer: R_test9@larsf Date: Sun Sep 21 15:28:36
> | 2008
> |
> | | Mapset: larsf Login of Creator: larsf
> | |
> | | Location: Frosta
> | |
> | | DataBase: /home/larsf/grassgis
> | |
> | | Title: ( R_test9 )
> | |
> | | Timestamp: none
> |
> | ----------------------------------------------------------------------
> |-- ----|
> |
> | Type of Map: raster Number of Categories: 7
> |
> | | Data Type: CELL
> | |
> | | Rows: 3360
> | |
> | | Columns: 3720
> | |
> | | Total Cells: 12499200
> | |
> | | Projection: UTM (zone 32)
> | |
> | | N: 7062800 S: 7046000 Res: 5
> | |
> | | E: 596600 W: 578000 Res: 5 |
> |
> | Range of data: min = 7 max = 7
> |
> | Data Description:
> | | generated by r.mapcalc
> |
> | Comments:
> | | if(R_VannA@larsf == 1, 7, R_paramNew@larsf)
>
> +------------------------------------------------------------------------
>-- --+

_______________________________________________
grass-user mailing list
grass-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-user

3

Lars Forseth wrote:

The NULL's were biting me :wink:

Did this; R_VannB@larsf = if(isnull(R_VannA@larsf), 0, 1)

And lo and behold this did work; R_test14 = if(R_VannB@larsf, 7,
R_paramNew@larsf)

And produced the wanted map!

But could someone enlighten me as to what exactly is the reason why this
worked?

Nearly all operators and functions return null if any of their
arguments are null.

The main exceptions are:

"isnull(x)" returns 0 if x is non-null and 1 if x is null.

"if(1,a,b)" returns a even if b is null.

"if(0,a,b)" returns b even if a is null.

"eval(a,b,c,...,x)" returns x even if any or all of a,b,c,... are null.

"0 &&& x" and "x &&& 0" return 0 even if x is null.

"1 ||| x" and "x ||| 1" return 1 even if x is null.

"graph(x,x1,y1,...,xn,yn)" will return null if x is null or any of the
x[i] are null, or if a "relevant" y[i] is null, but not where a y[i]
which isn't used in the calculation is null.

[The &&& and ||| operators are a relatively recent addition; the older
&& and || operators follow the usual behaviour for nulls, i.e. they
return null if either operand is null.]

The behaviour is quite intuitive if you view null as representing an
unknown quantity. This explains why both "x == null()" and "x != null()"
are always null, regardless of whether x is null or non-null.

null() is an unknown value, so it's unknown whether or not it's equal
to any particular (known) value. If x is also null, then you have two
unknown values, and it's unknown whether or not the two are equal.

In many respects, the behaviour is similar to that of NaN in
floating-point arithmetic.

With the exception of the cases listed above, r.mapcalc doesn't try to
detect "special" cases where the result can be deduced even if one of
the operands is null. E.g. if x is null, "x * 0" is null rather than
zero, "x == x" is null rather than 1, "x != x" is null rather than 0,
etc.

If you have a "boolean" map where the values are either null or 1, you
normally want to convert null to 0 before performing any other
processing. You can either replace the nulls in an existing map using
e.g. "r.null null=0", or adjust the tests within r.mapcalc, e.g.:

  R_test14 = if(not(isnull(R_VannB@larsf)), 7, R_paramNew@larsf)

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