[GRASS-dev] r.mapcal rand() strangeness

Maciej Sieczka tutey at o2.pl
Mon Feb 25 16:41:28 EST 2008


Glynn Clements pisze:
> Maciej Sieczka wrote:

>> The following silently forces 2147483647 on anything bigger:
>>
>> $ r.mapcalc 'map=9999999999999999999'
>> $ r.info -rt map
>> min=2147483647
>> max=2147483647
>> datatype=CELL

> That is down to the atoi() function used to parse integers.

Sorry if I got it wrong - is r.mapcalc limited to 32 bit integers? If 
so, does it have to?

> I could change it to use strtol(), which sets errno to ERANGE on
> overflow.

Would that make r.mapcalc accept integers bigger than int32?

>> Trying the same by adding a decimal place separator, results in 
>> something bizzare to me:
>>
>> $ r.mapcalc 'map=99999999999999999999999.0'
>> $ r.info -rt map
>> min=99999999999999991611392.000000
>> max=99999999999999991611392.000000
>> datatype=DCELL
> 
> The presence of a dot causes the value to be parsed as a
> floating-point value using atof().
> 
> A double has a precision of ~16 decimal digits, which matches what you
> see above.

So processing any number with more than 16 decimal diggits in r.mapcalc 
must yield such "strange" values? And how many diggits after the decimal 
separator are safe?

>> At the edges of Float32 range other strange things happen:
>>
>> Shouldn't the below be rather 
>> min,max=340282000000000000000000000000000000000.0?:
>>
>> $ r.mapcalc 'map=3.40282e38'
>> $ r.info -rt map
>> min=340282000000000014192072600942972764160.000000
>> max=340282000000000014192072600942972764160.000000
>> datatype=DCELL
>>
>> The the one below should be -340282000000000000000000000000000000000.0 I 
>> guess:
>>
>> $ r.mapcalc 'map=-3.40282e38'
>> $ r.info -rt map
>> min=-340282000000000014192072600942972764160.000000
>> max=-340282000000000014192072600942972764160.000000
>> datatype=DCEL

> Again, you're running into the limits of the precision of the double
> type.
> 
> Also, bear in mind that the conversion from a binary FP value to a
> string of decimal digits is being performed by r.info; r.mapcalc has
> no control over that part of the process.

Does that mean that actually in the raster the value is something else 
than what r.info reports?

> Essentially, the %f format behaves badly for large numbers. If you use
> %e or %g, it will use exponential notation, with a specific relative
> precision rather than a specific absolute precision.

>> At bigger Float64 values r.mapcalc crashes:
>>
>> $r.mapcalc 'map=3.40282e110'
>>
>> *** stack smashing detected ***: r.mapcalc terminated
>> Aborted (core dumped)

> Right. It's sprintf()ing the value into a 64-byte buffer using %f.
> 
> I'll change it to use %g:

That fixes it. Thanks. I have noticed that r.null has the same problem 
though. More modules could? To reproduce, please run the attached script 
  (having set region big enough) and proceed utnil the last, Float64. 
step. r.null with crash with "*** stack smashing detected ***".

I wrote the script to validate Martin's patch for r.out.gdal [1] in 
extreme conditions but I'm running into issues, thus this thread.

[1]http://trac.osgeo.org/grass/attachment/ticket/67/r-out-gdal-no-data.diff

Maciek
-------------- next part --------------
A non-text attachment was scrubbed...
Name: datatypes_grass_gdal.sh
Type: application/x-shellscript
Size: 1745 bytes
Desc: not available
Url : http://lists.osgeo.org/pipermail/grass-dev/attachments/20080225/681715f3/datatypes_grass_gdal.bin


More information about the grass-dev mailing list