# Can anyone shed any light on how the Shaded Relief Map of Spearfish
# was produced, which appears on the abstract page of the document,
# Introduction to GRASS 4 by James Westerveldt, June 11, 1991?
#
# +---------------------------------------------------------------------+
# | John Miller | Internet Mail - ccjwm@jcu.edu.au |
This map was created with the GRASS shell script shade.rel.sh - distributed
with 4.0.
A slightly different version of this script created recently allows 3
light sources to be defined - red, green, blue. The script (I call it
shade.clr.sh) follows:
--------------------- CUT HERE ---------------------
#!/bin/sh
# set nsres and ewres
eval `g.region -g`
echo ""
echo Please provide the altitude of the sun in degrees above the
echo horizon and the azimuth of the red, green, and blue lights
echo in degrees to the east of 'north (N:0 E:90 S:180 W:270)'
echo Might we suggest 60 for red, 180 for green and 300 for blue.
echo ""
gotit=0
while test $gotit -eq 0
do
echo -n "altitude: "
read alt
if test $alt -gt 0 -a $alt -lt 90
then
gotit=1
else
echo Sorry, altitude must be greater than 0 and less than 90
fi
done
gotit=0
while test $gotit -eq 0
do
echo -n "red light azimuth: "
read raz
if test $raz -ge 0 -a $raz -lt 360
then
gotit=1
else
echo Sorry, azimuth must be greater than -1 and less than 360
fi
done
gotit=0
while test $gotit -eq 0
do
echo -n "green light azimuth: "
read gaz
if test $gaz -ge 0 -a $gaz -lt 360
then
gotit=1
else
echo Sorry, azimuth must be greater than -1 and less than 360
fi
done
gotit=0
while test $gotit -eq 0
do
echo -n "blue light azimuth: "
read baz
if test $baz -ge 0 -a $baz -lt 360
then
gotit=1
else
echo Sorry, azimuth must be greater than -1 and less than 360
fi
done
echo ""
g.ask type=old element=cell desc=raster prompt="Enter elevation file" unixfile=/tmp/$$
eval `cat /tmp/$$`
rm -f /tmp/$$
if [ ! "$file" ]
then
exit 0
fi
elev="'${fullname}'"
echo "$elev"
raz=`expr $raz - 90`
gaz=`expr $gaz - 90`
baz=`expr $baz - 90`
echo ""
echo Running r.mapcalc, please stand by.
echo Your new map will be named shade. Please consider renaming.
echo ""
r.mapcalc << EOF
shade = eval( \\
x=($elev[-1,-1] + 2*$elev[0,-1] + $elev[1,-1] \\
-$elev[-1,1] - 2*$elev[0,1] - $elev[1,1])/(8.*$ewres) , \\
y=($elev[-1,-1] + 2*$elev[-1,0] + $elev[-1,1] \\
-$elev[1,-1] - 2*$elev[1,0] - $elev[1,1])/(8.*$nsres) , \\
slope=90.-atan(sqrt(x*x + y*y)), \\
a=round(atan(x,y)), \\
aspect=if(x||y,if(a,a,360)), \\
rang = sin($alt)*sin(slope) + cos($alt)*cos(slope) * cos($raz-aspect), \\
red = int(if(rang < 0,0,4.9*rang)), \\
gang = sin($alt)*sin(slope) + cos($alt)*cos(slope) * cos($gaz-aspect), \\
green = int(if(gang < 0,0,4.9*gang)), \\
bang = sin($alt)*sin(slope) + cos($alt)*cos(slope) * cos($baz-aspect), \\
blue = int(if(bang < 0,0,4.9*bang)), \\
1 + red + 5 * green + 25 * blue )
EOF
r.colors shade color=grey
echo ""
echo shade.relief map created and named shade. Consider renaming
r.colors shade color=rules 2>&1 > /dev/null << EOF
1 0 0 0
5 255 0 0
6 0 64 0
10 255 64 0
11 0 128 0
15 255 128 0
16 0 192 0
20 255 192 0
21 0 255 0
25 255 255 0
26 0 0 64
30 255 0 64
31 0 64 64
35 255 64 64
36 0 128 64
40 255 128 64
41 0 192 64
45 255 192 64
46 0 255 64
50 255 255 64
51 0 0 128
55 255 0 128
56 0 64 128
60 255 64 128
61 0 128 128
65 255 128 128
66 0 192 128
70 255 192 128
71 0 255 128
75 255 255 128
76 0 0 192
80 255 0 192
81 0 64 192
85 255 64 192
86 0 128 192
90 255 128 192
91 0 192 192
95 255 192 192
96 0 255 192
100 255 255 192
101 0 0 255
105 255 0 255
106 0 64 255
110 255 64 255
111 0 128 255
115 255 128 255
116 0 192 255
120 255 192 255
121 0 255 255
125 255 255 255
EOF
--------------------- CUT HERE ---------------------