[GRASS-SVN] r60555 - grass-addons/grass7/raster/r.roughness.vector
svn_grass at osgeo.org
svn_grass at osgeo.org
Wed May 28 11:22:53 PDT 2014
Author: hellik
Date: 2014-05-28 11:22:53 -0700 (Wed, 28 May 2014)
New Revision: 60555
Added:
grass-addons/grass7/raster/r.roughness.vector/r.roughness.vector.html
grass-addons/grass7/raster/r.roughness.vector/r.roughness.vector.py
Log:
r.roughness.vector: add pythonized grass6 addon (on behalf of Carlos Grohmann) - step 2
Added: grass-addons/grass7/raster/r.roughness.vector/r.roughness.vector.html
===================================================================
--- grass-addons/grass7/raster/r.roughness.vector/r.roughness.vector.html (rev 0)
+++ grass-addons/grass7/raster/r.roughness.vector/r.roughness.vector.html 2014-05-28 18:22:53 UTC (rev 60555)
@@ -0,0 +1,55 @@
+<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
+<html><head><title>r.roughness.window.vector</title>
+<meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
+<link rel="stylesheet" href="grassdocs.css" type="text/css"></head>
+
+<body style="background-color: white;">
+<img style="width: 76px; height: 91px;" src="grass_logo.png" alt="GRASS logo"><hr align="center" noshade="noshade" size="6"><h2>NAME</h2>
+<em><b>r.roughness.vector.py</b></em> - Calculates surface roughness in a moving-window, as the orientation of vectors normal to surface planes.
+<h2>KEYWORDS </h2>
+raster, elevation, slope, aspect
+<h2>SYNOPSIS</h2>
+<span style="font-weight: bold;">r.roughness.window.vector</span><br style="font-weight: bold;">
+<span style="font-weight: bold;">r.roughness.window.vector help</span><br>
+<span style="font-weight: bold;">r.roughness.window.vector</span> <b>map</b>=<em>string</em>
+[<b>slope</b>=<em>string</em>] [<b>aspect</b>=<em>stri</em>ng]
+[<b>window</b>=<em>integer</em>] [<span style="font-weight: bold;">strength</span>=<em>string</em>]
+[<span style="font-weight: bold;">fisher</span>=<em>string</em>] [<b>compass</b>=<em>string</em>]
+[<b>colatitude</b>=<em>string</em>] [<span style="font-weight: bold;">xcos</span><b></b>=<em>string</em>] [<span style="font-weight: bold;">ycos</span>=<em>string</em>] [<span style="font-weight: bold;">zcos</span>=<em>string</em>] [--<b>overwrite</b>] [--<b>verbose</b>] [--<b>quiet</b>]
+<h3>Flags:</h3>
+<dl><dt><b>--overwrite</b></dt>
+<dd>Allow output files to overwrite existing files</dd>
+<dt><b>--verbose</b></dt>
+<dd>Verbose module output</dd>
+<dt><b>--quiet</b></dt>
+<dd>Quiet module output</dd>
+</dl><h3>Parameters:</h3>
+<dl><dt><b>map</b>=<em>string</em></dt>
+<dd>Input elevation raster map</dd>
+<dt><b>slope</b>=<em>string</em></dt>
+<dd>Input slope map</dd>
+<dt><b>aspect</b>=<em>integer</em></dt>
+<dd>Input aspect map</dd><dt><b>window</b>=<em>integer</em></dt>
+<dd>Moving-window size (uses <span style="font-style: italic;">r.neighbors</span>)</dd><dt><b>strength</b>=<em>string (optional)</em></dt>
+<dd>Output "vector strength" map </dd><dt><b>fisher</b>=<em>string (optional)</em><em></em></dt>
+<dd>Output "Fisher's K parameter" map </dd><dt><b>compass</b>=<em>string</em><em> (optional)</em></dt>
+<dd>Input compass aspect map</dd><dt><b>colatitude</b>=<em>string</em><em> (optional)</em></dt><dd>Input colatitude map</dd><dt><b>xcos</b>=<em>string</em><em> (optional)</em></dt><dd>Input <span style="font-style: italic;">x</span> directional cosine map </dd><dt><b>ycos</b>=<em>string</em><em> (optional)</em></dt><dd>Input <span style="font-style: italic;">y</span> directional cosine map </dd><dt><b>zcos</b>=<em>string</em><em> (optional)</em></dt><dd>Input <span style="font-style: italic;">z</span> directional cosine map </dd></dl><h2>DESCRIPTION</h2>
+In this script surface roughness is taken as the dispersion of vectors
+normal to surface areas (pixels). Normal vectors are defined by
+slope and aspect.<br><br>This script will create several temporary
+maps, for the directional cosines in each direction (x,y,z), for the
+sum of these cosines and vector strengh.<br><br>The options <span style="font-style: italic;">compass</span>, <span style="font-style: italic;">colatitude</span>, <span style="font-style: italic;">xcos</span>, <span style="font-style: italic;">ycosm</span> and <span style="font-style: italic;">zcos</span>
+are created as temporary files each time the script is run. If the user
+wants to create several map (with different window sizes, for
+instance), it is recommended to create those maps with <span style="font-style: italic;">r.mapcalc</span> and use them as input:<br><br><span style="font-style: italic;">r.mapcalc compass = "if(aspect==0,0,if(aspect < 90, 90-aspect, 360+90-aspect))"</span><br style="font-style: italic;"><br style="font-style: italic;"><span style="font-style: italic;">r.mapcalc colatitude = "90 - slope"</span><br style="font-style: italic;"><br style="font-style: italic;"><span style="font-style: italic;">r.mapcalc xcos = "sin(colatitude)*cos(compass)"</span><br style="font-style: italic;"><br style="font-style: italic;"><span style="font-style: italic;">r.mapcalc ycos = "sin(colatitude)*sin(compass)"</span><br style="font-style: italic;"><br style="font-style: italic;"><span style="font-style: italic;">r.mapcalc zcos = "cos(colatitude)"</span><br style="font-style: italic;"><br><br><br>If the user does not specify the output maps names, they will be set to <span style="font-style: i
talic;"><br><br>INPUT_MAP_vector_strength_NxN</span> <br>and<br><span style="font-style: italic;">INPUT_MAP_fisher_K_NxN</span><br><br>where N is the window size.<p>
+</p><h2>REFERENCES</h2><p>Hobson, R.D., 1972. Surface roughness in topography: quantitative approach. In: Chorley, R.J. (ed) <span style="font-style: italic;">Spatial analysis in geomorphology</span>. Methuer, London, p.225-245.<br></p><p>McKean, J. & Roering, J., 2004. <span style="font-style: italic;">Objective landslide detection and surface morphology mapping using high-resolution airborne laser altimetry. Geomorphology</span>, 57:331-351</p>
+<h2>AUTHORS</h2>
+Carlos Henrique Grohmann - Institute of Energy and Environment, University of São Paulo, Brazil.<br>
+Helmut Kudrnovsky - http://www.alectoria.at
+<p>
+
+<i>Last changed: $Date: 2006/04/20 02:43:23 $</i>
+</p><hr><p><a href="index.html">Main
+index</a> - <a href="raster.html">raster index</a>
+- <a href="full_index.html">Full index</a></p>
+</body></html>
\ No newline at end of file
Added: grass-addons/grass7/raster/r.roughness.vector/r.roughness.vector.py
===================================================================
--- grass-addons/grass7/raster/r.roughness.vector/r.roughness.vector.py (rev 0)
+++ grass-addons/grass7/raster/r.roughness.vector/r.roughness.vector.py 2014-05-28 18:22:53 UTC (rev 60555)
@@ -0,0 +1,352 @@
+#!/usr/bin/env python
+#-*- coding:utf-8 -*-
+#
+#
+############################################################################
+#
+# MODULE: r.roughness.vector.py
+# AUTHOR(S): Carlos H. Grohmann <carlos dot grohmann at gmail dot com >
+# Helmut Kudrnovsky <alectoria at gmx dot at>
+#
+# PURPOSE: Calculates surface roughness from DEMs.
+# Python version of r.roughness.vector.sh
+#
+# In this script surface roughness is taken as the dispersion
+# of vectors normal to surface areas (pixels). Normal vectors
+# are defined by slope and aspect.
+# Reference:
+# Hobson, R.D., 1972. Surface roughness in topography:
+# quantitative approach. In: Chorley, R.J. (ed) Spatial
+# analysis in geomorphology. Methuer, London, p.225-245.
+#
+# This script will create several temporary maps, for the
+# directional cosines in each direction (x,y,z), for the sum
+# of these cosines and vector strengh.
+#
+# If the user does not specify the output map name, it will be
+# set to INPUT_MAP_roughness_vector_NxN
+# where N is the window size.
+#
+# COPYRIGHT: (C) 2014 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: Calculates surface roughness with a moving-window approach
+#% keywords: raster
+#% keywords: terrain
+#% keywords: aspect
+#% keywords: slope
+#% keywords: roughness
+#%end
+
+#%option G_OPT_R_ELEV
+#% key: elevation
+#% description: Name of elevation raster map
+#% required: yes
+#%end
+
+#%option
+#% key: slope
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: Input raster slope map
+#% required : yes
+#%end
+#%option
+#% key: aspect
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: Input raster aspect map
+#% required : yes
+#%end
+#%option
+#% key: window
+#% type: integer
+#% description: Window size (3,5,7,...,25)
+#% required : no
+#% answer : 3
+#%end
+#%option
+#% key: strength
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: Output vector strength map
+#% required : no
+#%end
+#%option
+#% key: fisher
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: Output Fischer's K parameter map
+#% required : no
+#%end
+#%option
+#% key: compass
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: OPTIONAL INPUT: Compass aspect values (longitude)
+#% required : no
+#%end
+#%option
+#% key: colatitude
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: OPTIONAL INPUT: Colatitude values (90 - slope)
+#% required : no
+#%end
+#%option
+#% key: xcos
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: OPTIONAL INPUT: X directional cosine map
+#% required : no
+#%end
+#%option
+#% key: ycos
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: OPTIONAL INPUT: Y directional cosine map
+#% required : no
+#%end
+#%option
+#% key: zcos
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: OPTIONAL INPUT: Z directional cosine map
+#% required : no
+#%end
+#
+
+import sys
+import atexit
+import grass.script as grass
+
+# cleaning up tem files
+def cleanup():
+ rasts = ['aspect_compass','colat_angle','cosine_x','cosine_y','cosine_z','sum_Xcosine','sum_Ycosine','sum_Zcosine']
+ grass.run_command('g.remove', flags = 'f', rast = rasts, quiet = True)
+
+def main():
+
+#set options
+ elevmap = options['elevation'] # input
+ slope = options['slope'] # input
+ aspect = options['aspect'] # input
+ window = options['window'] # input
+ strength = options['strength'] # output
+ fisher = options['fisher'] # output
+ compass = options['compass'] # temporary
+ colatitude = options['colatitude'] # temporary
+ xcos = options['xcos'] # temporary
+ ycos = options['ycos'] # temporary
+ zcos = options['zcos'] # temporary
+
+
+# check if input files exist
+ grass.message( "----" )
+ grass.message( "Check if input files exist ..." )
+
+ find_elev = grass.find_file(elevmap, element = 'cell')
+ if find_elev['name'] == "":
+ print "Map %s not found! Aborting." % elevmap
+ sys.exit()
+
+ find_slope = grass.find_file(slope, element = 'cell')
+ if find_slope['name'] == "":
+ print "Map %s not found! Aborting." % slope
+ sys.exit()
+
+ find_aspect = grass.find_file(aspect, element = 'cell')
+ if find_aspect['name'] == "":
+ print "Map %s not found! Aborting." % aspect
+ sys.exit()
+
+#########################################################################################################
+
+#ACERTAR O NOME DOS ARQUIVOS - TIRAR TUDO DESDE O @!!!!
+
+#########################################################################################################
+
+
+
+# give default names to outputs, in case the user doesn't provide them
+ grass.message( "----" )
+ grass.message( "Define default output names when not defined by user ..." )
+
+ if strength == "":
+ strength = "%s_vector_strength_%sx%s" % (find_elev['name'],window,window)
+
+
+ if fisher == "":
+ fisher = "%s_fisher_1K_%sx%s" % (find_elev['name'],window,window)
+
+#####################
+# calculate compass-oriented aspect and colatitude
+# (temp rasters)
+
+# correct aspect angles from cartesian (GRASS default) to compass angles
+# if(A==0,0,if(A < 90, 90-A, 360+90-A))
+
+ grass.message( "----" )
+ grass.message( "Calculate compass aspect values ..." )
+
+ if compass == "":
+ aspect_compass = 'aspect_compass'
+# print "Calculating compass aspect values (longitude)"
+# aspect_compass = grass.tempfile()
+ grass.mapcalc("${out} = if(${rast1}==0,0,if(${rast1} < 90, 90-${rast1}, 360+90-${rast1}))",
+ out = aspect_compass,
+ rast1 = aspect)
+ else:
+ grass.message( "Using previous calculated compass aspect values (longitude)" )
+ # print "Using previous calculated compass aspect values (longitude)"
+ aspect_compass = compass
+
+# calculates colatitude (90-slope)
+
+ grass.message( "----" )
+ grass.message( "Calculate colatitude ..." )
+
+ if colatitude == "":
+ colat_angle = 'colat_angle'
+# print "Calculating colatitude (90-slope)"
+# colat_angle = grass.tempfile()
+ grass.mapcalc("${out} = 90 - ${rast1}",
+ out = colat_angle,
+ rast1 = slope)
+ else:
+ grass.message( "Using previous calculated colatitude values" )
+ # print "Using previous calculated colatitude values"
+ colat_angle = colatitude
+
+#####################
+# calculate direction cosines
+# direction cosines relative to axis oriented north, east and up
+# direction cosine calculation according to McKean & Roering (2004), Geomorphology, 57:331-351.
+
+ grass.message( "----" )
+ grass.message( "Calculate direction cosines ..." )
+
+# X cosine
+ if xcos == "":
+ cosine_x = 'cosine_x'
+# print "Calculating X direction cosine"
+# cosine_x = grass.tempfile()
+ grass.mapcalc("${out} = sin(${rast1}) * cos(${rast2})",
+ out = 'cosine_x',
+ rast1 = aspect_compass,
+ rast2 = colat_angle)
+ else:
+ grass.message( "Using previous calculated X direction cosine value" )
+ # print "Using previous calculated X direction cosine values"
+ cosine_x = xcos
+
+# Y cosine
+ if ycos == "":
+ cosine_y = 'cosine_y'
+# print "Calculating Y direction cosine"
+# cosine_y = grass.tempfile()
+ grass.mapcalc("${out} = sin(${rast1}) * sin(${rast2})",
+ out = 'cosine_y',
+ rast1 = aspect_compass,
+ rast2 = colat_angle)
+ else:
+ grass.message( "Using previous calculated Y direction cosine values" )
+ # print "Using previous calculated Y direction cosine values"
+ cosine_y = ycos
+
+# Z cosine
+ if zcos == "":
+ cosine_z = 'cosine_z'
+# print "Calculating Y direction cosine"
+# cosine_z = grass.tempfile()
+ grass.mapcalc("${out} = cos(${rast1})",
+ out = 'cosine_z',
+ rast1 = aspect_compass)
+ else:
+ grass.message( "Using previous calculated Y direction cosine values" )
+ # print "Using previous calculated Y direction cosine values"
+ cosine_z = zcos
+
+
+# calculate SUM of direction cosines
+
+ grass.message( "----" )
+ grass.message( "Calculate sum of direction cosines ..." )
+
+ grass.message( "Calculating sum of X direction cosines ..." )
+ # print "Calculating sum of X direction cosines"
+# sum_Xcosine = grass.tempfile()
+ grass.run_command("r.neighbors",
+ input=cosine_x,
+ output='sum_Xcosine',
+ method='sum',
+ size=window,
+ overwrite=True)
+
+ grass.message( "Calculating sum of Y direction cosines ..." )
+# print "Calculating sum of Y direction cosines"
+# sum_Ycosine = grass.tempfile()
+ grass.run_command("r.neighbors",
+ input=cosine_y,
+ output='sum_Ycosine',
+ method='sum',
+ size=window,
+ overwrite=True)
+
+ grass.message( "Calculating sum of Z direction cosines ..." )
+# print "Calculating sum of Z direction cosines"
+# sum_Zcosine = grass.tempfile()
+ grass.run_command("r.neighbors",
+ input=cosine_z,
+ output='sum_Zcosine',
+ method='sum',
+ size=window,
+ overwrite=True)
+
+#####################
+# calculate vector strength
+
+ grass.message( "----" )
+ grass.message( "Calculate vector strength ..." )
+
+# print "Calculating vector strength"
+# print strength
+ grass.mapcalc("${out} = sqrt(exp(${rast1},2) + exp(${rast2},2) + exp(${rast3},2))",
+ out = strength,
+ rast1 = 'sum_Xcosine',
+ rast2 = 'sum_Ycosine',
+ rast3 = 'sum_Zcosine')
+
+# calculate Inverted Fisher's K parameter
+# k=1/((N-1)/(N-R))
+
+ grass.message( "----" )
+ grass.message( "Calculate inverted Fisher's K parameter ..." )
+
+ w = int(window)
+# print "Calculating Inverted Fisher's K parameter"
+ grass.mapcalc("${out} = ($w * $w - ${rast1}) / ($w * $w - 1)",
+ out = fisher,
+ rast1 = strength,
+ w = int(window))
+
+# calculations done
+
+ grass.message( "----" )
+ grass.message( "Result maps:" )
+ grass.message( strength )
+ grass.message( fisher )
+ grass.message( "Calculations done." )
+ grass.message( "----" )
+
+# this "if" condition instructs execution of code contained in this script, *only* if the script is being executed directly
+if __name__ == "__main__": # this allows the script to be used as a module in other scripts or as a standalone script
+ options, flags = grass.parser() #
+ atexit.register(cleanup)
+ sys.exit(main()) #
More information about the grass-commit
mailing list