[GRASS-CVS] markus: grass6/scripts/r.regression.line
description.html, 1.2, 1.2.6.1 r.regression.line, 1.15, 1.15.2.1
grass at intevation.de
grass at intevation.de
Tue Dec 4 15:41:43 EST 2007
Author: markus
Update of /grassrepository/grass6/scripts/r.regression.line
In directory doto:/tmp/cvs-serv8110
Modified Files:
Tag: releasebranch_6_3
description.html r.regression.line
Log Message:
sync'ed to HEAD
Index: description.html
===================================================================
RCS file: /grassrepository/grass6/scripts/r.regression.line/description.html,v
retrieving revision 1.2
retrieving revision 1.2.6.1
diff -u -d -r1.2 -r1.2.6.1
--- description.html 30 Aug 2005 14:57:07 -0000 1.2
+++ description.html 4 Dec 2007 20:41:41 -0000 1.2.6.1
@@ -2,20 +2,43 @@
<EM>r.regression.line</EM> Calculates linear regression from two raster maps,
according to the formula y = a + b*x, where x and y represent raster maps.
-Optionally saves regression coefficients to an ASCII file.
+Optionally saves regression coefficients to an ASCII file.
+The result includes the following coefficients:
+offset (a) and gain (b), residuals (R),
+number of elements (N), medians (medX, medY), standard deviations
+(sdX, sdY), and the F test for testing the significance of the
+regression model as a whole (F).
<br>
<H2>EXAMPLE</H2>
Comparison of the old and the new DEM in Spearfish:
-<PRE>
-g.region rast=elevation.10m
-r.regression.line elevation.dem map2=elevation.10m
-</PRE>
+<div class="code"><pre>
+g.region rast=elevation.10m -p
+r.regression.line map1=elevation.dem map2=elevation.10m
+</pre></div>
+<p>
+
+Using the script style flag AND <em>eval</em> to make results
+available in the shell:
+<div class="code"><pre>
+g.region rast=elevation.10m -p
+eval `r.regression.line -g map1=elevation.dem map2=elevation.10m`
+echo $a
+479.615
+
+echo $b
+0.645631
+
+echo $R
+0.804441
+</pre></div>
+
<H2>AUTHOR</H2>
-Dr. Agustin Lobo - alobo at ija.csic.es<BR>
-Updated to GRASS 5.7 Michael Barton, Arizona State University
+Dr. Agustin Lobo - alobo at ija.csic.es<BR>
+Updated to GRASS 5.7 Michael Barton, Arizona State University<BR>
+Script style output Markus Neteler
<p><i>Last changed: $Date$</i>
Index: r.regression.line
===================================================================
RCS file: /grassrepository/grass6/scripts/r.regression.line/r.regression.line,v
retrieving revision 1.15
retrieving revision 1.15.2.1
diff -u -d -r1.15 -r1.15.2.1
--- r.regression.line 29 Mar 2007 09:34:39 -0000 1.15
+++ r.regression.line 4 Dec 2007 20:41:41 -0000 1.15.2.1
@@ -19,6 +19,10 @@
#% description: Calculates linear regression from two raster maps: y = a + b*x
#% keywords: raster, statistics
#%End
+#%flag
+#% key: g
+#% description: Print in shell script style
+#%End
#%option
#% key: map1
#% type: string
@@ -63,6 +67,21 @@
LC_NUMERIC=C
export LC_NUMERIC
+ok=yes
+eval `g.findfile element=cell file=$GIS_OPT_MAP1`
+if [ -z "$name" ] ; then
+ g.message -e "Raster map <$GIS_OPT_MAP1> not found"
+ ok=no
+fi
+eval `g.findfile element=cell file=$GIS_OPT_MAP2`
+if [ -z "$name" ] ; then
+ g.message -e "Raster map <$GIS_OPT_MAP2> not found"
+ ok=no
+fi
+if [ "$ok" = "no" ] ; then
+ exit 1
+fi
+
TMP="`g.tempfile pid=$$`"
if [ $? -ne 0 ] || [ -z "$TMP" ] ; then
g.message -e "Unable to create temporary files"
@@ -71,7 +90,7 @@
#calculate regression equation
-r.stats -c input=$GIS_OPT_MAP1,$GIS_OPT_MAP2 > "$TMP"
+r.stats -cnA input=$GIS_OPT_MAP1,$GIS_OPT_MAP2 > "$TMP"
awk '{tot += $3;sumX +=$1 * $3; sumsqX +=$1*$1*$3;sumY +=$2 * $3; sumsqY +=$2*$2*$3;\
sumXY +=$1*$2*$3;\
}\
@@ -81,32 +100,42 @@
mediaY=sumY/tot;sumsqY=sumsqY/tot;varY=sumsqY-(mediaY^2);sdY=varY^0.5;\
A=mediaY - B*mediaX; F= R^2/(1-R^2/tot-2);\
print A, B, R, tot, F, mediaX, sdX, mediaY, sdY}' "$TMP" > "$TMP"b
-#cat "$TMP"b
+
+echo "a b R N F medX sdX medY sdY" | tr -s ' ' '\n' > "$TMP"d
+cat "$TMP"b | tr -s ' ' '\n' > "$TMP"e
#send output to screen or text file
if [ -z "$GIS_OPT_OUTPUT" ]
then
- RESULTADO=`cat "$TMP"b`
- echo "y = a + b*x"
- echo " a: offset"
- echo " b: gain"
- echo " R: sumXY - sumX*sumY/tot"
- echo " N: number of elements"
- echo " medX, medY: Medians"
- echo " sdX, sdY: Standard deviations"
- echo "a b R N F medX sdX medY sdY"
- echo ${RESULTADO}
+ if [ $GIS_FLAG_G -eq 1 ] ; then
+ paste -d'=' "$TMP"d "$TMP"e
+ else
+ echo "y = a + b*x"
+ echo " a: offset"
+ echo " b: gain"
+ echo " R: sumXY - sumX*sumY/tot"
+ echo " N: number of elements"
+ echo " medX, medY: Medians"
+ echo " sdX, sdY: Standard deviations"
+ echo "a b R N F medX sdX medY sdY"
+ RESULTADO=`cat "$TMP"b`
+ echo ${RESULTADO}
+ fi
else
- echo "y = a + b*x"
- echo " a: offset"
- echo " b: gain"
- echo " R: sumXY - sumX*sumY/tot"
- echo " N: number of elements"
- echo " medX, medY: Medians"
- echo " sdX, sdY: Standard deviations"
- echo "a b R N F medX sdX medY sdY" > "$TMP"c
- cat "$TMP"b >> "$TMP"c
+ if [ $GIS_FLAG_G -eq 1 ] ; then
+ paste -d'=' "$TMP"d "$TMP"e >> "$TMP"c
+ else
+ echo "y = a + b*x" > "$TMP"c
+ echo " a: offset" >> "$TMP"c
+ echo " b: gain" >> "$TMP"c
+ echo " R: sumXY - sumX*sumY/tot" >> "$TMP"c
+ echo " N: number of elements" >> "$TMP"c
+ echo " medX, medY: Medians" >> "$TMP"c
+ echo " sdX, sdY: Standard deviations" >> "$TMP"c
+ echo "a b R N F medX sdX medY sdY" >> "$TMP"c
+ cat "$TMP"b >> "$TMP"c
+ fi
mv "$TMP"c "$GIS_OPT_OUTPUT"
cat "$GIS_OPT_OUTPUT"
fi
-rm -f "$TMP" "$TMP"b "$TMP"c
+rm -f "$TMP" "$TMP"b "$TMP"c "$TMP"d "$TMP"e
More information about the grass-commit
mailing list