Creates a multidirectional, oblique-weighted, shaded-relief image from an input DEM image.

r.hillshade.multi

#!/bin/sh -x
############################################################################
#
# MODULE:       r.hillshade.multi
#
# AUTHOR(S):    Emmanuel Sambale [email protected] [email protected]
#               with comments and improvement from the GRASS user mailing list.
#
# PURPOSE:      Creates a multidirectional, oblique-weighted, shaded-relief image
#               from an input DEM image.  Original ARC 6.0.1 GRID procedure was 
#               from Mark, Robert. 1992. Multidirectional, oblique-weighted, 
#               shaded-relief image of the Island of Hawaii. U.S. Geological 
#               Survey. Menlo Park, CA 94025. Available online:  
#
# COPYRIGHT:    (C) 2006 by the GRASS Development Team
#
#               This program is free software under the GNU General Public
#               License (>=v2). Read the file COPYING that comes with GRASS
#               for details.
#
#############################################################################

#%Module
#% description: Creates a multidirectional, oblique-weighted, shaded-relief image from an input DEM image.  The input DEM should cropped to the area of interest only.
#%End
#%flag
#%  key: r
#%  description: Remove temporary maps
#%END
#%option
#% key: input
#% type: string
#% gisprompt: old,cell,raster
#% description: Name of DEM
#% required : yes
#%End
#%option
#% key: window
#% type: integer
#% description: Window size for aspect resampling default is 5.  Window size must be odd.
#% answer : 5
#% required : yes
#%END
#%option
#% key: zscale
#% type: integer
#% description: Scale factor for converting horizontal units to vertical.
#% answer: 1 
#% required : yes
#%END

if  [ -z "$GISBASE" ]
then
	echo ""
	echo "You must be in GRASS GIS to run this program"
	echo ""
	exit 1
fi

if   [ "$1" != "@ARGS_PARSED@" ]
then
	exec g.parser "$0" "$@"
fi


unset LC_ALL
export LC_NUMERIC=C

#set to current input map region
g.region rast=$GIS_OPT_input -p

# Compute hillshade at azimuth 225, 270, 315 and 360 at 30 degrees sun illumination angle
echo ""
echo "Computing hillshade at azimuth 225, 270, 315 and 360 at 30 degrees sun illumination angle ..."
echo "The azimuth of the sun in degrees to the east of north (a value between 0 and 360 degrees) "
r.relief input=$GIS_OPT_input output=shade225 altitude=30 azimuth=225 scale=$GIS_OPT_zscale
r.relief input=$GIS_OPT_input output=shade270 altitude=30 azimuth=270 scale=$GIS_OPT_zscale
r.relief input=$GIS_OPT_input output=shade315 altitude=30 azimuth=315 scale=$GIS_OPT_zscale
r.relief input=$GIS_OPT_input output=shade360 altitude=30 azimuth=360 scale=$GIS_OPT_zscale

# Resample DEM map
echo ""
echo "Resampling DEM map from the given window size $GIS_OPT_window"
echo ""
#g.region res=100
r.resample in=$GIS_OPT_input out=res1000
#g.region res=30
r.neighbors in=res1000 out=res1 method=average size="$GIS_OPT_window"

# compute aspect map
echo ""
echo "Computing aspect map..."
echo "Aspect is calculated counterclockwise from east"
r.slope.aspect elevation=res1 format=degrees aspect=aspect1 min_slope=0.0
#"re-orient aspect"
r.mapcalc "aspect = aspect1 + 90"
# compute weights 225, 270, 315 and 360

echo ""
echo "computing weights..."
echo ""
r.mapcalc "w225 = (cos(aspect - 225) + 1) / 2"
r.mapcalc "w270 = (cos(aspect - 270) + 1) / 2"
r.mapcalc "w315 = (cos(aspect - 315) + 1) / 2"
r.mapcalc "w360 = (cos(aspect) + 1) / 2"

#compute weighted
r.mapcalc "weightedshade = (((w225 * shade225) + (w270 * shade270) + (w315 * shade315) + (w360 * shade360))/2)"

if [ $GIS_FLAG_r -eq 1 ] ; then
  echo ""
  echo "Temporary files deleted ...."
  g.remove rast=shade225,shade270,shade315,shade360,w225,w270,w315,w360,aspect,res1,res1000
  fi

echo ""
echo "Weighted hillshade image created raster name weightedshade."
echo ""

#display
r.colors map=weightedshade color=grey
d.mon select=x0
d.erase
d.rast weightedshade