-
-
Save amcox/b54d65f08fe4448840462901fd6350d4 to your computer and use it in GitHub Desktop.
Creates a multidirectional, oblique-weighted, shaded-relief image from an input DEM image.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| #!/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 | |
| #%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 | |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment