[GRASS-SVN] r38695 - in grass-addons/raster: . r.roughness

svn_grass at osgeo.org svn_grass at osgeo.org
Tue Aug 11 11:05:24 EDT 2009


Author: guano
Date: 2009-08-11 11:05:22 -0400 (Tue, 11 Aug 2009)
New Revision: 38695

Added:
   grass-addons/raster/r.roughness/
   grass-addons/raster/r.roughness/r.roughness.sh
   grass-addons/raster/r.roughness/r.roughness.window.area
   grass-addons/raster/r.roughness/r.roughness.window.vector
   grass-addons/raster/r.roughness/r.roughness.window.vector.html
Log:
Surface roughness modules (scripts)

Added: grass-addons/raster/r.roughness/r.roughness.sh
===================================================================
--- grass-addons/raster/r.roughness/r.roughness.sh	                        (rev 0)
+++ grass-addons/raster/r.roughness/r.roughness.sh	2009-08-11 15:05:22 UTC (rev 38695)
@@ -0,0 +1,247 @@
+#!/bin/sh
+#
+############################################################################
+#
+# MODULE:	r.roughness
+# AUTHOR(S):	Carlos H. Grohmann <carlos dot grohmann at gmail dot com >
+# PURPOSE:	Calculates surface roughness from DEMs. (uses r.surf.area)
+#		In this script surface roughness is used in the sense of 
+#		Hobson (1972), who describes it as the ratio between surface 
+#		(real) area and flat (plan) area of square cells; in this 
+#		approach, flat surfaces would present values close to 1, 
+#		whilst in irregular ones the ratio shows a curvilinear 
+#		relationship which asymptotically approaches infinity as the 
+#		real areas increases.
+#		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 square sub-regions with size defined by
+#		the option GRID. In each sub-region, the real and planar areas
+#		will be calculated by r.surf.area, and the results (points at 
+#		the center of sub-regions) will be interpolated with v.surf.rst.
+#		The user also can set the tension and smooth parameters.
+#
+# COPYRIGHT:	(C) 2006-2009 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
+#%End
+#%option
+#% key: map
+#% gisprompt: old,cell,raster
+#% type: string
+#% description: Input raster surface
+#% required : yes
+#%end
+#%option
+#% key: grid
+#% type: integer
+#% description: Grid size in meters
+#% required : no
+#% answer : 1000
+#%end
+#%option
+#% key: rough
+#% gisprompt: old,cell,raster
+#% type: string
+#% description: Output raster map 
+#% required : no
+#%end
+#%option
+#% key: tension
+#% type: double
+#% description: Spline tension parameter
+#% required : no
+#% answer : 40.
+#%end
+#%option
+#% key: smooth
+#% type: double
+#% description: Spline smoothing parameter
+#% required : no
+#% answer : 0.1
+#%end
+
+
+if test "$GISBASE" = ""; then
+ echo "You must be in GRASS GIS to run this program." >&2
+ exit 1
+fi   
+
+# save command line
+if [ "$1" != "@ARGS_PARSED@" ] ; then
+    CMDLINE="`basename $0`"
+    for arg in "$@" ; do
+        CMDLINE="$CMDLINE \"$arg\""
+    done
+    export CMDLINE
+    exec g.parser "$0" "$@"
+fi
+
+PROG=`basename $0`
+
+
+#### check if we have awk
+if [ ! -x "`which awk`" ] ; then
+    echo "$PROG: awk required, please install awk or gawk first" 2>&1
+    exit 1
+fi
+
+# setting environment, so that awk works properly in all languages
+unset LC_ALL
+export LC_NUMERIC=C
+
+eval `g.gisenv`
+: ${GISBASE?} ${GISDBASE?} ${LOCATION_NAME?} ${MAPSET?}
+LOCATION=$GISDBASE/$LOCATION_NAME/$MAPSET
+
+TMP_ascii="`g.tempfile pid=$$`"
+if [ $? -ne 0 ] || [ -z "${TMP_ascii}" ] ; then
+    echo "ERROR: unable to create temporary files" 1>&2
+    exit 1
+fi
+
+
+#vars
+
+elev=$GIS_OPT_MAP
+grid=$GIS_OPT_GRID
+#rough=$GIS_OPT_ROUGH
+
+#check if input file exists
+eval `g.findfile element=cell file=$elev`
+if [ -z "$name" ] ; then
+   echo "ERROR: map <$elev> not found."
+   exit 1
+fi
+
+if [ "$GIS_OPT_MAP" = "$GIS_OPT_ROUGH" ]; then
+	echo ""
+	echo "Input elevation map and output roughness map must have different names"
+	exit 1
+fi
+
+if [ -z "$GIS_OPT_ROUGH" ]; then
+    ROUGHNESS="${elev}_roughness_${grid}"
+else
+    ROUGHNESS="$GIS_OPT_ROUGH"
+fi
+
+
+#######################################################################
+cleanup() 
+{
+    g.region region=oldregion 
+    rm -f $GISDBASE/$LOCATION_NAME/$MAPSET/windows/oldregion
+    rm -f $TMP ${TMP} TMP_*
+    g.remove vect=TMP_vect > /dev/null
+}
+
+# what to do in case of user break:
+exitprocedure()
+{
+    echo "User break!"
+    cleanup
+    exit 1
+}
+# shell check for user break (signal list: trap -l)
+trap "exitprocedure" 2 3 15
+
+
+#######################################################################
+
+########################################################################
+# get region limits
+maxnorth="`g.region -p | grep north | sed -e s/.*:\ *//`"
+maxsouth="`g.region -p | grep south | sed -e s/.*:\ *//`"
+maxwest="` g.region -p | grep west | sed -e s/.*:\ *//`"
+maxeast="` g.region -p | grep east | sed -e s/.*:\ *//`"
+
+g.region save=oldregion
+
+
+ns_dist="`echo $maxnorth $maxsouth |awk '{printf("%f", $1 - $2);}'`"
+ew_dist="`echo $maxeast $maxwest |awk '{printf("%f", $1 - $2);}'`"
+rows="`echo $ns_dist $grid |awk '{printf("%.0f", $1 / $2);}'`"
+cols="`echo $ew_dist $grid |awk '{printf("%.0f", $1 / $2);}'`"
+
+########################################################################
+
+north=$maxnorth
+west=$maxwest
+south="`echo $north $grid |awk '{printf("%f", $1 - $2);}'`"
+east="`echo $west $grid |awk '{printf("%f", $1 + $2);}'`"
+
+# number of region
+no_of_region="0"
+
+### rows N-S
+while [ `echo $south $maxsouth |awk '{printf("%d", $1 >= $2);}'` = 1 ];
+do 
+    echo "north -> south"    # 
+    # columns W-E
+    while [ `echo $east $maxeast |awk '{printf("%d", $1<= $2);}'` = 1 ];
+    do
+        echo "west -> east"  # 
+       
+        g.region n=$north s=$south w=$west e=$east
+
+        dx="`echo $east $west |awk '{printf("%f", $1 - $2);}'`"
+        dy="`echo $north $south |awk '{printf("%f", $1 - $2);}'`"
+        coord_x="`echo $west $dx |awk '{printf("%f", $1 + $2);}'`"
+        coord_y="`echo $north $dy |awk '{printf("%f", $1 - $2);}'`"
+
+        planarea="`r.surf.area input=$elev | grep Current | sed -e s/.*:\ *//`"
+        realarea="`r.surf.area input=$elev | grep Estimated | sed -e s/.*:\ *//`"
+
+echo "$coord_x $coord_y $realarea $planarea" | awk '{printf "%d %d %f\n", $1, $2, $3 / $4}'>> $TMP_ascii
+
+        west=$east
+        east="`echo $west $grid |awk '{printf("%f", $1 + $2);}'`"
+        
+        no_of_region="`echo $no_of_region |awk '{printf("%.0f", $1 + 1);}'`"
+        regions_total="`echo $rows $cols |awk '{printf("%.0f", $1 * $2);}'`"
+        echo "--------- REGION NUMBER $no_of_region OF $regions_total ----------"
+
+    done
+    
+    north=$south
+    south="`echo $north $grid |awk '{printf("%f", $1 - $2);}'`";
+    
+    # go west
+    west=$maxwest
+    east="`echo $west $grid |awk '{printf("%f", $1 + $2);}'`";
+done;
+
+
+
+g.region region=oldregion;
+v.in.ascii input=$TMP_ascii output=TMP_vect format=point fs=space skip=0 x=1 y=2 z=3 cat=0 -z;
+v.surf.rst input=TMP_vect layer=0 elev=$ROUGHNESS zmult=1.0 tension=$GIS_OPT_TENSION smooth=$GIS_OPT_SMOOTH --overwrite;
+
+r.colors "$ROUGHNESS" color=rainbow
+
+# record metadata
+r.support "$ROUGHNESS" title="Relief roughness of \"$GIS_OPT_MAP\"" history=""
+r.support "$ROUGHNESS" history="grid size: $grid"
+
+
+### cleaning
+cleanup
+
+echo ""
+if [ -n "$GIS_OPT_ROUGH" ] ; then
+    echo "Surface roughness map created and named [$ROUGHNESS]."
+else
+    echo "Surface roughness map created and named [$ROUGHNESS]. Consider renaming."
+fi
+echo "Done."
+exit 0


Property changes on: grass-addons/raster/r.roughness/r.roughness.sh
___________________________________________________________________
Added: svn:executable
   + *

Added: grass-addons/raster/r.roughness/r.roughness.window.area
===================================================================
--- grass-addons/raster/r.roughness/r.roughness.window.area	                        (rev 0)
+++ grass-addons/raster/r.roughness/r.roughness.window.area	2009-08-11 15:05:22 UTC (rev 38695)
@@ -0,0 +1,186 @@
+#!/bin/sh
+#
+############################################################################
+#
+# MODULE:	r.roughness
+# AUTHOR(S):	Carlos H. Grohmann <carlos dot grohmann at gmail dot com >
+# PURPOSE:	Calculates surface roughness from DEMs. (uses r.surf.area)
+#		In this script surface roughness is used in the sense of 
+#		Hobson (1972), who describes it as the ratio between surface 
+#		(real) area and flat (plan) area of square cells; in this 
+#		approach, flat surfaces would present values close to 1, 
+#		whilst in irregular ones the ratio shows a curvilinear 
+#		relationship which asymptotically approaches infinity as the 
+#		real areas increases.
+#		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 a map of surface area for each pixel
+#		using slope and trigonometry. A map of plan area is create as
+#		ns_res*ew_res. Using the specified moving window size, maps of
+#		the sum of areas are created, and then the ratio is calculated
+#		with r.mapcalc.
+#
+#		If the user does not specify the output map name, it will be
+#		set to INPUT_MAP_roughness_NxN
+#		where N is the window size.
+#
+# COPYRIGHT:	(C) 2007 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
+#%End
+#%option
+#% key: map
+#% gisprompt: old,cell,raster
+#% type: string
+#% description: Input raster elevation map
+#% required : yes
+#%end
+#%option
+#% key: slope
+#% gisprompt: old,cell,raster
+#% type: string
+#% description: Input raster slope map
+#% required : yes
+#%end
+#%option
+#% key: window
+#% type: integer
+#% description: Window size (3,5,7,...,25)
+#% required : no
+#% answer : 3
+#%end
+#%option
+#% key: rough
+#% gisprompt: old,cell,raster
+#% type: string
+#% description: Output raster map 
+#% required : no
+#%end
+
+
+if test "$GISBASE" = ""; then
+ echo "You must be in GRASS GIS to run this program." >&2
+ exit 1
+fi   
+
+# save command line
+if [ "$1" != "@ARGS_PARSED@" ] ; then
+    CMDLINE="`basename $0`"
+    for arg in "$@" ; do
+        CMDLINE="$CMDLINE \"$arg\""
+    done
+    export CMDLINE
+    exec g.parser "$0" "$@"
+fi
+PROG=`basename $0`
+
+# setting environment, so that awk works properly in all languages
+unset LC_ALL
+LC_NUMERIC=C
+export LC_NUMERIC
+
+eval `g.gisenv`
+: ${GISBASE?} ${GISDBASE?} ${LOCATION_NAME?} ${MAPSET?}
+LOCATION=$GISDBASE/$LOCATION_NAME/$MAPSET
+
+program=`basename $0`
+
+map=$GIS_OPT_MAP
+slope=$GIS_OPT_SLOPE
+window=$GIS_OPT_WINDOW
+
+
+#check if input file exists
+
+eval `g.findfile element=cell file=$map`
+if [ -z "$name" ] ; then
+   g.message -e  "Map <$map> not found! Aborting."
+   exit 1
+fi
+
+eval `g.findfile element=cell file=$slope`
+if [ -z "$name" ] ; then
+   g.message -e  "Map <$slope> not found! Aborting."
+   exit 1
+fi
+
+if [ "$GIS_OPT_MAP" = "$GIS_OPT_ROUGH" ]; then
+	echo ""
+	echo "Input elevation map and output roughness map must have different names"
+	exit 1
+fi
+
+if [ -z "$GIS_OPT_ROUGH" ]; then
+    ROUGHNESS="${map%%@*}_area_ratio_${window}x${window}"
+else
+    ROUGHNESS="$GIS_OPT_ROUGH"
+fi
+
+
+#######################################################################
+# what to do in case of user break:
+exitprocedure()
+{
+ g.message -e 'User break!'
+ #delete any TMP files:
+ g.remove rast=$TMP1,$TMP2,$TMP3,$TMP4 > /dev/null
+
+}
+# shell check for user break (signal list: trap -l)
+trap "exitprocedure" 2 3 15
+
+
+#######################################################################
+
+########################################################################
+# get region resolution
+nsres="`g.region -p | grep nsres | sed -e s/.*:\ *//`"
+ewres="`g.region -p | grep ewres | sed -e s/.*:\ *//`"
+
+########################################################################
+
+
+#make temp rasters:
+
+TMP1=surfarea_$$
+TMP2=planarea_$$
+TMP3=sum_surfarea_$$
+TMP4=sum_planarea_$$
+
+r.mapcalc $TMP1 = "$nsres*($nsres/cos($slope))"
+r.mapcalc $TMP2 = "$nsres*$ewres"
+
+r.neighbors -a input=$TMP1 output=$TMP3 method=sum size=$window
+r.neighbors -a input=$TMP2 output=$TMP4 method=sum size=$window
+
+r.mapcalc $ROUGHNESS = "$TMP3 / $TMP4"
+
+r.colors "$ROUGHNESS" color=rainbow
+
+# record metadata
+r.support "$ROUGHNESS" title="Relief roughness of \"$GIS_OPT_MAP\"" history=""
+r.support "$ROUGHNESS" history="grid size: $grid"
+
+
+#cleanup
+g.remove rast=$TMP1,$TMP2,$TMP3,$TMP4 > /dev/null
+
+
+echo ""
+if [ -n "$GIS_OPT_ROUGH" ] ; then
+    echo "Surface roughness map created and named [$ROUGHNESS]."
+else
+    echo "Surface roughness map created and named [$ROUGHNESS]. Consider renaming."
+fi
+echo "Done."
+exit 0


Property changes on: grass-addons/raster/r.roughness/r.roughness.window.area
___________________________________________________________________
Added: svn:executable
   + *

Added: grass-addons/raster/r.roughness/r.roughness.window.vector
===================================================================
--- grass-addons/raster/r.roughness/r.roughness.window.vector	                        (rev 0)
+++ grass-addons/raster/r.roughness/r.roughness.window.vector	2009-08-11 15:05:22 UTC (rev 38695)
@@ -0,0 +1,308 @@
+#!/bin/sh
+#
+############################################################################
+#
+# MODULE:	r.roughness
+# AUTHOR(S):	Carlos H. Grohmann <carlos dot grohmann at gmail dot com >
+# PURPOSE:	Calculates surface roughness from DEMs.
+#		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) 2007 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
+#%End
+#%option
+#% key: map
+#% gisprompt: old,cell,raster
+#% type: string
+#% description: Input raster elevation map
+#% required : yes
+#%end
+#%option
+#% key: slope
+#% gisprompt: old,cell,raster
+#% type: string
+#% description: Input raster slope map
+#% required : yes
+#%end
+#%option
+#% key: aspect
+#% gisprompt: old,cell,raster
+#% type: string
+#% 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
+#% gisprompt: old,cell,raster
+#% type: string
+#% description: Output vector strength map 
+#% required : no
+#%end
+#%option
+#% key: fisher
+#% gisprompt: old,cell,raster
+#% type: string
+#% description: Output Fischer's K parameter map 
+#% required : no
+#%end
+#%option
+#% key: compass
+#% gisprompt: old,cell,raster
+#% type: string
+#% description: OPTIONAL: Compass aspect values (longitude)
+#% required : no
+#%end
+#%option
+#% key: colatitude
+#% gisprompt: old,cell,raster
+#% type: string
+#% description: OPTIONAL: Colatitude values (90 - slope)
+#% required : no
+#%end
+#%option
+#% key: xcos
+#% gisprompt: old,cell,raster
+#% type: string
+#% description: OPTIONAL: X directional cosine map
+#% required : no
+#%end
+#%option
+#% key: ycos
+#% gisprompt: old,cell,raster
+#% type: string
+#% description: OPTIONAL: Y directional cosine map
+#% required : no
+#%end
+#%option
+#% key: zcos
+#% gisprompt: old,cell,raster
+#% type: string
+#% description: OPTIONAL: Z directional cosine map
+#% required : no
+#%end
+#
+if test "$GISBASE" = ""; then
+ echo "You must be in GRASS GIS to run this program." >&2
+ exit 1
+fi   
+
+# save command line
+if [ "$1" != "@ARGS_PARSED@" ] ; then
+    CMDLINE="`basename $0`"
+    for arg in "$@" ; do
+        CMDLINE="$CMDLINE \"$arg\""
+    done
+    export CMDLINE
+    exec g.parser "$0" "$@"
+fi
+PROG=`basename $0`
+
+# setting environment, so that awk works properly in all languages
+unset LC_ALL
+LC_NUMERIC=C
+export LC_NUMERIC
+
+eval `g.gisenv`
+: ${GISBASE?} ${GISDBASE?} ${LOCATION_NAME?} ${MAPSET?}
+LOCATION=$GISDBASE/$LOCATION_NAME/$MAPSET
+
+program=`basename $0`
+
+map=$GIS_OPT_MAP
+slope=$GIS_OPT_SLOPE
+aspect=$GIS_OPT_ASPECT
+window=$GIS_OPT_WINDOW
+
+
+#check if input file exists
+
+eval `g.findfile element=cell file=$map`
+if [ -z "$name" ] ; then
+   g.message -e  "Map <$map> not found! Aborting."
+   exit 1
+fi
+
+eval `g.findfile element=cell file=$slope`
+if [ -z "$name" ] ; then
+   g.message -e  "Map <$slope> not found! Aborting."
+   exit 1
+fi
+
+eval `g.findfile element=cell file=$aspect`
+if [ -z "$name" ] ; then
+   g.message -e  "Map <$aspect> not found! Aborting."
+   exit 1
+fi
+
+if [ -z "$GIS_OPT_STRENGTH" ]; then
+    STRENGTH="${map%%@*}_vector_strength_${window}x${window}"
+else
+    STRENGTH="$GIS_OPT_STRENGTH"
+fi
+
+if [ -z "$GIS_OPT_FISHER" ]; then
+    FISHER="${map%%@*}_fisher_K_${window}x${window}"
+else
+    FISHER="$GIS_OPT_FISHER"
+fi
+
+#######################################################################
+# what to do in case of user break:
+exitprocedure()
+{
+ g.message -e 'User break!'
+ #delete any TMP files:
+ g.remove rast=$TMP1,$TMP3,$TMP4,$TMP5,$TMP6,$TMP7,$TMP8 > /dev/null
+
+}
+# shell check for user break (signal list: trap -l)
+trap "exitprocedure" 2 3 15
+
+
+#######################################################################
+#some temp rasters:
+
+TMP6=sum_dir_cosine_x_$$
+TMP7=sum_dir_cosine_y_$$
+TMP8=sum_dir_cosine_z_$$
+
+#correct aspect angles from cartesian (GRASS default) to compass angles
+#   if(A==0,0,if(A < 90, 90-A, 360+90-A))
+if [ -z "$GIS_OPT_COMPASS" ]; then
+    echo "Calculating compass aspect values (longitude)"
+    TMP1=aspect_compass_$$
+    r.mapcalc $TMP1 = "if($aspect==0,0,if($aspect < 90, 90-$aspect, 360+90-$aspect))"
+else
+    echo "Using previous calculated compass aspect values (longitude)"
+    TMP1="$GIS_OPT_COMPASS"
+fi
+
+#Calculates colatitude (90-slope)
+if [ -z "$GIS_OPT_COLATITUDE" ]; then
+   echo "Calculating colatitude (90-slope)"
+   TMP2=normal_angle_$$
+   r.mapcalc $TMP2 = "90 - $slope"
+else
+    echo "Using previous calculated colatitude values"
+    TMP2="$GIS_OPT_COLATITUDE"
+fi
+
+#direction cosines relative to axis oriented north, east and up
+#direction cosine calculation according to McKean & Roering (2004), Geomorphology, 57:331-351.
+
+if [ -z "$GIS_OPT_XCOS" ]; then
+echo "Calculating X direction cosine"
+    TMP3=dir_cosine_x_$$
+    r.mapcalc $TMP3 = "sin($TMP2)*cos($TMP1)"
+else
+    echo "Using previous calculated X directional cosine values"
+    TMP3="$GIS_OPT_XCOS"
+fi
+
+if [ -z "$GIS_OPT_YCOS" ]; then
+echo "Calculating Y direction cosine"
+    TMP4=dir_cosine_y_$$
+    r.mapcalc $TMP4 = "sin($TMP2)*sin($TMP1)"
+else
+    echo "Using previous calculated Y directional cosine values"
+    TMP4="$GIS_OPT_YCOS"
+fi
+
+if [ -z "$GIS_OPT_ZCOS" ]; then
+echo "Calculating Z direction cosine"
+    TMP5=dir_cosine_z_$$
+    r.mapcalc $TMP5 = "cos($TMP2)"
+else
+    echo "Using previous calculated Z directional cosine values"
+    TMP5="$GIS_OPT_ZCOS"
+fi  
+
+
+echo "Calculating sum of X direction cosines"
+r.neighbors input=$TMP3 output=$TMP6 method=sum size=$window --overwrite
+
+echo "Calculating sum of Y direction cosines"
+r.neighbors input=$TMP4 output=$TMP7 method=sum size=$window --overwrite
+
+echo "Calculating sum of Z direction cosines"
+r.neighbors input=$TMP5 output=$TMP8 method=sum size=$window --overwrite
+
+echo "Calculating vector strength"
+r.mapcalc $STRENGTH = "sqrt(exp($TMP6,2) + exp($TMP7,2) + exp($TMP8,2))"
+
+
+echo "Calculating Inverted Fisher's K parameter"
+# k=1/((N-1)/(N-R))
+r.mapcalc $FISHER = "($window * $window - $STRENGTH) / ($window * $window - 1)"
+
+#cleanup
+g.remove rast=$TMP6,$TMP7,$TMP8 > /dev/null
+
+if [ -z "$GIS_OPT_COMPASS" ]; then
+g.remove rast=$TMP1 > /dev/null
+fi
+
+if [ -z "$GIS_OPT_COLATITUDE" ]; then
+g.remove rast=$TMP2 > /dev/null
+fi
+
+if [ -z "$GIS_OPT_XCOS" ]; then
+g.remove rast=$TMP3 > /dev/null
+fi
+
+if [ -z "$GIS_OPT_YCOS" ]; then
+g.remove rast=$TMP4 > /dev/null
+fi
+
+if [ -z "$GIS_OPT_ZCOS" ]; then
+g.remove rast=$TMP5 > /dev/null
+fi
+
+echo ""
+if [ -n "$GIS_OPT_STRENGTH" ] ; then
+    echo "Surface roughness map created and named [$STRENGTH]."
+else
+    echo "Surface roughness map created and named [$STRENGTH]. Consider renaming."
+fi
+
+echo ""
+if [ -n "$GIS_OPT_FISHER" ] ; then
+    echo "Surface roughness map created and named [$FISHER]."
+else
+    echo "Surface roughness map created and named [$FISHER]. Consider renaming."
+fi
+
+
+echo "Done."
+exit 0
+
+


Property changes on: grass-addons/raster/r.roughness/r.roughness.window.vector
___________________________________________________________________
Added: svn:executable
   + *

Added: grass-addons/raster/r.roughness/r.roughness.window.vector.html
===================================================================
--- grass-addons/raster/r.roughness/r.roughness.window.vector.html	                        (rev 0)
+++ grass-addons/raster/r.roughness/r.roughness.window.vector.html	2009-08-11 15:05:22 UTC (rev 38695)
@@ -0,0 +1,50 @@
+<!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.window.vector</b></em> - Calculates surface roughness in a moving-window, as the orientation of vectors normal to surface planes.
+<h2>KEYWORDS&nbsp;</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&nbsp;help</span><br>
+<span style="font-weight: bold;">r.roughness.window.vector</span>&nbsp;<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&nbsp;</dd><dt><b>ycos</b>=<em>string</em><em> (optional)</em></dt><dd>Input <span style="font-style: italic;">y</span> directional cosine map&nbsp;</dd><dt><b>zcos</b>=<em>string</em><em> (optional)</em></dt><dd>Input <span style="font-style: italic;">z</span> directional cosine map&nbsp;</dd></dl><h2>DESCRIPTION</h2>
+In this script surface roughness is taken as the dispersion of vectors
+normal to surface areas (pixels). Normal vectors&nbsp;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 &lt; 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 -&nbsp;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: italic;"><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. &amp; 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>AUTHOR</h2>Carlos Henrique Grohmann<br>Institute of Geosciences, University of São Paulo, Brazil.
+<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



More information about the grass-commit mailing list