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

GRASS GRASS user
1

Hello,

I found an article creating a multidirectional, oblique-weighted,
shaded-relief image from the USGS
(http://pubs.usgs.gov/of/1992/of92-422/of92-422.pdf).

I basically creates a hillshade using different illumination angles.
This is especially useful in very mountainous areas. I tried to
recreate this with GRASS but having difficulty with r.mapcalc
equations.

Heres what I've made so far:

# Multidirectional, oblique-weighted, shaded-relief image

# Compute hillshade at azimuth 225, 270, 315 and 360 at 30 degrees sun
illumination angle
r.shaded.relief map=DEM_90@PERMANENT shadedmap=shade225 altitude=30
azimuth=225 zmult=1 scale=1
r.shaded.relief map=DEM_90@PERMANENT shadedmap=shade270 altitude=30
azimuth=270 zmult=1 scale=1
r.shaded.relief map=DEM_90@PERMANENT shadedmap=shade315 altitude=30
azimuth=315 zmult=1 scale=1
r.shaded.relief map=DEM_90@PERMANENT shadedmap=shade360 altitude=30
azimuth=360 zmult=1 scale=1

# compute aspect map
r.slope.aspect elevation=DEM_90@PERMANENT format=degrees prec=float
aspect=aspect zfactor=1.0 min_slp_allowed=0.0

# Resample aspect map
# compute focalmean
???

# compute weights 225, 270, 315 and 360
# this is the basic formula how in r,mapcalc?
r.mapcalc
w225 = sin2(aspect - 225°)
w270 = sin2(aspect - 270°)
w315 = sin2(aspect - 315°)
w360 = sin2(aspect - 360°)

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

Any ideas?

thanks,

maning
--
|---------|----------------------------------------------------------|
| __.-._ |"Ohhh. Great warrior. Wars not make one great." -Yoda |
| '-._"7' |"Freedom is still the most radical idea of all" -N.Branden|
| /'.-c |Linux registered user #402901, http://counter.li.org/ |
| | /T |www.esambale.wikispaces.com|
| _)_/LI |www.geocities.com/esambale/philbiodivmap/philbirds.html |
|---------|----------------------------------------------------------|

2

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

On 6/4/06, maning sambale <emmanuel.sambale@gmail.com> wrote:

Hello,

I found an article creating a multidirectional, oblique-weighted,
shaded-relief image from the USGS
(http://pubs.usgs.gov/of/1992/of92-422/of92-422.pdf).

I basically creates a hillshade using different illumination angles.
This is especially useful in very mountainous areas. I tried to
recreate this with GRASS but having difficulty with r.mapcalc
equations.

Heres what I've made so far:

# Multidirectional, oblique-weighted, shaded-relief image

# Compute hillshade at azimuth 225, 270, 315 and 360 at 30 degrees sun
illumination angle
r.shaded.relief map=DEM_90@PERMANENT shadedmap=shade225 altitude=30
azimuth=225 zmult=1 scale=1
r.shaded.relief map=DEM_90@PERMANENT shadedmap=shade270 altitude=30
azimuth=270 zmult=1 scale=1
r.shaded.relief map=DEM_90@PERMANENT shadedmap=shade315 altitude=30
azimuth=315 zmult=1 scale=1
r.shaded.relief map=DEM_90@PERMANENT shadedmap=shade360 altitude=30
azimuth=360 zmult=1 scale=1

# compute aspect map
r.slope.aspect elevation=DEM_90@PERMANENT format=degrees prec=float
aspect=aspect zfactor=1.0 min_slp_allowed=0.0

# Resample aspect map
# compute focalmean
???

# compute weights 225, 270, 315 and 360
# this is the basic formula how in r,mapcalc?
r.mapcalc
w225 = sin2(aspect - 225°)
w270 = sin2(aspect - 270°)
w315 = sin2(aspect - 315°)
w360 = sin2(aspect - 360°)

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

Any ideas?

thanks,

maning
--
|---------|----------------------------------------------------------|
| __.-._ |"Ohhh. Great warrior. Wars not make one great." -Yoda |
| '-._"7' |"Freedom is still the most radical idea of all" -N.Branden|
| /'.-c |Linux registered user #402901, http://counter.li.org/ |
| | /T |www.esambale.wikispaces.com|
| _)_/LI |www.geocities.com/esambale/philbiodivmap/philbirds.html |
|---------|----------------------------------------------------------|

_______________________________________________
grassuser mailing list
grassuser@grass.itc.it
http://grass.itc.it/mailman/listinfo/grassuser

--
David Finlayson

3

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@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:
# 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