[GRASSLIST:3619] Re: very low fp values problem

[Glynn Clements]

Maciek Sieczka wrote:

> I'm trying to reduce vignetting on my aerial photos and luckily I'm
> approaching some results. By the way a strange thing occured when working
> with rasters with very low floating point values.
> 
> if fixed already in 5.3 or please let me know and sorry for bugging!
> 
> the problematic raster's name is 'kor_lx'
> detailed description is at the bottom and a screendump is attached
> 
> the issue step by step:
> 
>     r.info -r kor_lx
>     min=0.000000
>     max=0.000000

Am I correct in assuming that the values in the map are all of a
magnitude of ten-to-the-minus-N where N is a "large" value (e.g. 10 or
more)?

If so, then the problem is essentially that they will end up getting
approximated to zero by the use of printf("%f") or similar.

r.info uses "%f", which only shows 6 digits. The r.univar script uses
"r.stats -1n" to get the actual cell values, and "r.univar -1" uses
"%.10f", which only shows the first 10 digits after the decimal point.

The simplest workaround would be to scale the values by a large
constant such that printf("%f") retains sufficient accuracy.

> but I wanted to know the exact values so I tried:
> 
>     r.univar kor_lx
> 
>     This module calculates univariate statistics...
>     Calculation for full image kor_lx at vignet...
>     (no base mask map set, ignoring NULL cells)
>     Reading raster map...
>     r.stats:  100%
>     Calculating statistics...
> 
> but it failed:
> 
>     Number of cells (excluding NULL cells): 3407232
>     Minimum: 0
>     Maximum: 0
>     Range:  0
>     Arithmetic mean:  0
>     Variance:  0
>     Standard deviation: 0
>     awk: cmd. line:17: (FILENAME=- FNR=3407232) fatal: division by zero
>     attempted

That is due to this line:

print "Variation coefficient:",100*(sqrt((sum2 - sum*sum/N)/N))/(sum/N),"%"

"sum" is zero, so "sum/N" is zero, so it ends up dividing by zero.

> so I tried copying the map category values to their grey scale equivalents:
> 
>     r.mapcalc 'kor256=#kor_lx'
>      100%
> 
> and strange here the strange is that the max value is not 255 as it should:
> 
>     r.info -r kor256
>     min=0
>     max=254

The maximum value is mapped to an intensity value of exactly 255, but
rounding error in the lookup probably results in a value of 254.9999,
which gets truncated to 254.

I suspect that G__interpolate_color_rule() should probably be rounding
the result to the nearest integer, rather than truncating it.

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