[GRASS-user] script for a multidirectional, oblique-weighted, shaded-relief image

[maning sambale]

> FOCALMEAN is just a smoothing kernel using the mean value of a 3x3 neighborhood
>
> 3 3 2
> 1 4 8 = (3 + 3 + 2 + 1 + 4 + 8 + 9 + 5 + 3) / 9 = 4.2
> 9 5 3
>
> I think you can get the same effect by using r.neighbors:
>
> r.neighbors in=<ingrid> out=<outgrid> method=average size=3

Thanks!

Here is the script I made.  This is my second attempt to create a
script so please bear with me if it might be a bit crappy script.  I
created this based from an arc grid formula from
http://pubs.usgs.gov/of/1992/of92-422/of92-422.pdf  I'm not sure if I
made them correctly.  Comments are very much welcome.  Do you think it
can be useful to others?

You can find the page for the script here:
http://esambale.wikispaces.com/hillshade+multi

The script is below
#!/bin/sh
############################################################################
#
# MODULE:       r.hillshade.multi
#
# AUTHOR(S):    Emmanuel Sambale esambale at yahoo.com emmanuel.sambale at 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:
#               http://pubs.usgs.gov/of/1992/of92-422/of92-422.pdf
#
# 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.
#%End
#%flag
#%  key: r
#%  description: remove temporary maps
#%END
#%option
#% key: input
#% type: string
#% gisprompt: old,cell,raster
#% description: Name of DEM map
#% required : yes
#%End
#%option
#% key: window
#% type: integer
#% description: window size for aspect resampling default is 5
#% 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 ""
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"
echo ""
r.neighbors in=$GIS_OPT_input out=res1 method=average size="$GIS_OPT_window"

# compute aspect map
echo ""
echo "Computing aspect map..."
echo ""
r.slope.aspect elevation=res1 format=degrees prec=float aspect=aspect
zfactor=1.0 min_slp_allowed=0.0

# 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
  fi

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

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

#####################################################################

Cheers,

Maning