← Back to branch summary

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

« back to all changes in this revision

Viewing changes to scripts/r.plane/r.plane

  • 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/r.plane/r.plane
1
 
#!/bin/sh
2
 
#
3
 
############################################################################
4
 
#
5
 
# MODULE:       r.plane for GRASS 5.7; based on r.plane for GRASS 5
6
 
# AUTHOR(S):    CERL?; updated to GRASS 5.7 by Michael Barton
7
 
#               Dec 2004: Alessandro Frigeri & Ivan Marchesini
8
 
#               Modified to produce floating and double values maps
9
 
# PURPOSE:      Creates a raster plane map from user specified inclination and azimuth
10
 
# COPYRIGHT:    (C) 2004 by the GRASS Development Team
11
 
#
12
 
#               This program is free software under the GNU General Public
13
 
#               License (>=v2). Read the file COPYING that comes with GRASS
14
 
#               for details.
15
 
#
16
 
#############################################################################
17
 
 
18
 
#%Module
19
 
#%  description: Creates raster plane map given dip (inclination), aspect (azimuth) and one point.
20
 
#%  keywords: raster, elevation
21
 
#%End
22
 
#%option
23
 
#% key: name
24
 
#% type: string
25
 
#% gisprompt: new,cell,raster
26
 
#% description: Name of raster plane to be created
27
 
#% required : yes
28
 
#%end
29
 
#%option
30
 
#% key: dip
31
 
#% type: double
32
 
#% gisprompt: -90-90
33
 
#% answer: 0.0
34
 
#% description: Dip of plane. Value must be between -90 and 90 degrees
35
 
#% required : yes
36
 
#%end
37
 
#%option
38
 
#% key: azimuth
39
 
#% type: double
40
 
#% gisprompt: 0-360
41
 
#% answer: 0.0
42
 
#% description: Azimuth of the plane. Value must be between 0 and 360 degrees
43
 
#% required : yes
44
 
#%end
45
 
#%option
46
 
#% key: easting
47
 
#% type: double
48
 
#% description: Easting coordinate of a point on the plane
49
 
#% required : yes
50
 
#%end
51
 
#%option
52
 
#% key: northing
53
 
#% type: double
54
 
#% description: Northing coordinate of a point on the plane
55
 
#% required : yes
56
 
#%end
57
 
#%option
58
 
#% key: elevation
59
 
#% type: double
60
 
#% description: Elevation coordinate of a point on the plane
61
 
#% required : yes
62
 
#%end
63
 
#%option
64
 
#% key: type
65
 
#% type: string 
66
 
#% options: int,float,double
67
 
#% description: Type of the raster map to be created
68
 
#% required : yes
69
 
#%end
70
 
 
71
 
if  [ -z "$GISBASE" ] ; then
72
 
 echo "You must be in GRASS GIS to run this program." >&2
73
 
 exit 1
74
 
fi   
75
 
 
76
 
if [ "$1" != "@ARGS_PARSED@" ] ; then
77
 
  exec g.parser "$0" "$@"
78
 
fi
79
 
 
80
 
PROG=`basename "$0"`
81
 
 
82
 
#### check if we have awk
83
 
if [ ! -x "`which awk`" ] ; then
84
 
    g.message -e "awk required, please install awk or gawk first"
85
 
    exit 1
86
 
fi
87
 
 
88
 
#### setup temporary file
89
 
TEMPFILE="`g.tempfile pid=$$`"
90
 
if [ $? -ne 0 ] || [ -z "$TEMPFILE" ] ; then
91
 
    g.message -e "unable to create temporary files" 
92
 
    exit 1
93
 
fi
94
 
 
95
 
#### trap ctrl-c so that we can clean up tmp
96
 
trap 'rm -f "$TEMPFILE"*' 2 3 15
97
 
 
98
 
# setting environment, so that awk works properly in all languages
99
 
unset LC_ALL
100
 
LC_NUMERIC=C
101
 
export LC_NUMERIC
102
 
 
103
 
### setup enviro vars ###
104
 
dip="$GIS_OPT_DIP"
105
 
az="$GIS_OPT_AZIMUTH"
106
 
ea="$GIS_OPT_EASTING"
107
 
no="$GIS_OPT_NORTHING"
108
 
el="$GIS_OPT_ELEVATION"
109
 
name="$GIS_OPT_NAME"
110
 
type="$GIS_OPT_TYPE"
111
 
 
112
 
eval `g.region -g`
113
 
 
114
 
 
115
 
### test input values ###
116
 
diptest=`echo "$dip" | awk '{ printf("%8d", int($1 + 0.5))'}`
117
 
if [ "$diptest" -lt -90 -o "$diptest" -gt 90 ]
118
 
then
119
 
        g.message -e "Sorry, dip must be greater than -90 and less than 90.\
120
 
          Please enter a valid value."
121
 
        exit 1
122
 
fi
123
 
 
124
 
aztest=`echo "$az" | awk '{ printf("%8d", int($1 + 0.5))'}`
125
 
if  [ "$aztest" -lt 0 -o "$aztest" -ge 360 ]
126
 
then
127
 
        g.message -e "Sorry, azimuth must be no less than 0 and less than 360"
128
 
        exit 1
129
 
fi
130
 
 
131
 
if  [ $type = "int" ]
132
 
then
133
 
        g.message "Preparing to produce a CELL map.."
134
 
        typeflag=''
135
 
        ctype=cell
136
 
fi
137
 
 
138
 
if  [ $type = "float" ]
139
 
then
140
 
        g.message "Preparing to produce an FCELL map.."
141
 
        typeflag=-f
142
 
        ctype=fcell
143
 
fi
144
 
 
145
 
if  [ $type = "double" ]
146
 
then
147
 
        g.message "Preparing to produce a DCELL map.."
148
 
        typeflag=-d
149
 
        ctype=dcell
150
 
fi
151
 
 
152
 
# these tests will fail for negative lat/lon?
153
 
eaint=`echo "$ea" | awk '{ printf "%d",int($1 + 0.5) }'`
154
 
if test "$eaint" -le `echo "$e" | cut -d. -f1` -a "$eaint" -ge `echo "$w" | cut -d. -f1`
155
 
then
156
 
        gotit=1
157
 
#       echo "east: $ea"
158
 
else
159
 
        g.message -e "Sorry, point must be within current region"
160
 
        g.message -e "Current region:"
161
 
        g.message -e "west: $w   east: $e"
162
 
        exit 1
163
 
fi
164
 
 
165
 
noint=`echo "$no" | awk '{ printf "%d", int($1 + 0.5) }'`
166
 
if test "$noint" -gt `echo "$s" | cut -d. -f1` -a "$noint" -lt `echo "$n" | cut -d. -f1`
167
 
then
168
 
        gotit=1
169
 
#       echo "north: $no"
170
 
else
171
 
        g.message -e "Sorry, point must be within current region"
172
 
        g.message -e "Current region:"
173
 
        g.message -e "south: $s   north: $n"
174
 
        exit 1
175
 
fi
176
 
 
177
 
 
178
 
### now the actual algorithm in awk (stored in a temporary file) ###
179
 
cat > "$TEMPFILE" << EOF
180
 
 
181
 
{
182
 
if (NR==1) {
183
 
# print file header
184
 
  rows = (north-south) / nsres
185
 
  cols = (east-west) / ewres
186
 
    printf("east:     %d\n",east)
187
 
    printf("west:     %d\n",west)
188
 
    printf("south:    %d\n",south)
189
 
    printf("north:    %d\n",north)
190
 
    printf("cols:     %d\n",cols)
191
 
    printf("rows:     %d\n",rows)
192
 
    cells=rows*cols
193
 
    z=1
194
 
    }
195
 
 
196
 
if (NR==2) {
197
 
    pi=3.14159265358979323846
198
 
    a2=(az*pi)/180
199
 
    dip2=(dip*pi)/180 
200
 
    tandip=(sin(dip2)/cos(dip2))
201
 
    northc=north-(0.5*nsres)
202
 
    southc=south+(0.5*nsres)
203
 
    eastc=east-(0.5*ewres)
204
 
    westc=west+(0.5*ewres)
205
 
   
206
 
    for (y=northc; y >= southc; y=y-nsres) {
207
 
       for (x=westc; x <= eastc; x=x+ewres) {
208
 
         dx=(ea-x)
209
 
         dy=(y-no)
210
 
         dist = sqrt((dx*dx) + (dy*dy))
211
 
         if (dist==0) {
212
 
           new_elev[z]=el
213
 
           }
214
 
         else {
215
 
           gamma = atan2((dx/dist),(dy/dist))
216
 
           epsilon=a2-gamma
217
 
           d=dist*cos(epsilon)
218
 
           h=(d*sin(dip2)/cos(dip2))    
219
 
           #print h
220
 
           new_elev[z]=el-h
221
 
           z++
222
 
           }
223
 
         }
224
 
     }
225
 
}
226
 
if (NR>=2) {
227
 
}
228
 
}
229
 
END {  
230
 
      for (z=1; z <= cells; z++) {
231
 
        if (typeflag==""){
232
 
          printf " %d",new_elev[z]
233
 
          }
234
 
        else {
235
 
          printf " %f",new_elev[z]
236
 
          }
237
 
        }
238
 
}
239
 
 
240
 
EOF
241
 
 
242
 
 
243
 
### execute awk and remove temporary file ###
244
 
 
245
 
awk  -f "$TEMPFILE" -v east=$e west=$w north=$n south=$s ea=$ea no=$no typeflag=$typeflag \
246
 
  nsres=$nsres ewres=$ewres  dip=$dip az=$az el=$el "$TEMPFILE" > "$TEMPFILE".2
247
 
 
248
 
rm -f "$TEMPFILE"
249
 
 
250
 
g.message "Running r.in.ascii, please stand by.."
251
 
r.in.ascii `echo $typeflag` i="$TEMPFILE".2 o=$name
252
 
 
253
 
r.support "$name" hist="$PROG name=$name dip=$dip azimuth=$az easting=$ea \\"
254
 
r.support "$name" hist="northing=$no elevation=$el type=$type"
255
 
 
256
 
g.message "Done."
257
 
g.message message="Raster map <$name> generated by r.plane \
258
 
  at point $ea E, $no N, elevation $el with dip=$dip degrees and \
259
 
  aspect=$az degrees ccw from north."
260
 
#echo "$dat, user: $user"
261
 
 
262
 
rm -f "$TEMPFILE".2
263