← Back to branch summary

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

« back to all changes in this revision

Viewing changes to scripts/r.shaded.relief/r.shaded.relief

  • 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.shaded.relief/r.shaded.relief
1
 
#!/bin/sh
2
 
 
3
 
############################################################################
4
 
#
5
 
# MODULE:       r.shaded.relief
6
 
# AUTHOR(S):    CERL
7
 
#               updates: Michael Barton, 2004
8
 
#               updates: Gordon Keith, 2003
9
 
#               updates: Andreas Lange, 2001
10
 
#               updates: David Finlayson, 2001
11
 
#               updates: Markus Neteler, 2001, 1999
12
 
# PURPOSE:      Creates shaded relief map from raster elevation map (DEM)
13
 
# COPYRIGHT:    (C) 1999 - 2004 by the GRASS Development Team
14
 
#
15
 
#               This program is free software under the GNU General Public
16
 
#               License (>=v2). Read the file COPYING that comes with GRASS
17
 
#               for details.
18
 
#
19
 
#############################################################################
20
 
#
21
 
#   July 2007 - allow input from other mapsets (Brad Douglas)
22
 
#
23
 
#   May 2005 - fixed wrong units parameter (Markus Neteler)
24
 
#
25
 
#   September 2004 - Added z exaggeration control (Michael Barton) 
26
 
#   April 2004 - updated for GRASS 5.7 by Michael Barton 
27
 
#
28
 
#   9/2004 Adds scale factor input (as per documentation); units set scale only if specified for lat/long regions
29
 
#    Also, adds option of controlling z-exaggeration.
30
 
#
31
 
#   6/2003 fixes for Lat/Long Gordon Keith <gordon.keith@csiro.au>
32
 
#   If n is a number then the ewres and nsres are mulitplied by that scale
33
 
#    to calculate the shading.
34
 
#   If n is the letter M (either case) the number of metres is degree of
35
 
#    latitude is used as the scale.
36
 
#   If n is the letter f then the number of feet in a degree is used.
37
 
#   It scales latitude and longitude equally, so it's only approximately
38
 
#   right, but for shading its close enough. It makes the difference
39
 
#   between an unusable and usable shade.
40
 
#
41
 
#   10/2001 fix for testing for dashes in raster file name
42
 
#        by Andreas Lange <andreas.lange@rhein-main.de>
43
 
#   10/2001 added parser support - Markus Neteler
44
 
#   9/2001 fix to keep NULLs as is (was value 22 before) - Markus Neteler
45
 
#   1/2001 fix for NULL by David Finlayson <david_finlayson@yahoo.com>
46
 
#   11/99 updated $ewres to ewres() and $nsres to nsres()
47
 
#       updated number to FP in r.mapcalc statement Markus Neteler
48
 
 
49
 
#% Module
50
 
#%  description: Creates shaded relief map from an elevation map (DEM).
51
 
#%  keywords: raster, elevation
52
 
#% End
53
 
#% option
54
 
#%  key: map
55
 
#%  type: string
56
 
#%  gisprompt: old,cell,raster
57
 
#%  description: Input elevation map
58
 
#%  required : yes
59
 
#% end
60
 
#% option
61
 
#%  key: shadedmap
62
 
#%  type: string
63
 
#%  gisprompt: new,cell,raster
64
 
#%  description: Output shaded relief map name
65
 
#%  required : no
66
 
#% end
67
 
#% option
68
 
#%  key: altitude
69
 
#%  type: double
70
 
#%  description: Altitude of the sun in degrees above the horizon
71
 
#%  required : no
72
 
#%  options : 0-90
73
 
#%  answer: 30
74
 
#% end
75
 
#% option
76
 
#%  key: azimuth
77
 
#%  type: double
78
 
#%  description: Azimuth of the sun in degrees to the east of north
79
 
#%  required : no
80
 
#%  options : 0-360
81
 
#%  answer: 270
82
 
#% end
83
 
#% option
84
 
#%  key: zmult    
85
 
#%  type: double
86
 
#%  description: Factor for exaggerating relief
87
 
#%  required : no
88
 
#%  answer: 1
89
 
#% end
90
 
#% option
91
 
#%  key: scale
92
 
#%  type: double
93
 
#%  description: Scale factor for converting horizontal units to elevation units
94
 
#%  required : no
95
 
#%  answer: 1
96
 
#% end
97
 
#% option
98
 
#%  key: units
99
 
#%  type: string
100
 
#%  description: Set scaling factor (applies to lat./long. locations only, none: scale=1)
101
 
#%  required : no
102
 
#%  options: none,meters,feet
103
 
#%  answer: none
104
 
#% end
105
 
 
106
 
if  [ -z "$GISBASE" ] ; then
107
 
    echo "You must be in GRASS GIS to run this program." >&2
108
 
 exit 1
109
 
fi
110
 
 
111
 
# save command line
112
 
if [ "$1" != "@ARGS_PARSED@" ] ; then
113
 
    CMDLINE=`basename "$0"`
114
 
    for arg in "$@" ; do
115
 
        CMDLINE="$CMDLINE \"$arg\""
116
 
    done
117
 
    export CMDLINE
118
 
    exec g.parser "$0" "$@"
119
 
fi
120
 
 
121
 
# set nsres and ewres (unneeded because functions in r.mapcalc)
122
 
#nsres=`g.region -g | grep 'nsres=' | sed s/nsres=//`
123
 
#ewres=`g.region -g | grep 'ewres=' | sed s/ewres=//`
124
 
 
125
 
PROG=`basename "$0"`
126
 
 
127
 
# setting environment, so that awk works properly in all languages
128
 
unset LC_ALL
129
 
LC_NUMERIC=C
130
 
export LC_NUMERIC
131
 
 
132
 
alt="$GIS_OPT_ALTITUDE"
133
 
az="$GIS_OPT_AZIMUTH"
134
 
elev="$GIS_OPT_MAP"
135
 
zmult="$GIS_OPT_ZMULT"
136
 
scale="$GIS_OPT_SCALE"
137
 
 
138
 
eval `g.findfile element=cell file=$elev`
139
 
if [ -z "$name" ] ; then
140
 
   g.message -e  "Map <$elev> not found."
141
 
   exit 1
142
 
fi
143
 
 
144
 
if [ "$GIS_OPT_MAP" = "$GIS_OPT_SHADEDMAP" ]; then
145
 
        g.message -e "Input elevation map and output relief map must have different names"
146
 
        exit 1
147
 
fi
148
 
 
149
 
if [ -z "$GIS_OPT_SHADEDMAP" ]; then
150
 
    ELEVSHADE="${elev%%@*}.shade"
151
 
else
152
 
    ELEVSHADE="$GIS_OPT_SHADEDMAP"
153
 
fi
154
 
 
155
 
# test for overwrite as r.mapcalc doesn't
156
 
eval `g.findfile element=cell file="$ELEVSHADE" mapset=.`
157
 
if [ -n "$name" ] ; then
158
 
   # map already exists
159
 
   if [ -z "$GRASS_OVERWRITE" ] || [ "$GRASS_OVERWRITE" -ne 1 ] ; then
160
 
      g.message -e "option <shadedmap>: <$ELEVSHADE> exists."
161
 
      exit 1
162
 
   fi
163
 
fi
164
 
 
165
 
# allow - and . chars in the map name
166
 
elev="\"$elev\""
167
 
elev_out="\"$ELEVSHADE\""
168
 
 
169
 
#LatLong locations only:
170
 
if [ "$GIS_OPT_UNITS" = "meters" ] ; then
171
 
   #scale=111120
172
 
   scale=`echo $scale | awk '{printf("%f", $1 * 1852*60 )}'`
173
 
fi
174
 
 
175
 
#LatLong locations only:
176
 
if [ "$GIS_OPT_UNITS" = "feet" ] ; then
177
 
   #scale=364567.2
178
 
   scale=`echo $scale | awk '{printf("%f", $1 * 6076.12*60 )}'`
179
 
fi              
180
 
 
181
 
#correct azimuth to East (GRASS convention):
182
 
#  this seems to be backwards, but in fact it works so leave it.
183
 
az=`echo "$az" | awk '{print $1 - 90.0}'`
184
 
 
185
 
g.message "Calculating shading, please stand by."
186
 
 
187
 
# Note: no space allowed after \\:
188
 
 
189
 
r.mapcalc << EOF
190
 
$elev_out = eval( \\
191
 
 x=($zmult*$elev[-1,-1] + 2*$zmult*$elev[0,-1] + $zmult*$elev[1,-1] \\
192
 
   -$zmult*$elev[-1,1] - 2*$zmult*$elev[0,1] - $zmult*$elev[1,1])/(8.*ewres()*$scale) , \\
193
 
 y=($zmult*$elev[-1,-1] + 2*$zmult*$elev[-1,0] + $zmult*$elev[-1,1] \\
194
 
   -$zmult*$elev[1,-1] - 2*$zmult*$elev[1,0] - $zmult*$elev[1,1])/(8.*nsres()*$scale) , \\
195
 
 slope=90.-atan(sqrt(x*x + y*y)), \\
196
 
 a=round(atan(x,y)), \\
197
 
 a=if(isnull(a),1,a), \\
198
 
 aspect=if(x!=0||y!=0,if(a,a,360.)), \\
199
 
 cang = sin($alt)*sin(slope) + cos($alt)*cos(slope) * cos($az-aspect), \\
200
 
 if(cang < 0.,0.,100.*cang), \\
201
 
 if(isnull(cang), null(), 100.*cang))
202
 
EOF
203
 
 
204
 
if [ $? -ne 0 ] ; then 
205
 
   g.message -e "In calculation, script aborted."
206
 
   exit 1
207
 
fi
208
 
 
209
 
r.colors "$ELEVSHADE" color=grey
210
 
 
211
 
# record metadata
212
 
r.support "$ELEVSHADE" title="Shaded relief of \"$GIS_OPT_MAP\"" history=""
213
 
r.support "$ELEVSHADE" history="r.shaded.relief settings:"
214
 
r.support "$ELEVSHADE" history="  altitude=$alt  azimuth=$GIS_OPT_AZIMUTH  zmult=$zmult  scale=$GIS_OPT_SCALE"
215
 
if [ -n "$GIS_OPT_UNITS" ] ; then
216
 
   r.support map="$ELEVSHADE" history="  units=$GIS_OPT_UNITS  (adjusted scale=$scale)"
217
 
fi
218
 
 
219
 
if [ -n "$GIS_OPT_SHADEDMAP" ] ; then
220
 
    g.message "Shaded relief map created and named <$ELEVSHADE>."
221
 
else
222
 
    g.message "Shaded relief map created and named <$ELEVSHADE>. Consider renaming."
223
 
fi
224
 
 
225
 
# write cmd history:
226
 
r.support $ELEVSHADE history="${CMDLINE}"
227
 
 
228
 
exit 0