[GRASS-user] mapcalc and cell coordinates

[Glynn Clements]

Moritz Lennert wrote:

> Well, this works:
> 
> r.mapcalc brol3=if(y()<=9521600 && x()<=427000,1,null())
> 
> but this doesn't:
> 
> r.mapcalc brol3=if(y()==9521600 && x()==427000,1,null())
> 
> which raises the issue of precision of the coordinates...

What are the region settings? Bear in mind that x() and y() return the
coordinates of the cell's centre, not a corner.

Also note that the calculation involves multiplying by the resolution. 
If the resolution cannot be exactly represented as both a binary
fraction and a decimal fraction, you will get rounding error.

[OTOH, if the resolution is an integer, you should be safe.]

More generally, floating point values usually need to be compared to
within some tolerance; see my answer to the original post.

Using == for floating-point values is extremely fragile. In some
situations, even numbers which should be exactly equal won't actually
compare equal (e.g. because the x86 has 80-bit floating-point
registers, but the value will be truncated to 64 bits if stored in a
"double").

Essentially, expect rounding error to occur unless you *know* that all
values involved, including all intermediate results, can be
represented exactly. Compiler optimisations can have an effect, e.g. 
-ffast-math (this option isn't enabled by any -O switch).

[BTW, never use Borland's compiler for code which does FP arithmetic. 
It "optimises" x/y to x*(1/y), which often introduces rounding error,
and this cannot be disabled.]

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