[GRASSLIST:3629] more r.mapcalc troubles

GRASS GRASS user
1

Hello Everyone!

Thanks for the support so far! Here is one more strange thing in r.mapcalc I
encountered during my low-level computing :).

I'm using Grass 5.03 installed from the binaries from the Grass www.

I made a following raster for my vignetting correction purposes:

    r.mapcalc 'kor_lx = ( 1218.01 / ( sqrt( (double(x()) - 941.0) ^2.0 +
    (double(y()) - 928.0)^2.0) ^2.0 + 1218.01^2.0) ) ^4.0'

the value range for 'kor_lx' file is (read from the color file in the 'colr'
folder):

    0.00000000000002026206628416476 - 0.000000000000454355718931865126

I checked if r.mapcalc hadles the 'kor_lx''s values properly by doing:

    r.mapcalc 'copy=kor_lx'

the data range in the 'colr' file for 'copy' was DIFFERENT (the min value is
higher than in 'kor_lx'):

    0.000000000000028123670807620238 - 0.000000000000454355718931865126

the difference is:

    0.000000000000028123670807620238 - 0.00000000000002026206628416476
    = 0.000000000000007861604523455478

how come?

The difference beetween 'kor_lx' and the 'copy' is very well visible on the
attached window dump. It's a "copy" partly drawn over the "kor_lx" (caught
in the 1/4 of drawing on Grass monitor - see the upper 1/4 of the picture).

Maciek Sieczka

(attachments)

rast_diff.tif (16 KB)

2

Maciek Sieczka wrote:

Thanks for the support so far! Here is one more strange thing in r.mapcalc I
encountered during my low-level computing :).

I'm using Grass 5.03 installed from the binaries from the Grass www.

I made a following raster for my vignetting correction purposes:

    r.mapcalc 'kor_lx = ( 1218.01 / ( sqrt( (double(x()) - 941.0) ^2.0 +
    (double(y()) - 928.0)^2.0) ^2.0 + 1218.01^2.0) ) ^4.0'

the value range for 'kor_lx' file is (read from the color file in the 'colr'
folder):

    0.00000000000002026206628416476 - 0.000000000000454355718931865126

I checked if r.mapcalc hadles the 'kor_lx''s values properly by doing:

    r.mapcalc 'copy=kor_lx'

the data range in the 'colr' file for 'copy' was DIFFERENT (the min value is
higher than in 'kor_lx'):

    0.000000000000028123670807620238 - 0.000000000000454355718931865126

FWIW, I can't reproduce that behaviour.

Having created an X/Y location with the following region:

  GRASS:~ > g.region -p
  projection: 0 (x,y)
  zone: 0
  north: 1824
  south: 0
  west: 0
  east: 1868
  nsres: 1
  ewres: 1
  rows: 1824
  cols: 1868

I ran the above r.mapcalc command, followed by "r.colors kor_lx color=grey",
and obtained the same range as yourself in the colour table:

  % 0.00000000000002026206628416476 0.000000000000454355718931865126
  0.00000000000002026206628416476:0 0.000000000000454355718931865126:255

However, if I run "r.mapcalc 'copy=kor_lx'", the copy has an identical
colour table:

  % 0.00000000000002026206628416476 0.000000000000454355718931865126
  0.00000000000002026206628416476:0 0.000000000000454355718931865126:255

Also, I get the same results with both 5.0.2 and a recent 5.3
snapshot.

BTW, if I do:

  r.mapcalc 'kor256=#kor_lx'
  r.info -r kor256

  GRASS:~ > r.mapcalc 'kor256=#kor_lx'
   100%
  GRASS:~ > r.info -r kor256
  min=0
  max=255

Presumably, you are using a different colour table; which one?

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

3

> I made a following raster for my vignetting correction purposes:
>
> r.mapcalc 'kor_lx = ( 1218.01 / ( sqrt( (double(x()) - 941.0) ^2.0
> + (double(y()) - 928.0)^2.0) ^2.0 + 1218.01^2.0) ) ^4.0'
>
> the value range for 'kor_lx' file is (read from the color file in the
> 'colr' folder):
>
> 0.00000000000002026206628416476 - 0.000000000000454355718931865126
>
> I checked if r.mapcalc hadles the 'kor_lx''s values properly by doing:
>
> r.mapcalc 'copy=kor_lx'
>
> the data range in the 'colr' file for 'copy' was DIFFERENT (the min
> value is higher than in 'kor_lx'):
>
> 0.000000000000028123670807620238 - 0.000000000000454355718931865126

FWIW, I can't reproduce that behaviour.

Having created an X/Y location with the following region:

GRASS:~ > g.region -p
projection: 0 (x,y)
zone: 0
north: 1824
south: 0
west: 0
east: 1868
nsres: 1
ewres: 1
rows: 1824
cols: 1868

I ran the above r.mapcalc command, followed by "r.colors kor_lx
color=grey",
and obtained the same range as yourself in the colour table:

% 0.00000000000002026206628416476 0.000000000000454355718931865126
0.00000000000002026206628416476:0 0.000000000000454355718931865126:255

However, if I run "r.mapcalc 'copy=kor_lx'", the copy has an identical
colour table:

% 0.00000000000002026206628416476 0.000000000000454355718931865126
0.00000000000002026206628416476:0 0.000000000000454355718931865126:255

Fixed that.
The reason was trivial - there was a mask made of the aerial photo borders
which I forgot
and I didn't know that r.mapcalc ignores the mask when no raster is given as
an argument i.e. "r.mapcalc 'raster=100'" (like in case of the formula on
the very top of the email).
So there was a difference beetween the 'kor_lx' and 'copy' because only
'copy' was masked.

Just to finish this toppic, could you only explain why it is so?

Also, I get the same results with both 5.0.2 and a recent 5.3
snapshot.

BTW, if I do:

r.mapcalc 'kor256=#kor_lx'
r.info -r kor256

GRASS:~ > r.mapcalc 'kor256=#kor_lx'
100%
GRASS:~ > r.info -r kor256
min=0
max=255

Presumably, you are using a different colour table; which one?

The reason was as above.

> Is there a way to learn about the value range etc. in case of such a low
> value raster in Grass 5.03?

I don't know. The range is actually stored as binary data, so the
information exists. However, r.info displays it using "%f", so small
values will be rounded to zero. I don't know of any other program
which will provide this information.

some solution: "r.colors raster_name color=rules" prints the value range
with a pretty good precission

for instance:

    r.colors kor_lx_bis color=rules
    Enter rules, "end" when done, "help" if you need it.
    fp: Data range is 0.0000000000000281236708076 to
0.0000000000004543557189319

Maciek

4

Maciek Sieczka wrote:

The reason was trivial - there was a mask made of the aerial photo borders
which I forgot
and I didn't know that r.mapcalc ignores the mask when no raster is given as
an argument i.e. "r.mapcalc 'raster=100'" (like in case of the formula on
the very top of the email).
So there was a difference beetween the 'kor_lx' and 'copy' because only
'copy' was masked.

Just to finish this toppic, could you only explain why it is so?

The mask is applied automatically when *reading* maps. If a mask
exists, any masked cells appear as nulls. This functionality is built
into the core functions which are used to read maps, and not
individual programs, so it applies consistently to all programs which
read raster maps (unless they explicitly use the _nomask versions,
which is rare).

The mask doesn't apply to writing maps (amongst other things, this
would make it rather awkward for programs such as r.fillnulls, r.grow2
etc which are meant to fill in null cells).

The fact that "r.mapcalc dst = src" results in dst being a masked copy
of src is largely intentional; this (or, alternatively, r.resample) is
a standard technique for creating a new map which "freezes" the
current region and the current mask into an existing map.

> > Is there a way to learn about the value range etc. in case of such a low
> > value raster in Grass 5.03?

> I don't know. The range is actually stored as binary data, so the
> information exists. However, r.info displays it using "%f", so small
> values will be rounded to zero. I don't know of any other program
> which will provide this information.

some solution: "r.colors raster_name color=rules" prints the value range
with a pretty good precission

for instance:

    r.colors kor_lx_bis color=rules
    Enter rules, "end" when done, "help" if you need it.
    fp: Data range is 0.0000000000000281236708076 to
0.0000000000004543557189319

That only works interactively. I have no idea how a script could
obtain the range if the values were small. We should consider
extending "r.info -r" to allow the use of %f/%g.

This issue does get discussed occasionally, but it's never resulted in
action. The main problem is that the issue is so widespread; lots of
code uses "%f", which falls down for very large or very small values.
OTOH, it's surprising (to me) how infrequently people actually get
bitten by it.

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