[GRASS-SVN] r33621 - grass/trunk/scripts/r.shaded.relief
svn_grass at osgeo.org
svn_grass at osgeo.org
Tue Sep 30 12:58:06 EDT 2008
Author: glynn
Date: 2008-09-30 12:58:06 -0400 (Tue, 30 Sep 2008)
New Revision: 33621
Added:
grass/trunk/scripts/r.shaded.relief/r.shaded.relief.py
Log:
Convert r.shaded.relief to Python
Added: grass/trunk/scripts/r.shaded.relief/r.shaded.relief.py
===================================================================
--- grass/trunk/scripts/r.shaded.relief/r.shaded.relief.py (rev 0)
+++ grass/trunk/scripts/r.shaded.relief/r.shaded.relief.py 2008-09-30 16:58:06 UTC (rev 33621)
@@ -0,0 +1,195 @@
+#!/usr/bin/env python
+
+############################################################################
+#
+# MODULE: r.shaded.relief
+# AUTHOR(S): CERL
+# parameters standardized: Markus Neteler, 2008
+# updates: Michael Barton, 2004
+# updates: Gordon Keith, 2003
+# updates: Andreas Lange, 2001
+# updates: David Finlayson, 2001
+# updates: Markus Neteler, 2001, 1999
+# Converted to Python by Glynn Clements
+# PURPOSE: Creates shaded relief map from raster elevation map (DEM)
+# COPYRIGHT: (C) 1999 - 2008 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.
+#
+#############################################################################
+#
+# July 2007 - allow input from other mapsets (Brad Douglas)
+#
+# May 2005 - fixed wrong units parameter (Markus Neteler)
+#
+# September 2004 - Added z exaggeration control (Michael Barton)
+# April 2004 - updated for GRASS 5.7 by Michael Barton
+#
+# 9/2004 Adds scale factor input (as per documentation); units set scale only if specified for lat/long regions
+# Also, adds option of controlling z-exaggeration.
+#
+# 6/2003 fixes for Lat/Long Gordon Keith <gordon.keith at csiro.au>
+# If n is a number then the ewres and nsres are mulitplied by that scale
+# to calculate the shading.
+# If n is the letter M (either case) the number of metres is degree of
+# latitude is used as the scale.
+# If n is the letter f then the number of feet in a degree is used.
+# It scales latitude and longitude equally, so it's only approximately
+# right, but for shading its close enough. It makes the difference
+# between an unusable and usable shade.
+#
+# 10/2001 fix for testing for dashes in raster file name
+# by Andreas Lange <andreas.lange at rhein-main.de>
+# 10/2001 added parser support - Markus Neteler
+# 9/2001 fix to keep NULLs as is (was value 22 before) - Markus Neteler
+# 1/2001 fix for NULL by David Finlayson <david_finlayson at yahoo.com>
+# 11/99 updated $ewres to ewres() and $nsres to nsres()
+# updated number to FP in r.mapcalc statement Markus Neteler
+
+#% Module
+#% description: Creates shaded relief map from an elevation map (DEM).
+#% keywords: raster, elevation
+#% End
+#% option
+#% key: input
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: Input elevation map
+#% required : yes
+#% end
+#% option
+#% key: output
+#% type: string
+#% gisprompt: new,cell,raster
+#% description: Output shaded relief map name
+#% required : no
+#% end
+#% option
+#% key: altitude
+#% type: integer
+#% description: Altitude of the sun in degrees above the horizon
+#% required : no
+#% options : 0-90
+#% answer: 30
+#% end
+#% option
+#% key: azimuth
+#% type: integer
+#% description: Azimuth of the sun in degrees to the east of north
+#% required : no
+#% options : 0-360
+#% answer: 270
+#% end
+#% option
+#% key: zmult
+#% type: double
+#% description: Factor for exaggerating relief
+#% required : no
+#% answer: 1
+#% end
+#% option
+#% key: scale
+#% type: double
+#% description: Scale factor for converting horizontal units to elevation units
+#% required : no
+#% answer: 1
+#% end
+#% option
+#% key: units
+#% type: string
+#% description: Set scaling factor (applies to lat./long. locations only, none: scale=1)
+#% required : no
+#% options: none,meters,feet
+#% answer: none
+#% end
+
+import sys
+import os
+import string
+import grass
+
+def main():
+ input = options['input']
+ output = options['output']
+ altitude = options['altitude']
+ azimuth = options['azimuth']
+ zmult = options['zmult']
+ scale = options['scale']
+
+ if not grass.findfile(elev)['file']:
+ grass.fatal("Map <%s> not found." % elev)
+
+ if input == output:
+ grass.fatal("Input elevation map and output relief map must have different names")
+
+ elev = input
+
+ if not output:
+ elevshade = elev.split('@')[0] + '.shade'
+ else:
+ elevshade = output
+
+ elev_out = "'%s'" % elevshade
+
+ alt = altitude
+ scale = float(scale)
+
+ # LatLong locations only:
+ if units == 'meters':
+ #scale=111120
+ scale *= 1852 * 60
+
+ # LatLong locations only:
+ if units == 'feet':
+ #scale=364567.2
+ scale *= 6076.12 * 60
+
+ #correct azimuth to East (GRASS convention):
+ az = azimuth - 90
+
+ grass.message("Calculating shading, please stand by.")
+
+ t = string.Template(r'$elev_out = eval(
+x=($zmult*$elev[-1,-1] + 2*$zmult*$elev[0,-1] + $zmult*$elev[1,-1] - \
+ $zmult*$elev[-1,1] - 2*$zmult*$elev[0,1] - $zmult*$elev[1,1])/(8.*ewres()*$scale),
+y=($zmult*$elev[-1,-1] + 2*$zmult*$elev[-1,0] + $zmult*$elev[-1,1] - \
+ $zmult*$elev[1,-1] - 2*$zmult*$elev[1,0] - $zmult*$elev[1,1])/(8.*nsres()*$scale),
+slope=90.-atan(sqrt(x*x + y*y)),
+a=round(atan(x,y)),
+a=if(isnull(a),1,a),
+aspect=if(x!=0||y!=0,if(a,a,360.)),
+cang = sin($alt)*sin(slope) + cos($alt)*cos(slope) * cos($az-aspect),
+if(cang < 0.,0.,100.*cang),
+if(isnull(cang), null(), 100.*cang))')
+ expr = t.substitute(alt = alt, az = az, elev = elev, elev_out = elev_out, scale = scale, zmult = zmult)
+ p = grass.feed_command('r.mapcalc')
+ p.stdin.write(expr)
+ p.stdin.close()
+
+ if p.wait() != 0:
+ grass.fatal("In calculation, script aborted.")
+
+ grass.run_command('r.colors', map = elevshade, color = 'grey')
+
+ # record metadata
+ grass.run_command('r.support', map = elevshade, title = 'Shaded relief of "%s"' % input, history = '')
+ grass.run_command('r.support', map = elevshade, history = "r.shaded.relief settings:")
+ t = string.Template("altitude=$alt azimuth=$azimuth zmult=$zmult scale=$scale")
+ parms = dict(alt = alt, azimuth = azimuth, zmult = zmult, scale = options['scale'])
+ grass.run_command('r.support', map = elevshade, history = t.substitute(parms))
+ if units:
+ grass.run_command('r.support', map = elevshade, history = " units=%s (adjusted scale=%s)" % (units, scale))
+
+ if output:
+ grass.message("Shaded relief map created and named <%s>." % elevshade)
+ else:
+ grass.message("Shaded relief map created and named <%s>. Consider renaming." % elevshade)
+
+ # write cmd history:
+ raster_history(elevshade)
+
+if __name__ == "__main__":
+ options, flags = grass.parser()
+ main()
Property changes on: grass/trunk/scripts/r.shaded.relief/r.shaded.relief.py
___________________________________________________________________
Name: svn:executable
+ *
More information about the grass-commit
mailing list