[GRASS-user] mapcalc and cell coordinates

G. Allegri giohappy at gmail.com
Mon Jul 14 09:18:48 EDT 2008


Thanks Glynn. Exaustive explanation :-)
So I just need to check if the difference abs($x-x()) is lower then
half the resolution...

2008/7/14 Glynn Clements <glynn at gclements.plus.com>:
>
> 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>
>


More information about the grass-user mailing list