[GRASS5] Strange r.mapcalc behaviour - bug?

[Markus Neteler]

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