[GRASS5] Re: [GRASS_DE] r.mapcalc round()

[Markus Neteler]

On Wed, May 16, 2001 at 03:18:40PM +0100, Glynn Clements wrote:
> Markus Neteler wrote:
> 
> > > a) the awk script "r.univar",
> > 
> > BTW: A C implementation would be better (maybe using s.univar). Would
> > be a nice exercise, but I don't have the time.
> 
> Simply re-implementing it in C should be trivial; however, that won't
> solve the problems when the variance is swamped by the rounding error.
> 
> Unless someone is planning on implementing a more robust variance
> algorithm, we should at least add a check for negative variance in
> r.univar. BTW, does s.univar suffer from the same problem?

I have tested this:
r.mapcalc plane=1.1
range: 1.1000000238 1.1000000238    # yes, we could use %f
r.to.sites -a in=plane out=plane    # without "-a" you get #1.1 in sites
                                      list which is forbidden by definition

head $LOCATION/site_lists/plane 
 name|plane
 desc|r.to.sites -a input=plane output=plane label=No Label
 form|||%
 labels|Easting|Northing|%No Label
 3570506.25|5764093.75|%1.10000002 
 3570518.75|5764093.75|%1.10000002 
 3570531.25|5764093.75|%1.10000002 
 3570543.75|5764093.75|%1.10000002 

-> precision problem as well

s.univar plane
        number of points 53760
                    mean 1.1
      standard deviation 1.45885e-13
coefficient of variation 1.32622e-11
                skewness -1
                kurtosis -2
         mean of squares 1.21
 mean of absolute values 1.1
                 minimum 1.1
          first quartile 1.1
                  median 1.1
          third quartile 1.1
                 maximum 1.1

-> far from being perfect

Unfortunately the precision problems seem to be spreaded in GRASS.

> > > > r.mapcalc test=1.1
> > > > CREATING SUPPORT FILES FOR test
> > > > range: 1.1000000238 1.1000000238
> > > 
> > > This is to be expected; r.mapcalc also uses "%.10f", which is more
> > > precision than a float has (float will be promoted to double when
> > > passed as an unprototyped argument, as is the case with printf()).
> > 
> > O.k. I think I understand the problem now. But how to deal with it?
> 
> Just use "%f" or "%.6f"
> 
> > r.mapcalc seems to generate FCELL only. Shouldn't the range printed with
> > float precision then?
> 
> Yes. float has just under 7 digits of precision (FLT_EPSILON is around
> 1.19e-7 for IEEE single precision), so "%.6f" should do it. I suspect
> that "%.10f" was "pulled out of a hat".

O.k. I have changed that to "%.6f" in CVS.
Now I get:

r.mapcalc plane=1.1
range: 1.1 1.1

which is looking better.

and:
r.mapcalc plane=1.123456789
range: 1.123457 1.123457


> > Sorry for insisting here, but above range will
> > be expected to be wrong by the average user not familiar with precision
> > issues.
> > 
> > And r.stats: maybe the printed precision should be dependent from the
> > input map type, means more decimals for DCELL maps that FCELL maps.
> 
> Maybe; although this won't help if a DCELL map is generated using
> float data somewhere in the process.
O.k., but we should address the problem when displaying a FCELL map.
In above example I get:

r.stats -1 plane
1.1234568357
1.1234568357
1.1234568357

instead of
1.1234568
1.1234568
1.1234568

Could the RASTER_MAP_TYPE be used somehow for the "fprintf()"s?

Regards

 Markus

---------------------------------------- 
If you want to unsubscribe from GRASS Development Team mailing list write to:
minordomo at geog.uni-hannover.de with
subject 'unsubscribe grass5'