[GRASS5] Strange r.mapcalc behaviour - bug?

Glynn Clements glynn.clements at virgin.net
Mon Mar 22 16:26:17 EST 2004

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);
	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);
	    G_update_fp_range (G_get_histogram_cat(i--,&histogram), &fprange);
    if(data_type == CELL_TYPE)
        G_write_range (name, &range);
        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

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>

More information about the grass-dev mailing list