stats/std dev

[Michael Shapiro]

|I wanted to know the standard deviation of the cell values of a raster map
|and didn't see a way to get it from the current tools, so I made a 
|shell script to do it using r.stats and awk. I didn't follow the rules
|by writing it in 'sh' (it is 'csh'), but it is easily converted.
|(the whole thing boils down to a one line command, really)
|Anyway, I thought I would share it to the mailing list to mix in with
|all the subscribe, unsubscribe, and test messages we've all been getting
|biffed with ;-}.
|
|--Chris Rewerts, Agricultural Engineering, Purdue University   
|
Here is a script in /bin/sh that does 
	min, max, mean, mode and standard deviation
-------------------------------------------------------------
:
while test $# != 0
do
	case "$1" in
	-z) z=z;shift;;
	-v) v=v;shift;;
	-zv|-vz) z=z;v=v;shift;;
	-|-*) oops=yes;break;;
	*) break;;
	esac
done
if test $# != 1 -o "$oops" = yes
then
	echo "Usage: `basename $0` [-vz] cellfile" >&2
	exit 1
fi
r.stats -c$z$v "$1" | awk '
   BEGIN{sum=0.0;sum2=0.0}
   NR==1{min=$1; max=$1}
	{sum += $1 * $2; sum2 += $1 * $1 * $2; N += $2}
	{if($1 > max) max = $1; if ($1 < min) min = $1}
	{if($2 > modecount) {mode=$1; modecount=$2}}
     END{
         print "min      ", min
	 print "max      ", max
	 print "mean     ", sum/N
	 print "mode     ", mode
	 print "variance ", (sum2 - sum*sum/N)/N
	 print "deviation", sqrt((sum2 - sum*sum/N)/N)
        }'