[GRASS-SVN] r59939 - grass-addons/grass6/raster/r.roughness
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu Apr 24 14:58:28 PDT 2014
Author: neteler
Date: 2014-04-24 14:58:27 -0700 (Thu, 24 Apr 2014)
New Revision: 59939
Added:
grass-addons/grass6/raster/r.roughness/r.roughness.window.area.sh
grass-addons/grass6/raster/r.roughness/r.roughness.window.vector.sh
Removed:
grass-addons/grass6/raster/r.roughness/r.roughness.window.area
grass-addons/grass6/raster/r.roughness/r.roughness.window.vector
Modified:
grass-addons/grass6/raster/r.roughness/Makefile
grass-addons/grass6/raster/r.roughness/r.roughness.sh
Log:
r.roughness.window.vector: installation fix
Modified: grass-addons/grass6/raster/r.roughness/Makefile
===================================================================
--- grass-addons/grass6/raster/r.roughness/Makefile 2014-04-24 14:11:03 UTC (rev 59938)
+++ grass-addons/grass6/raster/r.roughness/Makefile 2014-04-24 21:58:27 UTC (rev 59939)
@@ -1,6 +1,6 @@
MODULE_TOPDIR = ../..
-PGM=r.roughness.window.vector
+PGM=r.roughness.window.vector.sh
include $(MODULE_TOPDIR)/include/Make/Script.make
Modified: grass-addons/grass6/raster/r.roughness/r.roughness.sh
===================================================================
--- grass-addons/grass6/raster/r.roughness/r.roughness.sh 2014-04-24 14:11:03 UTC (rev 59938)
+++ grass-addons/grass6/raster/r.roughness/r.roughness.sh 2014-04-24 21:58:27 UTC (rev 59939)
@@ -2,7 +2,7 @@
#
############################################################################
#
-# MODULE: r.roughness
+# MODULE: r.roughness.sh
# 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
@@ -38,7 +38,7 @@
#% key: map
#% gisprompt: old,cell,raster
#% type: string
-#% description: Input raster surface
+#% description: Input raster elevation map
#% required : yes
#%end
#%option
@@ -85,7 +85,6 @@
export CMDLINE
exec g.parser "$0" "$@"
fi
-
PROG=`basename $0`
@@ -252,5 +251,6 @@
else
echo "Surface roughness map created and named [$ROUGHNESS]. Consider renaming."
fi
+
echo "Done."
exit 0
Deleted: grass-addons/grass6/raster/r.roughness/r.roughness.window.area
===================================================================
--- grass-addons/grass6/raster/r.roughness/r.roughness.window.area 2014-04-24 14:11:03 UTC (rev 59938)
+++ grass-addons/grass6/raster/r.roughness/r.roughness.window.area 2014-04-24 21:58:27 UTC (rev 59939)
@@ -1,186 +0,0 @@
-#!/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
Copied: grass-addons/grass6/raster/r.roughness/r.roughness.window.area.sh (from rev 59870, grass-addons/grass6/raster/r.roughness/r.roughness.window.area)
===================================================================
--- grass-addons/grass6/raster/r.roughness/r.roughness.window.area.sh (rev 0)
+++ grass-addons/grass6/raster/r.roughness/r.roughness.window.area.sh 2014-04-24 21:58:27 UTC (rev 59939)
@@ -0,0 +1,187 @@
+#!/bin/sh
+#
+############################################################################
+#
+# MODULE: roughness.window.area.sh
+# 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
Deleted: grass-addons/grass6/raster/r.roughness/r.roughness.window.vector
===================================================================
--- grass-addons/grass6/raster/r.roughness/r.roughness.window.vector 2014-04-24 14:11:03 UTC (rev 59938)
+++ grass-addons/grass6/raster/r.roughness/r.roughness.window.vector 2014-04-24 21:58:27 UTC (rev 59939)
@@ -1,308 +0,0 @@
-#!/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
-
-
Copied: grass-addons/grass6/raster/r.roughness/r.roughness.window.vector.sh (from rev 59870, grass-addons/grass6/raster/r.roughness/r.roughness.window.vector)
===================================================================
--- grass-addons/grass6/raster/r.roughness/r.roughness.window.vector.sh (rev 0)
+++ grass-addons/grass6/raster/r.roughness/r.roughness.window.vector.sh 2014-04-24 21:58:27 UTC (rev 59939)
@@ -0,0 +1,305 @@
+#!/bin/sh
+#
+############################################################################
+#
+# MODULE: r.roughness.window.vector.sh
+# 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
More information about the grass-commit
mailing list