[GRASS-SVN] r41381 - grass/trunk/scripts/r.shaded.relief
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu Mar 11 14:35:39 EST 2010
Author: martinl
Date: 2010-03-11 14:35:37 -0500 (Thu, 11 Mar 2010)
New Revision: 41381
Modified:
grass/trunk/scripts/r.shaded.relief/r.shaded.relief.py
Log:
r.shaded.relief: minor clean up
Modified: grass/trunk/scripts/r.shaded.relief/r.shaded.relief.py
===================================================================
--- grass/trunk/scripts/r.shaded.relief/r.shaded.relief.py 2010-03-11 13:33:29 UTC (rev 41380)
+++ grass/trunk/scripts/r.shaded.relief/r.shaded.relief.py 2010-03-11 19:35:37 UTC (rev 41381)
@@ -12,7 +12,7 @@
# 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
+# COPYRIGHT: (C) 1999 - 2008, 2010 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
@@ -51,21 +51,20 @@
#% Module
#% description: Creates shaded relief map from an elevation map (DEM).
#% keywords: raster
-#% keywords: elevation
#% End
#% option
#% key: input
#% type: string
#% gisprompt: old,cell,raster
-#% description: Input elevation map
+#% description: Name of input elevation map
#% required : yes
#% end
#% option
#% key: output
#% type: string
#% gisprompt: new,cell,raster
-#% description: Output shaded relief map name
-#% required : no
+#% description: Name for output shaded relief map
+#% required : yes
#% end
#% option
#% key: altitude
@@ -96,6 +95,7 @@
#% description: Scale factor for converting horizontal units to elevation units
#% required : no
#% answer: 1
+#% guisection: Scaling
#% end
#% option
#% key: units
@@ -104,6 +104,7 @@
#% required : no
#% options: none,meters,feet
#% answer: none
+#% guisection: Scaling
#% end
import sys
@@ -117,57 +118,46 @@
altitude = options['altitude']
azimuth = options['azimuth']
zmult = options['zmult']
- scale = options['scale']
+ scale = float(options['scale'])
units = options['units']
-
+ verbose_level = os.getenv('GRASS_VERBOSE')
+ if verbose_level is None:
+ verbose_level = 2
+
if not grass.find_file(input)['file']:
- grass.fatal(_("Map <%s> not found.") % input)
+ grass.fatal(_("Raster map <%s> not found") % input)
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=111120
scale *= 1852 * 60
# LatLong locations only:
if units == 'feet':
- #scale=364567.2
+ # scale=364567.2
scale *= 6076.12 * 60
#correct azimuth to East (GRASS convention):
# this seems to be backwards, but in fact it works so leave it.
az = float(azimuth) - 90
- grass.message(_("Calculating shading, please stand by."))
-
t = string.Template(
r'''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), \
+ x=($zmult*$input[-1,-1] + 2*$zmult*$input[0,-1] + $zmult*$input[1,-1] - \
+ $zmult*$input[-1,1] - 2*$zmult*$input[0,1] - $zmult*$input[1,1])/(8.*ewres()*$scale), \
+ y=($zmult*$input[-1,-1] + 2*$zmult*$input[-1,0] + $zmult*$input[-1,1] - \
+ $zmult*$input[1,-1] - 2*$zmult*$input[1,0] - $zmult*$input[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) \
+ cang = sin($altitude)*sin(slope) + cos($altitude)*cos(slope) * cos($az-aspect) \
)
- $elev_out = if(isnull(cang), null(), 100.*cang)''')
- expr = t.substitute(alt = alt, az = az, elev = elev, elev_out = elev_out, scale = scale, zmult = zmult)
+ $output = if(isnull(cang), null(), 100.*cang)''')
+ expr = t.substitute(altitude = altitude, az = az, input = input, output = output, scale = scale, zmult = zmult)
p = grass.feed_command('r.mapcalc')
p.stdin.write(expr)
p.stdin.close()
@@ -175,24 +165,22 @@
if p.wait() != 0:
grass.fatal(_("In calculation, script aborted."))
- grass.run_command('r.colors', map = elevshade, color = 'grey')
-
+ if verbose_level < 3:
+ grass.run_command('r.colors', map = output, color = 'grey', quiet = True)
+ else:
+ grass.run_command('r.colors', map = output, 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))
+ grass.run_command('r.support', map = output, title = 'Shaded relief of "%s"' % input, history = '')
+ grass.run_command('r.support', map = output, history = "r.shaded.relief settings:")
+ t = string.Template("altitude=$altitude azimuth=$azimuth zmult=$zmult scale=$scale")
+ parms = dict(altitude = altitude, azimuth = azimuth, zmult = zmult, scale = options['scale'])
+ grass.run_command('r.support', map = output, history = t.substitute(parms))
if units:
- grass.run_command('r.support', map = elevshade, history = " units=%s (adjusted scale=%s)" % (units, scale))
+ grass.run_command('r.support', map = output, 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:
- grass.raster_history(elevshade)
+ grass.raster_history(output)
if __name__ == "__main__":
options, flags = grass.parser()
More information about the grass-commit
mailing list