[GRASS-dev] Floating and integer numbers in r.mapcalc

Glynn Clements glynn at gclements.plus.com
Wed Sep 24 16:04:05 PDT 2014


Paulo van Breugel wrote:

> I have the following function (to calculate cell area)
> 
> PI=3.1415926536
> RES=0.0008333
> R=6371007
> r.mapcalc "test3 = (sin(y()+$RES/2) - sin(y()-$RES/2)) * ($RES * $PI/180) *
> $R^2"
> 
> This will give me a wrong results, apparently because R is an integer
> value. If I adapt the formula:
> 
> r.mapcalc "test3 = (sin(y()+$RES/2) - sin(y()-$RES/2)) * ($RES * $PI/180) *
> float($R)^2"
> 
> I get the correct result. From the help file I understand that if any of
> the arguments is a floating number the result should be a floating number.
> But that does not seem the case above so apparently I am missing something
> here. Can anybody explain the logic behind the above?

In the first case, 6371007^2 is calculated using integer arithmetic,
and the result (40589730194049) will exceed the range of an "int"
(2^31-1 = 2147483647), typically resulting in overflow (R^2 will
probably equal -2005720447). The value will be converted to "double"
for the multiplication, but by that point the error has already
occurred.

In the second case, R is converted to "float" (single-precision
floating-point), and the result of squaring it is well within the
range of a "float", although it will be rounded to 40589731037184.0
(using double() instead of float() will use double-precision and avoid
the rounding error).

-- 
Glynn Clements <glynn at gclements.plus.com>


More information about the grass-dev mailing list