[GRASSLIST:3614] very low fp values problem

Maciek_Sieczka · June 7, 2004, 7:59pm

hello fellow grassies!

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

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@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

and when run interactively:

r.univar

    Enter base raster map (or return to skip)
    Enter 'list' for a list of existing raster files
    Hit RETURN to cancel request
    > kor_lx
    <kor_lx>
    kor_lx@vignet

    Enter cover raster map
    Enter 'list' for a list of existing raster files
    Hit RETURN to cancel request
    > kor_lx
    <kor_lx>
    kor_lx@vignet

    This module calculates univariate statistics...
    Calculation for images kor_lx@vignet masked through kor_lx@vignet...
    (ignoring NULL cells)
    Reading raster maps...
    r.stats: 100%

ERROR: Map kor_lx contains only NULL data in current region.

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

GRASS:~ > r.info kor_lx

(attachments)

Glynn_Clements · June 8, 2004, 2:02am

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

Maciek_Sieczka · June 8, 2004, 4:26pm

thanks for the detailed explanation!

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)?

yes

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

could you please explain what is that 'printf("%f")'?

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.

absolutely! rounded is closer to the actual value than the trunctated one,
isn't it?

Is there a way to learn about the value range etc. in case of such a low
value raster in Grass 5.03?
Hamish suggested I should try compiling the new version of r.univar in
src/raster/r.univar2 - would that do? (But then I have to go for the Grass
5.3, right?)

Maciek

Glynn_Clements · June 9, 2004, 3:12am

Maciek Sieczka wrote:

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

could you please explain what is that 'printf("%f")'?

printf() is a standard ANSI C library function for formatting text. It
accepts a variable number arguments, the first of which is a "format
string", which may contain "conversion specifiers" (which start with a
"%" sign). It prints (on the terminal) the format string, but with
each conversion specifier replaced with a textual representation of
the corresponding argument.

For example:

printf("You scored %d out of %d\n", score, total);

If "score" had the value 8, and "total" had the value 10, it would
print:

You scored 8 out of 10

The "%f" conversion specifier causes it to display a floating-point
value in decimal notation, to six decimal places. "%.10f" is similar,
except that it uses ten decimal places.

That's probably more information than you wanted to know. What matters
is that "r.stats -1n" (which is how r.univar gets the data) only
displays values to ten decimal places (and many other programs may
only use six, as that is the default for printf's "%f" specifier,
which is the way that most programs convert floating-point values to
text).

So, if your values are of the order of ten-to-the-minus-N, where N is
more than ten, r.stats will simply output one long stream of
"0.0000000000" values, with all of the significant digits being lost.

While printf() can display values in exponential form (e.g. 1.2e-10)
using the "%e" specifier, and awk (which r.univar uses to perform the
actual calculation) can read that form, r.stats doesn't have an option
to use it.

Maybe someone will consider extending r.stats to allow the use of "%e"
(or "%g", which uses ordinary decimal notation for values between
0.0001 and 1000000, and exponential notation otherwise).

But, until then, r.univar will effectively round all of your data to
zero. So, you need to scale it upwards so that the significant digits
appear further to the left.

> 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.

absolutely! rounded is closer to the actual value than the trunctated one,
isn't it?

Is there a way to learn about the value range etc. in case of such a low
value raster in Grass 5.03?

I don't know. The range is actually stored as binary data, so the
information exists. However, r.info displays it using "%f", so small
values will be rounded to zero. I don't know of any other program
which will provide this information.

Hamish suggested I should try compiling the new version of r.univar in
src/raster/r.univar2 - would that do? (But then I have to go for the Grass
5.3, right?)

Yes. r.univar2 is a C program, and so doesn't suffer from the problem
which the r.univar script has (i.e. that the values are passed from
r.stats to awk as text, using decimal notation).

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