[GRASS-user] r.mapcalc: columns()/modulo - rounding issue, possible bugfeature ?

[Hamish]

peter.loewe wrote:
> I stumbled over this "funny" behavior of the col() function of
> r.mapcalc (GRASS6.2 and GRASS6.4, both on Linux). 
> 
> According to The Book, col() "returns current col of moving window".
> 
> I assume that for each row of a given raster layer col()
> should provide the values 1,2,3,4,5, etc. as it travels
> through the row. 
> 
> To produce a map of repeating value sequences (for example
> the numbers 1 2 3 4 ) the module operator can be used on the
> col()-function: 
> 
> col()  operator result
> 1      mod 4     = 1
> 2      mod 4     = 2
> 3      mod 4     = 3
> 4      mod 4     = 0
> 5      mod 4     = 1
> 6      mod 4     = 2
> etc
> 
> this approach was put into mapcalc. The "1 + "
> makes sure that the sequence  runs from 1-4 not 0 - 3:
> 
> r.mapcalc modulo="(1 + ( col() % 4 ) )"
> 
> A map was produced in a Spearfish location (col() should
> operate independently from the projection).
> 
> Here it comes: 
> For some (strange ?)reason, there are "ripples"
> in the "1 - 4" sequence.
> 
> Instead of going 1234 1234 1234 ad infinitum once in a
> while a number is omitted: 1234 1234 234  1234 134 etc.
> 
> Has anybody come across this behavior before ?
> Is this a feature/good-thing ?


For the the xmon looks ok to me,

#spearfish + grass 6.4rc
g.region -d
  #rows:       477
  #cols:       634
r.mapcalc modulo="(1 + ( col() % 4 ) )"

r.stats modulo -c  # count number of cells
 100%
1 75366
2 75843
3 75843
4 75366

# 75843 - 75366 = 477 = number of rows.
# zooming way in and looking with d.what.rast I see the left-most
#  column (western edge) is 2s and the right-most column (eastern edge)
#  is 3s.

r.univar modulo
sum: 756045

#rows:       477
#cols:       634

634 % 4 = 2
634 / 4 = 158  (truncated to int)

( 1*158 + 2*159 + 3*159 + 4*158 ) * 477  =  756045
which is the same as r.univar reports.


so only thing I see is that the "col()%4 +1" doesn't need the +1 as
col() starts at 1 not 0. For the missing value mid-stream I would guess
slight sampling/aliasing error?? try running "g.region rast=modulo" first
and see if you get the same result.


Hamish