#!/bin/sh -x
############################################################################
#
# MODULE:       r.hillshade.multi
#
# AUTHOR(S):    Emmanuel Sambale esambale@yahoo.com emmanuel.sambale@gmail.com
#               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

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.shaded.relief map=$GIS_OPT_input shadedmap=shade225 altitude=30 azimuth=225 zmult=1 scale=1
r.shaded.relief map=$GIS_OPT_input shadedmap=shade270 altitude=30 azimuth=270 zmult=1 scale=1
r.shaded.relief map=$GIS_OPT_input shadedmap=shade315 altitude=30 azimuth=315 zmult=1 scale=1
r.shaded.relief map=$GIS_OPT_input shadedmap=shade360 altitude=30 azimuth=360 zmult=1 scale=1

# 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 prec=float aspect=aspect1 zfactor=1.0 min_slp_allowed=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 = sin((aspect - 225)*sin(aspect - 225))"
r.mapcalc "w270 = sin((aspect - 270)*sin(aspect - 270))"
r.mapcalc "w315 = sin((aspect - 315)*sin(aspect - 315))"
r.mapcalc "w360 = sin((aspect)*sin(aspect))"


#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