[GRASS-user] mapcalc and cell coordinates

Glynn Clements glynn at gclements.plus.com
Tue Jul 15 04:52:40 EDT 2008


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


More information about the grass-user mailing list