r.mapcalc question

[Lars Schylberg]

> What I am doing is looking at multiple overlapping raster layers,
> and I want a resultant layer to have the modal value of the inputs,
> as long as they are *not* zero.

I wrote a mode-function for r.mapcalc a couple of years ago,
that might be of use for you. It can handle both multiple layers and
neigbourhood specifications.  It will be included in Grass4.2, but
somehow it was forgotten in grass4.1.  

Please contact me if you are interested and I will provide the code.

Another solution is to construct a mode operation directly in mapcalc.
That is a little more advanced and not as straight forward as my new
function, but I have included a little bourne shell script example
below that you could adopt for your needs.  This is example is with
neighbourhood but you could modify it for multilple layers as well.


Lars

Lars Schylberg                          Email: larss at fmi.kth.se  
Dept. of Geodesy and Photogrammetry 
Royal Institute of Technology (KTH)     Tel.   +46 8 790 86 33   
S-100 44  STOCKHOLM, SWEDEN             Fax.   +46 8 790 66 10

#-----------  Example of a mode filter in r.mapcalc ------------------
#
# Create rule for mapcalc
#
r.stats -cqz $IN | awk '{print $1}' > categories
#
echo $OUT '= \
eval( \' > rulefile
for i in  `cat categories`
do
echo '\
q1 = eval( if( '$IN'[-1,-1] == '$i' ) ) ,\
q2 = eval( if( '$IN'[-1,0] == '$i' ) ) ,\
q3 = eval( if( '$IN'[-1,1] == '$i' ) ) ,\
q4 = eval( if( '$IN'[0,-1] == '$i' ) ) ,\
q5 = eval( if( '$IN'[0,0] == '$i' ) ) ,\
q6 = eval( if( '$IN'[0,1] == '$i' ) ) ,\
q7 = eval( if( '$IN'[1,-1] == '$i' ) ) ,\
q8 = eval( if( '$IN'[1,0] == '$i' ) ) ,\
q9 = eval( if( '$IN'[1,1] == '$i' ) ) ,\
qtot'$i' =  q1 + q2 + q3 + q4 + q5 + q6 + q7 + q8 + q9,\' >> rulefile
done
echo 'qmax = max( \' >> rulefile
for i in  `cat categories`
do
echo 'qtot'$i',\' >> rulefile
done
echo '0), \'>> rulefile
echo 'modeclass = max( \'>> rulefile
for i in  `cat categories`
do
echo 'eval( if( qmax == qtot'$i', '$i' )) ,\'>> rulefile
done
echo '0), \'>> rulefile
echo 'modeclass) \
'>> rulefile
#
# Run mapcalc
#
r.mapcalc < rulefile
#