[GRASS5] Strange r.mapcalc behaviour - bug?

GRASS GRASS developer
1

Hi,

there seems to be a problem with / in r.mapcalc:

g.region rast=modis_vi250m_16_days_NDVI -p
projection: 99 (Transverse Mercator)
zone: 0
datum: rome40
ellipsoid: international
north: 5232499.37200598
south: 5050874.37200598
west: 1593200.7883525
east: 1801075.7883525
nsres: 125
ewres: 125
rows: 1453
cols: 1663

r.support -r modis_vi250m_16_days_NDVI
  Updating the stats for [modis_vi250m_16_days_NDVI]
   Updating the number of categories for [modis_vi250m_16_days_NDVI]

r.info -r modis_vi250m_16_days_NDVI
min=0
max=36748

#now I want to get the typical range for NDVI: -1.0 .. 1.0
#values above will be filtered away later:

r.mapcalc "modis_vi250m_16_days_NDVI.2=1. * modis_vi250m_16_days_NDVI / 10000."

r.support -r modis_vi250m_16_days_NDVI.2
  Updating the stats for [modis_vi250m_16_days_NDVI.2]
   Updating the number of categories for [modis_vi250m_16_days_NDVI.2]

r.info -r modis_vi250m_16_days_NDVI.2
min=0.000000
max=4.000000

#############
The same happens when not using r.support.

It seems to round to integer internally? Or I am missing something else.

Markus

2

Markus Neteler wrote:

there seems to be a problem with / in r.mapcalc:

g.region rast=modis_vi250m_16_days_NDVI -p
projection: 99 (Transverse Mercator)
zone: 0
datum: rome40
ellipsoid: international
north: 5232499.37200598
south: 5050874.37200598
west: 1593200.7883525
east: 1801075.7883525
nsres: 125
ewres: 125
rows: 1453
cols: 1663

r.support -r modis_vi250m_16_days_NDVI
  Updating the stats for [modis_vi250m_16_days_NDVI]
   Updating the number of categories for [modis_vi250m_16_days_NDVI]

r.info -r modis_vi250m_16_days_NDVI
min=0
max=36748

#now I want to get the typical range for NDVI: -1.0 .. 1.0
#values above will be filtered away later:

r.mapcalc "modis_vi250m_16_days_NDVI.2=1. * modis_vi250m_16_days_NDVI / 10000."

r.support -r modis_vi250m_16_days_NDVI.2
  Updating the stats for [modis_vi250m_16_days_NDVI.2]
   Updating the number of categories for [modis_vi250m_16_days_NDVI.2]

r.info -r modis_vi250m_16_days_NDVI.2
min=0.000000
max=4.000000

#############
The same happens when not using r.support.

It seems to round to integer internally? Or I am missing something else.

I can reproduce this, but if you examine the actual cell values you
will notice that r.mapcalc isn't rounding anything.

The range calculation performed by r.support involves iterating over
the map's categories:

    if(data_type == CELL_TYPE)
        G_init_range (&range);
    else
  G_init_fp_range (&fprange);

    i = G_get_histogram_num (&histogram);
    while (i >= 0){
  if(data_type == CELL_TYPE)
      G_update_range (G_get_histogram_cat(i--,&histogram), &range);
  else
      G_update_fp_range (G_get_histogram_cat(i--,&histogram), &fprange);
    }
    if(data_type == CELL_TYPE)
        G_write_range (name, &range);
    else
        G_write_fp_range (name, &fprange);

For FP maps, the categories are determined by the map's quantisation
rules. r.mapcalc doesn't explicitly set any quantisation rules, but it
appears that libgis creates a default f_quant file containing the
string "round", which presumably indicates rounding to the nearest
integer.

Presumably iterating over categories was an optimisation which was
sensible when GRASS was restricted to integer maps, but isn't sensible
for FP maps.

Comparing r.support to the 4.3 version:

    G_init_range (&range);
    i = G_get_histogram_num (&histogram);
    while (i >= 0)
  G_update_range (G_get_histogram_cat(i--,&histogram), &range);
    G_write_range (name, &range);

It appears to have been extended to handle FP maps by "brute force",
i.e. retaining the original algorithm, but with the individual
operations extended to handle FP maps, when it should have used a
completely different algorithm for FP maps.

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

3

On Mon, Mar 22, 2004 at 09:26:17PM +0000, Glynn Clements wrote:

Markus Neteler wrote:

> there seems to be a problem with / in r.mapcalc:
>
> g.region rast=modis_vi250m_16_days_NDVI -p
> projection: 99 (Transverse Mercator)
> zone: 0
> datum: rome40
> ellipsoid: international
> north: 5232499.37200598
> south: 5050874.37200598
> west: 1593200.7883525
> east: 1801075.7883525
> nsres: 125
> ewres: 125
> rows: 1453
> cols: 1663
>
> r.support -r modis_vi250m_16_days_NDVI
> Updating the stats for [modis_vi250m_16_days_NDVI]
> Updating the number of categories for [modis_vi250m_16_days_NDVI]
>
> r.info -r modis_vi250m_16_days_NDVI
> min=0
> max=36748
>
> #now I want to get the typical range for NDVI: -1.0 .. 1.0
> #values above will be filtered away later:
>
> r.mapcalc "modis_vi250m_16_days_NDVI.2=1. * modis_vi250m_16_days_NDVI / 10000."
>
> r.support -r modis_vi250m_16_days_NDVI.2
> Updating the stats for [modis_vi250m_16_days_NDVI.2]
> Updating the number of categories for [modis_vi250m_16_days_NDVI.2]
>
> r.info -r modis_vi250m_16_days_NDVI.2
> min=0.000000
> max=4.000000
>
>
> #############
> The same happens when not using r.support.
>
> It seems to round to integer internally? Or I am missing something else.

I can reproduce this, but if you examine the actual cell values you
will notice that r.mapcalc isn't rounding anything.

The range calculation performed by r.support involves iterating over
the map's categories:

    if(data_type == CELL_TYPE)
        G_init_range (&range);
    else
  G_init_fp_range (&fprange);

    i = G_get_histogram_num (&histogram);
    while (i >= 0){
  if(data_type == CELL_TYPE)
      G_update_range (G_get_histogram_cat(i--,&histogram), &range);
  else
      G_update_fp_range (G_get_histogram_cat(i--,&histogram), &fprange);
    }
    if(data_type == CELL_TYPE)
        G_write_range (name, &range);
    else
        G_write_fp_range (name, &fprange);

For FP maps, the categories are determined by the map's quantisation
rules. r.mapcalc doesn't explicitly set any quantisation rules, but it
appears that libgis creates a default f_quant file containing the
string "round", which presumably indicates rounding to the nearest
integer.

Presumably iterating over categories was an optimisation which was
sensible when GRASS was restricted to integer maps, but isn't sensible
for FP maps.

Comparing r.support to the 4.3 version:

    G_init_range (&range);
    i = G_get_histogram_num (&histogram);
    while (i >= 0)
  G_update_range (G_get_histogram_cat(i--,&histogram), &range);
    G_write_range (name, &range);

It appears to have been extended to handle FP maps by "brute force",
i.e. retaining the original algorithm, but with the individual
operations extended to handle FP maps, when it should have used a
completely different algorithm for FP maps.

Oh, I see. Does this mean that the result might be ok
but simply r.support's results wrong?

Say, If I simply
- r.in.gdal
- r.mapcalc

will I get correct results? Then I just drop r.support from
the above procedure. Maybe it's not needed at all, but just
damaging the range file.

Markus

4

Markus Neteler wrote:

Oh, I see. Does this mean that the result might be ok
but simply r.support's results wrong?

Yes. If the original map contained values up to and including 36748,
the result map will contain values up to and including 3.6748 (to
within the accuracy of double-precision FP).

Say, If I simply
- r.in.gdal
- r.mapcalc

will I get correct results?

Yes.

Then I just drop r.support from
the above procedure. Maybe it's not needed at all, but just
damaging the range file.

I'm not sure what purpose "r.support -r" serves. The core libgis
functions automatically update the range; AFAICT, the only reasons
that the range/fp_range files may not accurately reflect the range of
the data are:

1. The actual cell/fcell data files have been modified other than
via the standard functions (open_cell/put_row/close_cell).

2. The range has been explicitly modified via e.g. G_update_range()
(which is what "r.support -r" does).

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

5

On Tue, Mar 23, 2004 at 11:29:33AM +0000, Glynn Clements wrote:

I'm not sure what purpose "r.support -r" serves. The core libgis
functions automatically update the range; AFAICT, the only reasons
that the range/fp_range files may not accurately reflect the range of
the data are:

1. The actual cell/fcell data files have been modified other than
via the standard functions (open_cell/put_row/close_cell).

2. The range has been explicitly modified via e.g. G_update_range()
(which is what "r.support -r" does).

So maybe we remove it?

Maybe it was introduced once in the inter version to overcome above
mentioned problems, then introduced into cmd to have it available on
command line.

In case of removal it should be removed from both cmd and inter.

Markus

6

Markus Neteler wrote:

> I'm not sure what purpose "r.support -r" serves. The core libgis
> functions automatically update the range; AFAICT, the only reasons
> that the range/fp_range files may not accurately reflect the range of
> the data are:
>
> 1. The actual cell/fcell data files have been modified other than
> via the standard functions (open_cell/put_row/close_cell).
>
> 2. The range has been explicitly modified via e.g. G_update_range()
> (which is what "r.support -r" does).

So maybe we remove it?

Maybe it was introduced once in the inter version to overcome above
mentioned problems, then introduced into cmd to have it available on
command line.

In case of removal it should be removed from both cmd and inter.

I suspect that r.support has its place, but for consistency with the
way in which the range is normally calculated, the calculation should
iterate over the cell values, not the category list (for both integer
and FP).

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