----------
X-Sun-Data-Type: text
X-Sun-Data-Description: text
X-Sun-Data-Name: text
X-Sun-Content-Lines: 6


Here is a slightly enhanced shade.rel.sh.  This version
gives you the possibility to adjust the "vertical exaggeration",
and also gives the output a slightly more sensible name.

Simon Cox
----------
X-Sun-Data-Type: shell-script
X-Sun-Data-Description: shell-script
X-Sun-Data-Name: shade.rel.sh
X-Sun-Content-Lines: 106

:

# 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 sun in degrees to the east of
echo 'north (N:0 E:90 S:180 W:270)'
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 "azimuth: "
	read az
	if test $az -ge 0 -a $az -lt 360
	then
		gotit=1
	else
		echo Sorry, azimuth must be greater than -1  and less than 360
	fi
done

echo ""
echo "Please provide a scale factor for vertical exaggeration  "
echo "use a negative number for a compression factor"
echo ""
echo "Use the vertical exaggeration also to compensate for maps"
echo "where the horizontal and vertical scales are not comparable"
echo ""

gotit=0
while test $gotit -eq 0
do
	echo -n "scale factor:  "
	read scale
	if test $scale -ne 0
	then
		gotit=1
	else
		echo Sorry, scale factor must be non-zero
	fi
done

if test $scale -gt 0
then
	ewres=`expr $ewres/$scale`
	nsres=`expr $nsres/$scale`
else
	ewres=`expr -1*$scale*$ewres`
	nsres=`expr -1*$scale*$nsres`
fi

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}"
shade="${name}.shade"

echo "$elev"

az=`expr $az - 90`

# echo $az $alt $ewres $nsres

echo ""
echo Running r.mapcalc, please stand by.
echo Your new map will be named $shade.
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)), \\
 cang = sin($alt)*sin(slope) + cos($alt)*cos(slope) * cos($az-aspect), \\
 if(cang < 0,0,100*cang))
EOF

r.colors $shade color=grey

echo ""
echo shade.relief map created and named $shade.
----
______________________________________________________________________________
from May 1st (or thereabouts)        __  L     Dr Simon Cox                   
CSIRO Div. Exploration & Mining   ,~'  L_|\    VIEPS Dept of Earth Sciences,  
39 Fairway, PO Box 437,        ,-'         \   Monash University, Clayton     
Nedlands, WA  6009,  Australia (   soon     \  Vic 3168,  Australia           
      Phone +61 9 389 8421     X<~~~~~~~~~  /         Phone +61 3 905 5762    
      Fax   +61 9 389 1906      L,~'   "\_x/          Fax   +61 3 905 5062    
    s.cox@dem.csiro.au                    u  simon@artemis.earth.monash.edu.au
______________________________________________________________________________

