[GRASS5] Strange r.mapcalc behaviour - bug?

[Glynn Clements]

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