How can I control mapalgebra "if" statements on the base of the
coordinates values of a cell?
The base problem is to change the values of some specific cells on the
base of their coordinates.
I've tried using d.rast.edit but it seems a little "contorted" to me,
having to specify so many parameters to change just a value...
Giovanni
On Mon, 14 Jul 2008, G. Allegri wrote:
How can I control mapalgebra "if" statements on the base of the
coordinates values of a cell?
The base problem is to change the values of some specific cells on the
base of their coordinates.
I suspect this might not be possible with r.mapcalc, because each variable in an r.mapcalc expression is effectively an entire map and the expression used to generate the output map can't be changed based on the co-ordinates of individual cells in the input map. I.e., the same expression applies to each and every cell in the output map.
That's my understanding anyway, based on an answer Glynn gave me a few years ago:
http://lists.osgeo.org/pipermail/grass-user/2003-January/008182.html
Does that seem relevant to your case?
Paul
Yer, but as mapcalc has some built-in functions like x() and y() to
retrieve the x-y coordinates, I think it could be easy to trasform
them in projected coordinates.
I'm writing a few lines in C to solve it, but I think I'll take a look
to r.mapcalc code to see if it would be possible to add something like
x_coord() and y_coord()...
2008/7/14 Paul Kelly <paul-grass@stjohnspoint.co.uk>:
On Mon, 14 Jul 2008, G. Allegri wrote:
How can I control mapalgebra "if" statements on the base of the
coordinates values of a cell?
The base problem is to change the values of some specific cells on the
base of their coordinates.I suspect this might not be possible with r.mapcalc, because each variable
in an r.mapcalc expression is effectively an entire map and the expression
used to generate the output map can't be changed based on the co-ordinates
of individual cells in the input map. I.e., the same expression applies to
each and every cell in the output map.That's my understanding anyway, based on an answer Glynn gave me a few years
ago:
http://lists.osgeo.org/pipermail/grass-user/2003-January/008182.html
Does that seem relevant to your case?Paul
G. Allegri wrote:
How can I control mapalgebra "if" statements on the base of the
coordinates values of a cell?
The base problem is to change the values of some specific cells on the
base of their coordinates.
r.mapcalc "$outmap = if(abs(x() - $x) < $eps && abs(y() - $y) < $eps,$newval,$inmap)"
where everything preceded with a $ sign is a shell variable (or you
can just insert literal values into the expression).
Needless to say, this gets rather ugly (and inefficient) if you want
to change more than a few cells. In that situation, you would probably
be better off using r.in.xyz and r.patch.
--
Glynn Clements <glynn@gclements.plus.com>
On 14/07/08 12:17, G. Allegri wrote:
Yer, but as mapcalc has some built-in functions like x() and y() to
retrieve the x-y coordinates, I think it could be easy to trasform
them in projected coordinates.
I'm writing a few lines in C to solve it, but I think I'll take a look
to r.mapcalc code to see if it would be possible to add something like
x_coord() and y_coord()...2008/7/14 Paul Kelly <paul-grass@stjohnspoint.co.uk>:
On Mon, 14 Jul 2008, G. Allegri wrote:
How can I control mapalgebra "if" statements on the base of the
coordinates values of a cell?
The base problem is to change the values of some specific cells on the
base of their coordinates.I suspect this might not be possible with r.mapcalc, because each variable
in an r.mapcalc expression is effectively an entire map and the expression
used to generate the output map can't be changed based on the co-ordinates
of individual cells in the input map. I.e., the same expression applies to
each and every cell in the output map.That's my understanding anyway, based on an answer Glynn gave me a few years
ago:
http://lists.osgeo.org/pipermail/grass-user/2003-January/008182.html
Does that seem relevant to your case?
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...
Moritz
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...
So, it's possible to check above a range of coordinates but not on a
specific one...
But I don't understand a thing:
r.mapcalc/coor.c uses G_northing_to_row() and G_easting_to_row() as
well as r.what, for example. Why the issue about precision doesn't
happen in this case?
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@gclements.plus.com>
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@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@gclements.plus.com>
On 14/07/08 15:14, Glynn Clements wrote:
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.
Just for completeness: It is a cell map, but resolution was aligned on bounds, so something like 1.9876 and 2.0045
When I set resolution to 2 and align bounds to resolution and then use
r.mapcalc brol3=if(y()==9521601 && x()==427001,1,null())
it works
whereas
r.mapcalc brol3=if(y()==9521600 && x()==427000,1,null())
doesn't. So clearly an issue of falling right on the center of the pixel.
Also note that the calculation involves multiplying by the resolution.
You mean internally by mapcalc, or that I should include resolution into the calculation ?
If the resolution cannot be exactly represented as both a binary
fraction and a decimal fraction, you will get rounding error.
That and/or hitting the cell center was probably the problem.
Thanks for your explanations.
Moritz
Moritz Lennert wrote:
> Also note that the calculation involves multiplying by the resolution.
You mean internally by mapcalc, or that I should include resolution into
the calculation ?
Internally.
The point is that if you have e.g. a 100m x 100m square, and you set
rows = cols = 300, so each cell is 1/3rd of a metre, *all* of the cell
centres will suffer from rounding error, even when the row and column
are multiples of 6, as the calculation will multiply by (a binary
approximation to) 1/3, not divide by 3.
--
Glynn Clements <glynn@gclements.plus.com>
G. Allegri wrote:
Thanks Glynn. Exaustive explanation
So I just need to check if the difference abs($x-x()) is lower then
half the resolution...
In hindsight:
If you want to test for the specific cell in which the given
coordinates lie, you may be better off converting the coordinates to a
row/col based upon the region settings, then testing against the
results of row() and col() (but note: row() and col() are 1-based, not
0-based).
Testing against x() and y() has the potential for a value which falls
exactly on (or very close to) the boundary between cells to match
either 2 (or 4) cells, or none, due to rounding error. Particularly if
the resolution isn't exactly representable in both binary and decimal.
--
Glynn Clements <glynn@gclements.plus.com>
If you want to test for the specific cell in which the given
coordinates lie, you may be better off converting the coordinates to a
row/col based upon the region settings, then testing against the
results of row() and col() (but note: row() and col() are 1-based, not
0-based).Testing against x() and y() has the potential for a value which falls
exactly on (or very close to) the boundary between cells to match
either 2 (or 4) cells, or none, due to rounding error. Particularly if
the resolution isn't exactly representable in both binary and decimal.
You're rigth Glynn, but how to check the row and col number in bash?I
could calculate it with the following straigthforward way
$x/(east-west) = $col/cols --> $col = $x*cols/(east-west)
and then use "floor($col)" as the col value.
Yet, could the precision problem arise in this case too?
G. Allegri wrote:
> If you want to test for the specific cell in which the given
> coordinates lie, you may be better off converting the coordinates to a
> row/col based upon the region settings, then testing against the
> results of row() and col() (but note: row() and col() are 1-based, not
> 0-based).
>
> Testing against x() and y() has the potential for a value which falls
> exactly on (or very close to) the boundary between cells to match
> either 2 (or 4) cells, or none, due to rounding error. Particularly if
> the resolution isn't exactly representable in both binary and decimal.You're rigth Glynn, but how to check the row and col number in bash?I
could calculate it with the following straigthforward way$x/(east-west) = $col/cols --> $col = $x*cols/(east-west)
and then use "floor($col)" as the col value.
bash doesn't do floating-point, so you need to use bc or awk, e.g.:
eval `g.region -g`
col=`echo "1 + ($x - $w) / $ewres" | bc`
row=`echo "1 + ($y - $s) / $nsres" | bc`
r.mapcalc "$newmap = if(row() == $row && col() == $col,$newval,$oldmap)"
Yet, could the precision problem arise in this case too?
You could get some kind of precision problem whichever approach you
use.
However, comparing the results of row() and col() (which return an
integer) means that you will always set exactly one cell (assuming
that the coordinates are within the bounds of the map). If the
original coordinates are exactly on the border between two cells, you
might set the "wrong" one, but it will always set one cell, never both
or neither.
--
Glynn Clements <glynn@gclements.plus.com>
Ok, Glynn. I'll adopt this way.
bash doesn't do floating-point, so you need to use bc or awk, e.g.:
Yes, ok, I know it. It was just a "pseduo-code"