← Back to branch summary

~ubuntu-branches/ubuntu/vivid/grass/vivid-proposed

« back to all changes in this revision

Viewing changes to scripts/d.correlate/d.correlate

  • Committer: Package Import Robot
  • Author(s): Bas Couwenberg
  • Date: 2015-02-20 23:12:08 UTC
  • mfrom: (8.2.6 experimental)
  • Revision ID: package-import@ubuntu.com-20150220231208-1u6qvqm84v430b10
Tags: 7.0.0-1~exp1
* New upstream release.
* Update python-ctypes-ternary.patch to use if/else instead of and/or.
* Drop check4dev patch, rely on upstream check.
* Add build dependency on libpq-dev to grass-dev for libpq-fe.h.
* Drop patches applied upstream, refresh remaining patches.
* Update symlinks for images switched from jpg to png.

Show diffs side-by-side

added added

removed removed

Lines of Context:
scripts/d.correlate/d.correlate
1
 
#!/bin/sh
2
 
#
3
 
############################################################################
4
 
#
5
 
# MODULE:       d.correlate for GRASS 6; based on dcorrelate.sh for GRASS 4,5
6
 
# AUTHOR(S):    CERL - Michael Shapiro;
7
 
#                updated to GRASS 6 by Markus Neteler 5/2005
8
 
#                H.Bowman added trend line & model fix flag June 2010
9
 
# PURPOSE:      prints a graph of the correlation between data layers (in pairs)
10
 
#               derived from <grass5>/src.local/d.correlate.sh
11
 
# COPYRIGHT:    (C) 2005-2013 by the GRASS Development Team
12
 
#
13
 
#               This program is free software under the GNU General Public
14
 
#               License (>=v2). Read the file COPYING that comes with GRASS
15
 
#               for details.
16
 
#
17
 
#############################################################################
18
 
 
19
 
#%Module
20
 
#% description: Prints a graph of the correlation between raster maps (in pairs).
21
 
#% keywords: display, correlation, scatterplot
22
 
#%End
23
 
#%option
24
 
#% key: layer1
25
 
#% type: string
26
 
#% gisprompt: old,cell,raster
27
 
#% description: raster input map
28
 
#% required : yes
29
 
#%end
30
 
#%option
31
 
#% key: layer2
32
 
#% type: string
33
 
#% gisprompt: old,cell,raster
34
 
#% description: raster input map
35
 
#% required : yes
36
 
#%end
37
 
#%option
38
 
#% key: layer3
39
 
#% type: string
40
 
#% gisprompt: old,cell,raster
41
 
#% description: raster input map
42
 
#% required : no
43
 
#%end
44
 
#%option
45
 
#% key: layer4
46
 
#% type: string
47
 
#% gisprompt: old,cell,raster
48
 
#% description: raster input map
49
 
#% required : no
50
 
#%end
51
 
#%flag
52
 
#% key: t
53
 
#% description: Draw a trend line and calculate linear regression formula
54
 
#%end
55
 
 
56
 
if [ -z "$GISBASE" ] ; then
57
 
    echo "You must be in GRASS GIS to run this program." >&2
58
 
    exit 1
59
 
fi
60
 
 
61
 
if [ "$1" != "@ARGS_PARSED@" ] ; then
62
 
    CMDLINE=`basename "$0"`
63
 
    for arg in "$@" ; do
64
 
        CMDLINE="$CMDLINE \"$arg\""
65
 
    done
66
 
    export CMDLINE
67
 
    exec g.parser "$0" "$@"
68
 
fi
69
 
 
70
 
 
71
 
#### check if we have awk
72
 
if [ ! -x "`which awk`" ] ; then
73
 
    g.message -e "awk required, please install awk or gawk first"
74
 
    exit 1
75
 
fi
76
 
 
77
 
# setting environment, so that awk works properly in all languages
78
 
unset LC_ALL
79
 
LC_NUMERIC=C
80
 
export LC_NUMERIC
81
 
 
82
 
if [ $# -gt 4 ] ; then
83
 
    g.message -e "A maximum of four raster maps is allowed"
84
 
    exit 1
85
 
fi
86
 
 
87
 
if [ "$GIS_OPT_LAYER1" = "$GIS_OPT_LAYER2" ] ; then
88
 
   g.message -e "Nothing to see if the maps are the same."
89
 
   exit 1
90
 
fi
91
 
 
92
 
ok=yes
93
 
for map in "$GIS_OPT_LAYER1" "$GIS_OPT_LAYER2" $GIS_OPT_LAYER3 $GIS_OPT_LAYER4 ; do
94
 
    eval `g.findfile element=cell file="$map"`
95
 
    if [ -z "$name" ] ; then
96
 
        g.message -e "Raster map <$map> not found"
97
 
        ok=no
98
 
    fi
99
 
done
100
 
 
101
 
if [ "$ok" = "no" ] ; then
102
 
    exit 1
103
 
fi
104
 
 
105
 
d.erase
106
 
if [ $? -ne 0 ] ; then
107
 
    exit 1
108
 
fi
109
 
 
110
 
TMP1="`g.tempfile pid=$$`"
111
 
if [ $? -ne 0 ] || [ -z "$TMP1" ] ; then
112
 
    g.message -e "Unable to create temporary file"
113
 
    exit 1
114
 
fi
115
 
 
116
 
# avoid output that overwrites one another
117
 
#  (check what kind of monitor is selected first, skip for Xmons?)
118
 
GRASS_PNG_READ=TRUE
119
 
export GRASS_PNG_READ
120
 
 
121
 
 
122
 
# count how many maps given on the cmd line
123
 
ARGNUM=`echo "$CMDLINE" | tr -s ' ' '\n' | grep -c 'layer[0-9]='`
124
 
 
125
 
echo "CORRELATION" | d.text color=grey size=3 line=1
126
 
colors="red black blue green gray violet"
127
 
if [ "$ARGNUM" -eq 2 ] ; then
128
 
  topline=93
129
 
  line=4
130
 
else
131
 
  line=2
132
 
fi
133
 
 
134
 
# get max in case of two maps for x, y axis
135
 
eval `r.univar -g "$GIS_OPT_LAYER1"`
136
 
max_layer1="$max"
137
 
eval `r.univar -g "$GIS_OPT_LAYER2"`
138
 
max_layer2="$max"
139
 
 
140
 
# don't go too overboard with r.stats
141
 
if [ "$n" -gt 2048 ] ; then
142
 
   n=2048
143
 
fi
144
 
 
145
 
i_loop=0
146
 
for i_map in "$GIS_OPT_LAYER1" "$GIS_OPT_LAYER2" $GIS_OPT_LAYER3 $GIS_OPT_LAYER4 ; do
147
 
   i_loop=`expr $i_loop + 1`
148
 
   j_loop=0
149
 
   for j_map in "$GIS_OPT_LAYER1" "$GIS_OPT_LAYER2" $GIS_OPT_LAYER3 $GIS_OPT_LAYER4 ; do
150
 
     j_loop=`expr $j_loop + 1`
151
 
 
152
 
     if [ "$i_map" != "$j_map" -a "$i_loop" -le "$j_loop" ] ; then
153
 
        g.message -v "$i_map vs. $j_map ..."
154
 
        colorstmp1=`echo "$colors" | cut -d' ' -f1`
155
 
        colorstmp2=`echo "$colors" | cut -d' ' -f2-`
156
 
        colors=`echo "$colorstmp2 $colorstmp1"`
157
 
 
158
 
        if [ "$ARGNUM" -eq 2 ] ; then
159
 
           echo "$i_map" | \
160
 
               d.text color=`echo "$colors" | cut -d' ' -f1` size=4 at=1,"$topline"
161
 
           echo "$j_map" | \
162
 
               d.text color=`echo "$colors" | cut -d' ' -f1` size=4 at=60,"$line"
163
 
        else
164
 
           echo "$i_map vs. $j_map" | \
165
 
               d.text color=`echo "$colors" | cut -d' ' -f1` size=4 line="$line"
166
 
        fi
167
 
 
168
 
        line=`expr $line + 1`
169
 
 
170
 
        r.stats -cnA input="$i_map,$j_map" nsteps="$n" > "$TMP1"
171
 
 
172
 
        m=`awk '$1 > max1  {max1=$1} \
173
 
                $2 > max2  {max2=$2} \
174
 
                min1 == 0 || $1 < min1  {min1=$1} \
175
 
                min2 == 0 || $2 < min2  {min2=$2} \
176
 
                END \
177
 
                {print min1,max1,min2,max2}' "$TMP1"`
178
 
 
179
 
        g.message -d message="min1,max1,min2,max2: $m"
180
 
 
181
 
        m1=`echo "$m" | cut -d' ' -f1`
182
 
        m2=`echo "$m" | cut -d' ' -f2`
183
 
        m3=`echo "$m" | cut -d' ' -f3`
184
 
        m4=`echo "$m" | cut -d' ' -f4`
185
 
 
186
 
        # scaled map1 value is plotted as x, scaled map2 value is plotted as y
187
 
        awk '{print "move", \
188
 
                ($1-min1+1) * 100.0 / (max1-min1+1), \
189
 
                ($2-min2+1) * 100.0 / (max2-min2+1);
190
 
              print "draw", \
191
 
                ($1-min1+1) * 100.0 / (max1-min1+1), \
192
 
                ($2-min2+1) * 100.0 / (max2-min2+1) }' \
193
 
              min1="$m1" max1="$m2" min2="$m3" max2="$m4" "$TMP1" | \
194
 
            d.graph color=`echo "$colors" | cut -d' ' -f1`
195
 
 
196
 
 
197
 
        if [ "$ARGNUM" -eq 2 ] ; then
198
 
           d.graph << EOF
199
 
              size 2 2
200
 
              move 1 90
201
 
              text max: $max_layer1
202
 
              move 75 2
203
 
              text max: $max_layer2
204
 
EOF
205
 
        fi
206
 
    fi
207
 
   done
208
 
done
209
 
 
210
 
rm -f "$TMP1"
211
 
 
212
 
 
213
 
#### overrdraw a trend line, slope+offset, and R^2 value.
214
 
if [ "$GIS_FLAG_T" -eq 1 ] ; then
215
 
 
216
 
   if [ "$ARGNUM" -ne 2 ] ; then
217
 
      g.message -e 'Will only draw trend line if two map layers are given'
218
 
      exit 1
219
 
   fi
220
 
 
221
 
   # perform linear regression
222
 
   eval `r.regression.line -g map1="$GIS_OPT_LAYER2" map2="$GIS_OPT_LAYER1"`
223
 
   g.message -v "y = $b*x + $a"
224
 
   g.message -v "R^2 = `echo "$R" | awk '{printf("%.4g", $1 * $1)}'`"
225
 
 
226
 
 
227
 
   #### calc coords for trend line in map1,map2 space
228
 
   A0="$m1"
229
 
   B0=`echo "$m1 $a $b" | awk '{print ($1-$2)/$3}'`
230
 
   A1="$m2"
231
 
   B1=`echo "$m2 $a $b" | awk '{print ($1-$2)/$3}'`
232
 
 
233
 
   # scale to 0-100% space
234
 
   X0=`echo "$A0 $m1 $m2" | awk '{print ($1-$2+1) * 100.0 / ($3-$2+1)}'`
235
 
   Y0=`echo "$B0 $m3 $m4" | awk '{print ($1-$2+1) * 100.0 / ($3-$2+1)}'`
236
 
   X1=`echo "$A1 $m1 $m2" | awk '{print ($1-$2+1) * 100.0 / ($3-$2+1)}'`
237
 
   Y1=`echo "$B1 $m3 $m4" | awk '{print ($1-$2+1) * 100.0 / ($3-$2+1)}'`
238
 
   g.message -d "Raw plot line:  $X0,$Y0    $X1,$Y1"
239
 
 
240
 
   # d.graph doesn't like % coords that are >100 or <0, so we need to work backwards
241
 
   # through the scaling and linear regression formulas to recalc coords in-bounds.
242
 
   if [ `echo "$Y0" | grep -c '\-'` -ne 0 ] ; then
243
 
      g.message -d message="y-offset at x=0% is negative, recalc"
244
 
      Y0=0
245
 
      B0=`echo "$m3" | awk '{print $1 - 1}'`
246
 
      A0=`echo "$B0 $a $b" | awk '{print ($1 * $3) + $2}'`
247
 
      X0=`echo "$A0 $m1 $m2" | awk '{print 100.0 * ($1-$2+1) / ($3-$2+1)}'`
248
 
   elif [ `echo "$Y0" | awk '{if($1 > 100) {print 1} else {print 0}}'` -eq 1 ] ; then
249
 
      g.message -d message="y position at left side of plot is > 100%, recalc"
250
 
      Y0=100
251
 
      B0="$m4"
252
 
      A0=`echo "$B0 $a $b" | awk '{print ($1 * $3) + $2}'`
253
 
      X0=`echo "$A0 $m1 $m2" | awk '{print 100.0 * ($1-$2+1) / ($3-$2+1)}'`
254
 
   fi
255
 
   g.message -d "After Y0 fixes:  $X0,$Y0    $X1,$Y1"
256
 
 
257
 
   if [ `echo "$Y1" | grep -c '\-'` -ne 0 ] ; then
258
 
      g.message -d message="y-offset at x=100% is negative, recalc"
259
 
      Y1=0
260
 
      B1=`echo "$m3" | awk '{print $1 - 1}'`
261
 
      A1=`echo "$B1 $a $b" | awk '{print ($1 * $3) + $2}'`
262
 
      X1=`echo "$A1 $m1 $m2" | awk '{print 100.0 * ($1-$2+1) / ($3-$2+1)}'`
263
 
   elif [ `echo "$Y1" | awk '{if($1 > 100) {print 1} else {print 0}}'` -eq 1 ] ; then
264
 
      g.message -d message="y position at right side of plot is > 100%, recalc"
265
 
      Y1=100
266
 
      B1="$m4"
267
 
      A1=`echo "$B1 $a $b" | awk '{print ($1 * $3) + $2}'`
268
 
      X1=`echo "$A1 $m1 $m2" | awk '{print 100.0 * ($1-$2+1) / ($3-$2+1)}'`
269
 
   fi
270
 
 
271
 
   g.message -d "final line coords: ($X0, $Y0) to ($X1, $Y1)"
272
 
 
273
 
 
274
 
   # plot the trend line and coefficient text
275
 
   d.graph << EOF
276
 
      color red
277
 
      width 2
278
 
      move $X0 $Y0
279
 
      draw $X1 $Y1
280
 
      width 0
281
 
 
282
 
      color indigo
283
 
      size 2
284
 
      move 5 80
285
 
      text y = $b*x + $a
286
 
      move 5 75
287
 
      text R^2 = `echo "$R" | awk '{printf("%.4g", $1 * $1)}'`
288
 
# um, r.regression.line's R is the R of R^2 fame, right??
289
 
EOF
290
 
fi
291