A new shell script combines the benefits of shade-relief maps (see the
GRASS program shade.rel.sh) and aspect maps. Shade-relief maps provide
for the subtle shadings which we naturally perceive as depth. A drawback
is that depending on the placement of the illumination source, the eye and
brain see different things. Displaying an aspect map using a color table
that 1) gradually changes color and 2) begins and ends at the same color -
provides the viewer with better information about the aspect direction. It
has the unfortunate drawback of removing all slope information. That is,
a 0.001% slope will be represented identically to a 99% slope.
A new grass shell script allows for the simultaneous generation and
combination of 3 shade-relief maps. One is illuminated with red, another
green, and another blue. The user chooses the identical altitude for the
lights and independent azimuths (compas headings). At CERL the new script,
which uses r.mapcalc, is called shade.clr.sh.
For remote sites, the script 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