[GRASS-SVN] r42525 -
grass/branches/develbranch_6/scripts/d.correlate
svn_grass at osgeo.org
svn_grass at osgeo.org
Tue Jun 8 19:02:19 EDT 2010
Author: hamish
Date: 2010-06-08 19:02:18 -0400 (Tue, 08 Jun 2010)
New Revision: 42525
Modified:
grass/branches/develbranch_6/scripts/d.correlate/d.correlate
Log:
- add a flag to draw the trend line (still buggy!);
- better cmd line parsing;
- update from grass5 bg color;
- tune text rendering;
- use more than 255 dots, but not more than 2048
Modified: grass/branches/develbranch_6/scripts/d.correlate/d.correlate
===================================================================
--- grass/branches/develbranch_6/scripts/d.correlate/d.correlate 2010-06-08 21:06:13 UTC (rev 42524)
+++ grass/branches/develbranch_6/scripts/d.correlate/d.correlate 2010-06-08 23:02:18 UTC (rev 42525)
@@ -47,6 +47,10 @@
#% description: raster input map
#% required : no
#%end
+#%flag
+#% key: t
+#% description: Draw a trend line
+#%end
if [ -z "$GISBASE" ] ; then
echo "You must be in GRASS GIS to run this program." >&2
@@ -104,9 +108,9 @@
fi
# how many cmd line arguments?
-ARGNUM=`echo "$CMDLINE" | tr -s ' ' '\n' | wc -l | awk '{print $1 - 1}'`
+ARGNUM=`echo "$CMDLINE" | tr -s ' ' '\n' | grep -c 'layer[0-9]='`
-echo "CORRELATION" | d.text color=white size=4 line=1
+echo "CORRELATION" | d.text color=grey size=3 line=1
colors="red black blue green gray violet"
if [ $ARGNUM -eq 2 ] ; then
line=4
@@ -117,13 +121,17 @@
jloop=0
# get max in case of two maps for x, y axis
-eval `r.univar -g $GIS_OPT_LAYER1`
-max_layer1=$max
-eval `r.univar -g $GIS_OPT_LAYER2`
-max_layer2=$max
+eval `r.univar -g "$GIS_OPT_LAYER1"`
+max_layer1=`echo "$max" | awk '{printf("%f", $1)}'`
+eval `r.univar -g "$GIS_OPT_LAYER2"`
+max_layer2=`echo "$max" | awk '{printf("%f", $1)}'`
-for i in $GIS_OPT_LAYER1 $GIS_OPT_LAYER2 $GIS_OPT_LAYER3 $GIS_OPT_LAYER4
-do
+# don't go too overboard with r.stats
+if [ "$n" -gt 2048 ] ; then
+ n=2048
+fi
+
+for i in $GIS_OPT_LAYER1 $GIS_OPT_LAYER2 $GIS_OPT_LAYER3 $GIS_OPT_LAYER4 ; do
iloop=`expr $iloop + 1`
for j in $GIS_OPT_LAYER1 $GIS_OPT_LAYER2 $GIS_OPT_LAYER3 $GIS_OPT_LAYER4 ; do
jloop=`expr $jloop + 1`
@@ -133,27 +141,34 @@
colors=`echo $colorstmp2 $colorstmp1`
if [ $ARGNUM -eq 2 ] ; then
- echo $j | d.text color=`echo $colors | cut -d' ' -f1` size=4 at=0,9$line
+ echo $j | d.text color=`echo $colors | cut -d' ' -f1` size=4 at=1,9$line
echo $i | d.text color=`echo $colors | cut -d' ' -f1` size=4 at=60,0$line
else
echo $i $j | d.text color=`echo $colors | cut -d' ' -f1` size=4 line=$line
fi
+
line=`expr $line + 1`
- r.stats -cnA input=$i,$j > "$TMP1"
- m="`awk '$1>max1{max1=$1} $2>max2{max2=$2} min1==0||$1<min1{min1=$1} min2==0||$2<min2{min2=$2} END {print min1,max1,min2,max2}' $TMP1`"
- m1=`echo $m | cut -d' ' -f1`
- m2=`echo $m | cut -d' ' -f2`
- m3=`echo $m | cut -d' ' -f3`
- m4=`echo $m | cut -d' ' -f4`
- awk '{print "move",($1-min1+1)*100.0/(max1-min1+1),($2-min2+1)*100.0/(max2-min2+1);print "draw",($1-min1+1)*100.0/(max1-min1+1),($2-min2+1)*100.0/(max2-min2+1) }' min1=$m1 max1=$m2 min2=$m3 max2=$m4 "$TMP1" | d.graph color=`echo $colors | cut -d' ' -f1`
+ r.stats -cnA input=$i,$j nsteps=$n > "$TMP1"
+
+ m=`awk '$1>max1{max1=$1} $2>max2{max2=$2} min1==0||$1<min1{min1=$1} min2==0||$2<min2{min2=$2} END {print min1,max1,min2,max2}' "$TMP1"`
+ g.message -d message="min1,max1,min2,max2: $m"
+ m1=`echo "$m" | cut -d' ' -f1`
+ m2=`echo "$m" | cut -d' ' -f2`
+ m3=`echo "$m" | cut -d' ' -f3`
+ m4=`echo "$m" | cut -d' ' -f4`
+ awk '{print "move", ($1-min1+1)*100.0/(max1-min1+1), ($2-min2+1)*100.0/(max2-min2+1);
+ print "draw", ($1-min1+1)*100.0/(max1-min1+1), ($2-min2+1)*100.0/(max2-min2+1) }' \
+ min1=$m1 max1=$m2 min2=$m3 max2=$m4 "$TMP1" | \
+ d.graph color=`echo $colors | cut -d' ' -f1`
+
if [ $ARGNUM -eq 2 ] ; then
d.graph << EOF
-size 2 2
-move 0 92
-text $max_layer1
-move 90 2
-text $max_layer2
+ size 2 2
+ move 1 90
+ text max: $max_layer1
+ move 75 2
+ text max: $max_layer2
EOF
fi
fi
@@ -162,3 +177,74 @@
done
rm -f "$TMP1"
+
+#### work in progress / feel free to fix it !!! ####
+
+## FIXME: dots drawn on a normalized scale, line is drawing using absolute scale
+## use above awk dot drawing to output x,y scale factors?
+
+# overrdraw a trend line, slope+offset, and R^2 value.
+if [ "$GIS_FLAG_T" -eq 1 ] ; then
+ START_X=`echo "$m1" "$m2" | awk '{print 100./($2 - $1 +1)}'`
+ START_Y=`echo "$m3" "$m4" | awk '{print 100./($2 - $1 +1)}'`
+ # calc this time just to get scaled slope
+ END_X=`echo "$m1" "$m2" | awk '{print ($2-$1+1)*100./($2 - $1 +1)}'`
+ END_Y=`echo "$m3" "$m4" | awk '{print ($2-$1+1)*100./($2 - $1 +1)}'`
+ #echo 1. $START_X,$START_Y $END_X,$END_Y
+
+# an attempt to scale the y=a+bx slope value into scaled-to-xmon space: ???
+ TRUE_SLOPE=`echo "$m1 $m2 $m3 $m4" | awk '{print ($4-$3)*1.0/($2-$1)}'`
+ SCALED_SLOPE=`echo "$START_X $START_Y $END_X $END_Y" | awk '{print ($4-$2)*1.0/($3-$1)}'`
+
+ g.message -d "true slope: $TRUE_SLOPE"
+ g.message -d "scaled slope: $SCALED_SLOPE"
+
+ eval `r.regression.line -sg map1="$GIS_OPT_LAYER2" map2="$GIS_OPT_LAYER1"`
+ # redo it, just to confirm
+ END_X=100
+ END_L=`echo "$a $b" | awk '{print ( $1 + (100 * $2) )}'`
+# END_Y=`echo "$m3 $m4 $END_Ym" | awk '{print ($3-$1+1)*100./($2 - $1 +1)}'`
+
+ g.message -d message="endL=$END_L m3=$m3 m4=$m4 st_y=$START_Y b=$b"
+ #echo "$START_Y $m3 $m4 $SCALED_SLOPE $b"
+ #no END_Y=`echo "$START_Y $m2 $SCALED_SLOPE" | awk '{print ($1 + ($3 * $2))}'`
+
+### not correct...
+ g.message -i "FIXME: red slope line rendering is not yet correct."
+ #END_Y=`echo "$START_Y $m2 $b $TRUE_SLOPE $SCALED_SLOPE" | awk '{print $1 + ($2 * $3*$4*$5)}'`
+ END_Y=`echo "$END_L $m3 $m4" | awk '{print $1/($3-$2+1)}'`
+
+
+ #echo 2. $START_X,$START_Y $END_X,$END_Y $END_L
+
+#if [ 0 -eq 1 ] ; then
+ # d.graph doesn't like negative coords or coords > 100%
+ if [ `echo "$a" | grep -c '\-'` -ne 0 ] ; then
+ START_Y=0
+ START_X=`echo "$a $b" | awk '{print (-1 * $1 / $2)}'`
+ fi
+ if [ `echo "$END_Y" | awk '{if($1 > 100) {print 1} else {print 0}}'` -eq 1 ] ; then
+ END_Y=100
+ END_X=`echo "$a $b" | awk '{print ((100 - $1) / $2)}'`
+ fi
+ #TODO: ensure other X,Ys do not exceed [0,100]
+#fi
+
+ #echo 3. $START_X,$START_Y $END_X,$END_Y
+ d.graph << EOF
+ color red
+ width 2
+ move $START_X $START_Y
+ draw $END_X $END_Y
+ width 0
+
+ color indigo
+ size 2
+ move 5 80
+ text y = $b*x + $a
+ move 5 75
+ text R^2 = `echo "$R" | awk '{printf("%.4g", $1 * $1)}'`
+# um, r.regression.line's R is the R of R^2 fame, right??
+EOF
+fi
+
More information about the grass-commit
mailing list