[GRASS-SVN] r57776 - grass/branches/develbranch_6/scripts/d.correlate
svn_grass at osgeo.org
svn_grass at osgeo.org
Sat Sep 21 03:51:56 PDT 2013
Author: hamish
Date: 2013-09-21 03:51:56 -0700 (Sat, 21 Sep 2013)
New Revision: 57776
Modified:
grass/branches/develbranch_6/scripts/d.correlate/d.correlate
grass/branches/develbranch_6/scripts/d.correlate/description.html
Log:
fix linear regression trendline plotting over the top of the scatterplot; assorted minor whitespace & quoting cleanup
Modified: grass/branches/develbranch_6/scripts/d.correlate/d.correlate
===================================================================
--- grass/branches/develbranch_6/scripts/d.correlate/d.correlate 2013-09-21 03:27:42 UTC (rev 57775)
+++ grass/branches/develbranch_6/scripts/d.correlate/d.correlate 2013-09-21 10:51:56 UTC (rev 57776)
@@ -8,18 +8,17 @@
# H.Bowman added trend line & model fix flag June 2010
# PURPOSE: prints a graph of the correlation between data layers (in pairs)
# derived from <grass5>/src.local/d.correlate.sh
-# COPYRIGHT: (C) 2005 by the GRASS Development Team
+# COPYRIGHT: (C) 2005-2013 by the GRASS Development Team
#
# This program is free software under the GNU General Public
# License (>=v2). Read the file COPYING that comes with GRASS
# for details.
#
-# TODO GRASS 7: rename parameters to map1, map2 or better multiple map=map1[,map2[...]]
#############################################################################
#%Module
-#% description: Prints a graph of the correlation between data layers (in pairs).
-#% keywords: display, diagram
+#% description: Prints a graph of the correlation between raster maps (in pairs).
+#% keywords: display, correlation, scatterplot
#%End
#%option
#% key: layer1
@@ -51,7 +50,7 @@
#%end
#%flag
#% key: t
-#% description: Draw a trend line
+#% description: Draw a trend line and calculate linear regression formula
#%end
if [ -z "$GISBASE" ] ; then
@@ -68,7 +67,6 @@
exec g.parser "$0" "$@"
fi
-PROG=`basename "$0"`
#### check if we have awk
if [ ! -x "`which awk`" ] ; then
@@ -82,162 +80,203 @@
export LC_NUMERIC
if [ $# -gt 4 ] ; then
- g.message -e "max 4 layers allowed"
- exit 1
+ g.message -e "A maximum of four raster maps is allowed"
+ exit 1
fi
-TMP1="`g.tempfile pid=$$`"
+if [ "$GIS_OPT_LAYER1" = "$GIS_OPT_LAYER2" ] ; then
+ g.message -e "Nothing to see if the maps are the same."
+ exit 1
+fi
ok=yes
-for f in $GIS_OPT_LAYER1 $GIS_OPT_LAYER2 $GIS_OPT_LAYER3 $GIS_OPT_LAYER4
-do
- eval `g.findfile element=cell file=$f`
- if [ -z "$name" ] ; then
- g.message -e "$f not found"
+for map in "$GIS_OPT_LAYER1" "$GIS_OPT_LAYER2" $GIS_OPT_LAYER3 $GIS_OPT_LAYER4 ; do
+ eval `g.findfile element=cell file="$map"`
+ if [ -z "$name" ] ; then
+ g.message -e "Raster map <$map> not found"
ok=no
- fi
+ fi
done
if [ "$ok" = "no" ] ; then
- exit 1
+ exit 1
fi
d.erase
+if [ $? -ne 0 ] ; then
+ exit 1
+fi
+
+TMP1="`g.tempfile pid=$$`"
+if [ $? -ne 0 ] || [ -z "$TMP1" ] ; then
+ g.message -e "Unable to create temporary file"
+ exit 1
+fi
+
+# avoid output that overwrites one another
+# (check what kind of monitor is selected first, skip for Xmons?)
GRASS_PNG_READ=TRUE
export GRASS_PNG_READ
-if [ $? != 0 ] ; then
- exit 1
-fi
-# how many cmd line arguments?
+
+# count how many maps given on the cmd line
ARGNUM=`echo "$CMDLINE" | tr -s ' ' '\n' | grep -c 'layer[0-9]='`
echo "CORRELATION" | d.text color=grey size=3 line=1
colors="red black blue green gray violet"
-if [ $ARGNUM -eq 2 ] ; then
+if [ "$ARGNUM" -eq 2 ] ; then
+ topline=93
line=4
else
line=2
fi
-iloop=0
-jloop=0
# get max in case of two maps for x, y axis
eval `r.univar -g "$GIS_OPT_LAYER1"`
-max_layer1=`echo "$max" | awk '{printf("%f", $1)}'`
+max_layer1="$max"
eval `r.univar -g "$GIS_OPT_LAYER2"`
-max_layer2=`echo "$max" | awk '{printf("%f", $1)}'`
+max_layer2="$max"
# 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`
- if [ "$i" != "$j" -a $iloop -le $jloop ] ; then
- colorstmp1=`echo $colors | cut -d' ' -f1`
- colorstmp2=`echo $colors | cut -d' ' -f2-`
- colors=`echo $colorstmp2 $colorstmp1`
+i_loop=0
+for i_map in "$GIS_OPT_LAYER1" "$GIS_OPT_LAYER2" $GIS_OPT_LAYER3 $GIS_OPT_LAYER4 ; do
+ i_loop=`expr $i_loop + 1`
+ j_loop=0
+ for j_map in "$GIS_OPT_LAYER1" "$GIS_OPT_LAYER2" $GIS_OPT_LAYER3 $GIS_OPT_LAYER4 ; do
+ j_loop=`expr $j_loop + 1`
- if [ $ARGNUM -eq 2 ] ; then
- 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
+ if [ "$i_map" != "$j_map" -a "$i_loop" -le "$j_loop" ] ; then
+ g.message -v "$i_map vs. $j_map ..."
+ colorstmp1=`echo "$colors" | cut -d' ' -f1`
+ colorstmp2=`echo "$colors" | cut -d' ' -f2-`
+ colors=`echo "$colorstmp2 $colorstmp1"`
- line=`expr $line + 1`
+ if [ "$ARGNUM" -eq 2 ] ; then
+ echo "$i_map" | \
+ d.text color=`echo "$colors" | cut -d' ' -f1` size=4 at=1,"$topline"
+ echo "$j_map" | \
+ d.text color=`echo "$colors" | cut -d' ' -f1` size=4 at=60,"$line"
+ else
+ echo "$i_map vs. $j_map" | \
+ d.text color=`echo "$colors" | cut -d' ' -f1` size=4 line="$line"
+ fi
- r.stats -cnA input=$i,$j nsteps=$n > "$TMP1"
+ line=`expr $line + 1`
- 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`
+ r.stats -cnA input="$i_map,$j_map" nsteps="$n" > "$TMP1"
- if [ $ARGNUM -eq 2 ] ; then
- d.graph << EOF
- size 2 2
- move 1 90
- text max: $max_layer1
- move 75 2
- text max: $max_layer2
+ 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`
+
+ # scaled map1 value is plotted as x, scaled map2 value is plotted as y
+ 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 1 90
+ text max: $max_layer1
+ move 75 2
+ text max: $max_layer2
EOF
- fi
+ fi
fi
done
- jloop=0
done
-rm -f "$TMP1"
+rm -f "$TMP1" "$TMP1".graph
-#### 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.
+#### 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)}'`
+ if [ "$ARGNUM" -ne 2 ] ; then
+ g.message -e 'Will only draw trend line if two map layers are given'
+ exit 1
+ fi
- g.message -d "true slope: $TRUE_SLOPE"
- g.message -d "scaled slope: $SCALED_SLOPE"
-
+ # perform linear regression
eval `r.regression.line -g 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 -v "y = $b*x + $a"
+ g.message -v "R^2 = `echo "$R" | awk '{printf("%.4g", $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)}'`
+ #### calc coords for trend line in map1,map2 space
+ A0="$m1"
+ B0=`echo "$m1 $a $b" | awk '{print ($1-$2)/$3}'`
+ A1="$m2"
+ B1=`echo "$m2 $a $b" | awk '{print ($1-$2)/$3}'`
+ # scale to 0-100% space
+ X0=`echo "$A0 $m1 $m2" | awk '{print ($1-$2+1) * 100.0 / ($3-$2+1)}'`
+ Y0=`echo "$B0 $m3 $m4" | awk '{print ($1-$2+1) * 100.0 / ($3-$2+1)}'`
+ X1=`echo "$A1 $m1 $m2" | awk '{print ($1-$2+1) * 100.0 / ($3-$2+1)}'`
+ Y1=`echo "$B1 $m3 $m4" | awk '{print ($1-$2+1) * 100.0 / ($3-$2+1)}'`
+ g.message -d "Raw plot line: $X0,$Y0 $X1,$Y1"
- #echo 2. $START_X,$START_Y $END_X,$END_Y $END_L
+ # d.graph doesn't like % coords that are >100 or <0, so we need to work backwards
+ # through the scaling and linear regression formulas to recalc coords in-bounds.
+ if [ `echo "$Y0" | grep -c '\-'` -ne 0 ] ; then
+ g.message -d message="y-offset at x=0% is negative, recalc"
+ Y0=0
+ B0=`echo "$m3" | awk '{print $1 - 1}'`
+ A0=`echo "$B0 $a $b" | awk '{print ($1 * $3) + $2}'`
+ X0=`echo "$A0 $m1 $m2" | awk '{print 100.0 * ($1-$2+1) / ($3-$2+1)}'`
+ elif [ `echo "$Y0" | awk '{if($1 > 100) {print 1} else {print 0}}'` -eq 1 ] ; then
+ g.message -d message="y position at left side of plot is > 100%, recalc"
+ Y0=100
+ B0="$m4"
+ A0=`echo "$B0 $a $b" | awk '{print ($1 * $3) + $2}'`
+ X0=`echo "$A0 $m1 $m2" | awk '{print 100.0 * ($1-$2+1) / ($3-$2+1)}'`
+ fi
+ g.message -d "After Y0 fixes: $X0,$Y0 $X1,$Y1"
-#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)}'`
+ if [ `echo "$Y1" | grep -c '\-'` -ne 0 ] ; then
+ g.message -d message="y-offset at x=100% is negative, recalc"
+ Y1=0
+ B1=`echo "$m3" | awk '{print $1 - 1}'`
+ A1=`echo "$B1 $a $b" | awk '{print ($1 * $3) + $2}'`
+ X1=`echo "$A1 $m1 $m2" | awk '{print 100.0 * ($1-$2+1) / ($3-$2+1)}'`
+ elif [ `echo "$Y1" | awk '{if($1 > 100) {print 1} else {print 0}}'` -eq 1 ] ; then
+ g.message -d message="y position at right side of plot is > 100%, recalc"
+ Y1=100
+ B1="$m4"
+ A1=`echo "$B1 $a $b" | awk '{print ($1 * $3) + $2}'`
+ X1=`echo "$A1 $m1 $m2" | awk '{print 100.0 * ($1-$2+1) / ($3-$2+1)}'`
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
+ g.message -d "final line coords: ($X0, $Y0) to ($X1, $Y1)"
+
+
+ # plot the trend line and coefficient text
d.graph << EOF
color red
width 2
- move $START_X $START_Y
- draw $END_X $END_Y
+ move $X0 $Y0
+ draw $X1 $Y1
width 0
color indigo
Modified: grass/branches/develbranch_6/scripts/d.correlate/description.html
===================================================================
--- grass/branches/develbranch_6/scripts/d.correlate/description.html 2013-09-21 03:27:42 UTC (rev 57775)
+++ grass/branches/develbranch_6/scripts/d.correlate/description.html 2013-09-21 10:51:56 UTC (rev 57776)
@@ -2,59 +2,42 @@
<em>d.correlate</em> is a shell (sh(1)) script that
graphically displays the results of an
-
-<em><a href="r.stats.html">r.stats</a></em>
-
+ <em><a href="r.stats.html">r.stats</a></em>
run on two raster map layers. This shell script is useful
for highlighting the correlation (or lack of it) among data
layers (scattergram).
-
<p>
-
The results are displayed in the active display frame on
the user's graphics monitor. <em>d.correlate</em> erases
the active frame before displaying results.
+<p>
+The <b>layer</b><i>n</i> parameters supply the names of two to
+four existing raster map layers to be included in the correlation.
-<h2>OPTIONS</h2>
-<h3>Parameters:</h3>
-
-<dl>
-
-<dt><em>layer1 layer2 </em>[<em>layer3</em> [<em>layer4</em>]]
-
-<dd>The names of from two to four existing raster map layers
-to be included in the correlation.
-</dl>
-
<h2>NOTES</h2>
-This is a shell script that uses
-<em><a href="r.stats.html">r.stats</a></em>
-and the UNIX <em>awk</em> command
-to calculate the correlation among data layers,
-and uses
-<em><a href="d.text.html">d.text</a></em> and
+This is a shell script that uses <em><a href="r.stats.html">r.stats</a></em>
+and the UNIX <em>awk</em> command to calculate the correlation among data
+layers, and uses <em><a href="d.text.html">d.text</a></em> and
<em><a href="d.graph.html">d.graph</a></em> to display the results.
-
<p>
-
If three or four map layers are specified, the correlation
among each combination of two data layers is displayed.
-<h2>FILES</h2>
-This program is simply a shell script. Users are
-encouraged to make their own shell script programs using
-similar techniques. See <KBD>$GISBASE/scripts/d.correlate</KBD>.
+<h2>EXAMPLE</h2>
-<h2>SEE ALSO</h2>
+Compare LANDSAT-7 bands 2 and 3 in the North Carolina sample dataset,
+and over-plotting a trend line:
+<div class="code"><pre>
+ g.region rast=lsat7_2002_20
+ d.correlate -t layer1=lsat7_2002_20 layer2=lsat7_2002_30 --verbose
+</pre></div>
-The UNIX <em>awk</em> command.
+<h2>SEE ALSO</h2>
-<p>
-
<em>
<a href="d.text.html">d.text</a>,
<a href="d.graph.html">d.graph</a>,
@@ -63,13 +46,15 @@
<a href="r.stats.html">r.stats</a>
</em>
-<h2>AUTHOR</h2>
+<h2>AUTHORS</h2>
+
Michael Shapiro,
<a href="http://www.cecer.army.mil/">U.S.Army Construction Engineering
Research Laboratory</a>
<p>
Rewritten to GRASS 6 (from csh to sh) by Markus Neteler<br>
-Trend line and model fix option added by Hamish Bowman
+Linear regression model and trend line plotting by Hamish Bowman
-<p><i>Last changed: $Date$</i>
+<p>
+<i>Last changed: $Date$</i>
More information about the grass-commit
mailing list