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

GRASS GIS GRASS GIS developer
1

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?

Rgds

Paulo

2

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@gclements.plus.com>

3

Hi Glynn

Thanks for the explanation. Good to know that the problem lies with the exceeding the integer range… something I have to be careful about in the future.

···

On Thu, Sep 25, 2014 at 1:04 AM, Glynn Clements <glynn@gclements.plus.com> wrote:

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@gclements.plus.com>