[GRASS-SVN] r36410 - in grass-addons/LandDyn: devs_landcover_scripts
devs_landcover_scripts/r.landcover.update
devs_landcover_scripts/r.soil.fertility
devs_landcover_scripts/rules r.landscape.evol r.soildepth
svn_grass at osgeo.org
svn_grass at osgeo.org
Wed Mar 18 18:11:19 EDT 2009
Author: isaacullah
Date: 2009-03-18 18:11:19 -0400 (Wed, 18 Mar 2009)
New Revision: 36410
Added:
grass-addons/LandDyn/devs_landcover_scripts/r.soil.fertility/
grass-addons/LandDyn/devs_landcover_scripts/r.soil.fertility/r.soil.fertility
grass-addons/LandDyn/devs_landcover_scripts/rules/sfertil_colors.txt
Modified:
grass-addons/LandDyn/devs_landcover_scripts/r.landcover.update/r.landcover.update
grass-addons/LandDyn/r.landscape.evol/r.landscape.evol
grass-addons/LandDyn/r.soildepth/r.soildepth
Log:
Major updates to core scripts.
Modified: grass-addons/LandDyn/devs_landcover_scripts/r.landcover.update/r.landcover.update
===================================================================
--- grass-addons/LandDyn/devs_landcover_scripts/r.landcover.update/r.landcover.update 2009-03-18 22:02:21 UTC (rev 36409)
+++ grass-addons/LandDyn/devs_landcover_scripts/r.landcover.update/r.landcover.update 2009-03-18 22:11:19 UTC (rev 36410)
@@ -36,6 +36,22 @@
#%END
#%option
+#% key: sfertil
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: map of current soil fertility
+#% required : yes
+#%END
+
+#%option
+#% key: sdepth
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: map of current soil depths
+#% required : yes
+#%END
+
+#%option
#% key: max
#% type: string
#% gisprompt: string
@@ -94,6 +110,10 @@
impacts=$GIS_OPT_impacts
+sfertil=$GIS_OPT_sfertil
+
+sdepth=$GIS_OPT_sdepth
+
outmap=$GIS_OPT_outmap
max=$GIS_OPT_max
@@ -105,6 +125,8 @@
txtout=$outmap"luse_stats.txt"
+temp_rate="temp_rate"
+
temp_lc=$outmap"temp_landcover"
temp_reclass=$outmap"temp_landcover_reclass"
@@ -118,12 +140,15 @@
r.mask --quiet input=$inmap maskcats=*
-#updating raw landscape category numbers
+# calculating rate of regrowth based on curretn soil fertility and depths. Recoding fertility (0 to 100) and depth (0 to >= 1) with a power regression curve from 0 to 1, then taking the mean of the two as the regrowth rate
-#r.mapcalc "$temp_lc = if ($inmap == $max && isnull($impacts), $max, (if ($inmap < $max && isnull($impacts), ($inmap + 1), ($inmap - $impacts))))"
+r.mapcalc "temp_rate=if($sdepth <= 1, ( ( ( (-0.000118528 * (exp($sfertil,2))) + (0.0215056 * $sfertil) + 0.0237987 ) + ( ( -0.000118528 * (exp($sdepth,2))) + (0.0215056 * $sdepth) + 0.0237987 ) ) / 2 ), ( ( ( (-0.000118528 * (exp($sfertil,2))) + (0.0215056 * $sfertil) + 0.0237987 ) + 1) / 2 ) )"
-r.mapcalc "$temp_lc = if ($inmap == $max && isnull($impacts), $max, (if ($inmap < $max && isnull($impacts), ($inmap + 1), if ($inmap > $max, ($max - $impacts), if ($inmap < 0, 0, ($inmap - $impacts))))))"
+#updating raw landscape category numbers based on agent impacts and newly calculated regrowth rate
+
+r.mapcalc "$temp_lc = if ($inmap == $max && isnull($impacts), $max, (if ($inmap < $max && isnull($impacts), ($inmap + temp_rate, if ($inmap > $max, ($max - $impacts), if ($inmap < 0, 0, ($inmap - $impacts))))))"
+
#adding text descriptions to raw landscape categories and setting colors
cat $lc_rules | r.reclass input=$temp_lc output=$temp_reclass
@@ -150,7 +175,7 @@
echo "" >> $txtout
echo "Landcover class #, Landcover description, Area (sq. m)" >> $txtout
echo "" >> $txtout
-r.stats -a -l -n input=$prfx"_landuse1" fs=, nv=* nsteps=255 >> $txtout
+r.stats -a -l -n input=$prfx"_landuse1" fs=, nv=* nsteps=$max >> $txtout
fi
@@ -163,7 +188,7 @@
echo ""
-g.remove --quiet rast=MASK,$temp_reclass,$temp_lc
+g.remove --quiet rast=MASK,$temp_reclass,$temp_lc,$temp_rate
Added: grass-addons/LandDyn/devs_landcover_scripts/r.soil.fertility/r.soil.fertility
===================================================================
--- grass-addons/LandDyn/devs_landcover_scripts/r.soil.fertility/r.soil.fertility (rev 0)
+++ grass-addons/LandDyn/devs_landcover_scripts/r.soil.fertility/r.soil.fertility 2009-03-18 22:11:19 UTC (rev 36410)
@@ -0,0 +1,153 @@
+#!/bin/sh
+#
+############################################################################
+#
+# MODULE: r.soil.fertility
+# AUTHOR(S): Isaac Ullah, Michael Barton, Arizona State University
+# PURPOSE: Updates soil fertility thru time based on an impacts
+# map created by an agent base model.
+# ACKNOWLEDGEMENTS: National Science Foundation Grant #BCS0410269
+# COPYRIGHT: (C) 2007 by Isaac Ullah, Michael Barton, Arizona State University
+# This program is free software under the GNU General Public
+# License (>=v2). Read the file COPYING that comes with GRASS
+# for details.
+#
+#############################################################################
+
+
+#%Module
+#% description: Updates soil fertility thru time based on an impact map created by an agent base model.
+#%END
+
+#%option
+#% key: inmap
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: input soil fertility map (values coded 0-max)
+#% required : yes
+#%END
+
+#%option
+#% key: impacts
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: map of impacts on soil fertility values
+#% required : yes
+#%END
+
+#%option
+#% key: max
+#% type: string
+#% gisprompt: string
+#% description: maximum value for soil fertility maps (number for highest fertility)
+#% answer: 100
+#% required : yes
+#%END
+
+#%option
+#% key: outmap
+#% type: string
+#% gisprompt: string
+#% description: New soil fertility output map name (no prefix)
+#% answer: landcover
+#% required : yes
+#%END
+
+#%option
+#% key: sf_color
+#% type: string
+#% gisprompt: string
+#% description: path to color rules file for landcover map
+#% answer: /usr/local/grass-6.3.cvs/scripts/rules/sfertil_colors.txt
+#% required : yes
+#%END
+
+#%flag
+#% key: s
+#% description: -s Output text file of soil fertility stats from the simulation (will be named "prefix"_sfertil_stats.txt, and will be overwritten if you run the simulation again with the same prefix)
+#%END
+
+
+if [ -z "$GISBASE" ] ; then
+ echo "You must be in GRASS GIS to run this program." >&2
+ exit 1
+fi
+
+if [ "$1" != "@ARGS_PARSED@" ] ; then
+ exec g.parser "$0" "$@"
+fi
+
+
+
+#setting up variables for use later on
+
+inmap=$GIS_OPT_inmap
+
+impacts=$GIS_OPT_impacts
+
+outmap=$GIS_OPT_outmap
+
+max=$GIS_OPT_max
+
+sf_rules=$GIS_OPT_sf_rules
+
+sf_color=$GIS_OPT_sf_color
+
+
+txtout=$outmap"sfertil_stats.txt"
+
+outsf=$outmap
+
+
+#setting initial conditions of map area
+
+g.region rast=$inmap
+
+r.mask --quiet input=$inmap maskcats=*
+
+#updating raw soil fertility category numbers
+
+r.mapcalc "$out_sf = if ($inmap == $max && isnull($impacts), $max, (if ($inmap < $max && isnull($impacts), ($inmap + 1), if ($inmap > $max, ($max - $impacts), if ($inmap < 0, 0, ($inmap - $impacts))))))"
+
+ r.colors map=$outsf rules=$sf_color
+
+
+#checking total area of updated cells
+
+ temparea=`eval r.stats -n -a fs=- input=$impacts | cut -d'-' -f2`
+ echo "***********************"
+ echo "Total area of impacted zones = $temparea square meters"
+ echo "***********************"
+
+#creating optional output text file of stats
+
+if [ "$GIS_FLAG_s" -eq 1 ]; then
+
+echo "Stats for $prfx'_soil_fertility' " > $txtout
+echo "" >> $txtout
+echo "Total area of impacted zones = $temparea square meters" >> $txtout
+echo "" >> $txtout
+echo "Landcover class #, Landcover description, Area (sq. m)" >> $txtout
+echo "" >> $txtout
+r.stats -a -l -n input=$prfx"_landuse1" fs=, nv=* nsteps=$max >> $txtout
+
+fi
+
+
+
+echo ""
+echo "*************************"
+echo " Cleaning up"
+echo "*************************"
+echo ""
+
+
+g.remove --quiet rast=MASK
+
+
+
+echo ""
+echo "DONE!"
+echo ""
+echo ""
+
Added: grass-addons/LandDyn/devs_landcover_scripts/rules/sfertil_colors.txt
===================================================================
--- grass-addons/LandDyn/devs_landcover_scripts/rules/sfertil_colors.txt (rev 0)
+++ grass-addons/LandDyn/devs_landcover_scripts/rules/sfertil_colors.txt 2009-03-18 22:11:19 UTC (rev 36410)
@@ -0,0 +1,6 @@
+0 white
+20 grey
+40 yellow
+60 orange
+80 brown
+100 black
Modified: grass-addons/LandDyn/r.landscape.evol/r.landscape.evol
===================================================================
--- grass-addons/LandDyn/r.landscape.evol/r.landscape.evol 2009-03-18 22:02:21 UTC (rev 36409)
+++ grass-addons/LandDyn/r.landscape.evol/r.landscape.evol 2009-03-18 22:11:19 UTC (rev 36410)
@@ -115,6 +115,11 @@
#% answer: 1
#% guisection: Input
#%end
+#%flag
+#% key: p
+#% description: -p Keep all slope maps
+#% guisection: Input
+#%end
#%option
@@ -162,6 +167,14 @@
#% guisection: Variables
#%end
#%option
+#% key: alpha
+#% type: string
+#% description: Critical slope threshold for mass movement of sediment (in degrees above horizontal)
+#% answer: 40
+#% required : yes
+#% guisection: Variables
+#%end
+#%option
#% key: cutoff1
#% type: string
#% description: Flow accumultion breakpoint value for shift from diffusion to overland flow (number of cells)
@@ -172,12 +185,20 @@
#%option
#% key: cutoff2
#% type: string
-#% description: Flow accumultion breakpoint value for shift from overland flow to channelized flow (number of cells)
+#% description: Flow accumultion breakpoint value for shift from overland flow to rill/gully flow (number of cells)
#% answer: 50
#% required : yes
#% guisection: Variables
#%end
#%option
+#% key: cutoff3
+#% type: string
+#% description: Flow accumultion breakpoint value for shift from rill/gully flow to stream flow (number of cells)
+#% answer: 10000
+#% required : yes
+#% guisection: Variables
+#%end
+#%option
#% key: number
#% type: integer
#% description: number of iterations (cycles) to run
@@ -315,13 +336,14 @@
sdensity=$GIS_OPT_sdensity
C=$GIS_OPT_C
kappa=$GIS_OPT_kappa
+alpha=$GIS_OPT_alpha
Ba=$GIS_OPT_Ba
Bb=$GIS_OPT_Bb
-aplpha=$GIS_OPT_method
sigma=$GIS_OPT_sigma
nbhood=$GIS_OPT_nbhood
cutoff1=$GIS_OPT_cutoff1
cutoff2=$GIS_OPT_cutoff2
+cutoff3=$GIS_OPT_cutoff3
num_iters=$GIS_OPT_number
echo "total number of iterations to be run= $num_iters"
echo ""
@@ -518,27 +540,16 @@
echo "Amount of free RAM being allocated for this step: $mem Megabites"
- tmpflacc=$prefx".tmpflacc"
-
- r.terraflow --q elev=$old_dem filled=$prefx".filled" direction=$prefx".direct" swatershed=$prefx".sink" accumulation=$tmpflacc tci=$prefx".tci" d8cut=infinity memory=$mem STREAM_DIR=/var/tmp
-
- #changing raw flow accumulation map to number of square meters of upslope contributing area
- r.mapcalc "$flowacc=$tmpflacc*exp($res, 2)"
-
- g.remove --quiet rast=$prefx".filled",$prefx".direct",$prefx".sink",$prefx".tci",$tmpflacc,$tmpdirection
-
+ r.terraflow --quiet elev=$old_dem filled=$prefx".filled" direction=$prefx".direct" swatershed=$prefx".sink" accumulation=$flowacc tci=$prefx".tci" d8cut=infinity memory=$mem STREAM_DIR=/var/tmp
+ g.remove --quiet rast=$prefx".filled",$prefx".direct",$prefx".sink",$prefx".tci",$tmpdirection
rm -f /var/tmp/STREAM*
else
echo "using r.watershed to calculate overland flow accumulation per cell (number of cells uplsope from each cell)"
tmpflacc=$prefx"tmpflacc"
-
- r.watershed --quiet elevation=$old_dem visual=$tmpflacc memory=$mem
-
- #cchanging raw flow accumulation map to number of square meters of upslope contributing area
- r.mapcalc "$flowacc=$tmpflacc*exp($res, 2)"
-
+ r.watershed --quiet -f elevation=$old_dem accumulation=$tmpflacc convergence=5
+ r.mapcalc "$flowacc=abs($tmpflacc)"
g.remove --quiet rast=$tmpflacc
fi
@@ -554,7 +565,7 @@
r.mapcalc "$sflowtopo=$flowacc * sin($slope)"
elif [ "$GIS_FLAG_r" -eq 1 -a "$GIS_FLAG_d" -eq 0 -a "$GIS_FLAG_w" -eq 0 ]; then
echo ""
- echo "calculating for channelzed flow across the entire map"
+ echo "calculating for rill/gully flow across the entire map"
r.mapcalc "$sflowtopo=exp(($flowacc),1.6) * exp(sin($slope),1.3)"
@@ -562,7 +573,7 @@
echo ""
echo "calculating for diffusive flow across the entire map"
- r.mapcalc "$sflowtopo=$kappa * exp(($res),2) * sin($slope)"
+ r.mapcalc "$sflowtopo=$kappa * sin($slope)"
elif [ "$GIS_FLAG_d" -eq 1 -a "$GIS_FLAG_r" -eq 1 -a "$GIS_FLAG_w" -eq 0 ]; then
echo ""
@@ -596,16 +607,26 @@
maxflow=`eval r.info -r map=$flowacc | grep "max=" | cut -d"=" -f2`
echo ""
- echo "The raw max ($maxflow) flow accumulation, min cutoff $cutoff1, and middle cutoff $cutoff2 will be used to scale flow exponents 'm' and 'n' as flow progresses in the drainage network"
+ echo "The raw max ($maxflow) flow accumulation, creep cutoff $cutoff1, rill cutoff $cutoff2, and stream cutoff $cutoff3 will be used to scale flow exponents 'm' and 'n' as flow progresses in the drainage network"
echo ""
echo ""
echo "calculating with appropriate flow types on different landforms"
- r.mapcalc "$m=if($flowacc > $cutoff2, 1.6, if($flowacc > $cutoff1 && $flowacc < $cutoff2, (1+($flowacc * (0.6/($cutoff2 - $cutoff1)))), 1))"
+ # first we have to figurte out the correct value for m and n depending on the flow accumulation value at a cell, which is used here as an approximation for the position of that cell in the drainage network of the watershed. user defiend cutoffs are used as inflexion points (where the slope of the linear equation changes) in the otherwise linear growth of m and n from their minumum overland flow values (both 1) to their maximum values (2.5 and 2). below the first cutoff, m and n are 0 and 1, which is approximates soil creep)
- r.mapcalc "$n=if($flowacc > $cutoff2,1.3, if($flowacc > $cutoff1 && $flowacc < $cutoff2, (1+($flowacc * (0.3/($cutoff2 - $cutoff1)))), 1))"
+ r.mapcalc "$m=if($flowacc > $cutoff3, 2.5000000, if($flowacc > $cutoff2 && $flowacc <= $cutoff3, ( 1.6000000 + ( ($flowacc - $cutoff2) * (0.9000000 / ($cutoff3 - $cutoff2) ) ) ), if($flowacc > $cutoff1 && $flowacc <= $cutoff2, ( 1.0000000 + ( ($flowacc - $cutoff1) * (0.60000000 / ($cutoff2 - $cutoff1) ) ) ), 0.0000000) ) )"
- r.mapcalc "$sflowtopo=if ($flowacc >= $cutoff1, exp(($flowacc),$m) * exp(sin($slope),$n), ($kappa * exp(($res),2) *sin($slope)))"
+ r.mapcalc "$n=if($flowacc > $cutoff3, 2.0000000, if($flowacc > $cutoff2 && $flowacc <= $cutoff3, ( 1.3000000 + ( ($flowacc - $cutoff2) * ( 0.7000000 / ($cutoff3 - $cutoff2) ) ) ), if($flowacc > $cutoff1 && $flowacc <= $cutoff2, ( 1.0000000 + ( ($flowacc - $cutoff1) * ( 0.3000000 / ($cutoff2 - $cutoff1) ) ) ), 1.0000000) ) )"
+
+ # now that m and n are defined for each cell, we use the continuity equation to approximate transport capacity at each pointin the landscape. alpha is our critical threshold for mass movement of materials, so if slope exceeds this, we use the talus slop/landslide variant of the equation. If flow accumulation is low, we use the soil creep version of the equation by usign the diffusion coefficient kappa (in these areas, m=0 and n=1, as calculated above). For the rest of the drainage network we used the calculated values of m and n so that as flow accumulation increases we move from m and n valuse for sheetwash to rill/gully flow to stream flow.
+# with res mult
+# r.mapcalc "$sflowtopo=if($slope > $alpha, (exp(($flowacc*$res),$m) * exp(((tan($slope)) - (tan($alpha))),$n)), if($flowacc >= $cutoff1, (exp(($flowacc*$res),$m) * exp(sin($slope),$n)), ($kappa * exp(sin($slope),$n)) ))"
+
+# without alpha test but with res mult
+ r.mapcalc "$sflowtopo=if($flowacc >= $cutoff1, (exp(($flowacc*$res),$m) * exp(sin($slope),$n)), ($kappa * exp(sin($slope),$n)) )"
+
+# without res mult
+# r.mapcalc "$sflowtopo=if($slope > $alpha, (exp(($flowacc),$m) * exp(((tan($slope)) - (tan($alpha))),$n)), if($flowacc >= $cutoff1, (exp(($flowacc),$m) * exp(sin($slope),$n)), ($kappa * exp(sin($slope),$n)) ))"
fi
echo ""
@@ -614,10 +635,17 @@
echo "*************************"
echo ""
- #if the old soil is 0 or negative(could happen i suppose) then we make k-factor really really small to simulate erosion on bedrock
+
+# changed
+# #if the old soil is 0 or negative(could happen i suppose) then we make k-factor really really small to simulate erosion on bedrock
r.mapcalc "$qsx=if ($old_soil <= 0, ($R * 0.001 * $C * $sflowtopo * cos($aspect)), ($R * $K * $C * $sflowtopo * cos($aspect)))"
r.mapcalc "$qsy=if ($old_soil <= 0, ($R * 0.001 * $C * $sflowtopo * sin($aspect)), ($R * $K * $C * $sflowtopo * sin($aspect)))"
+#r.mapcalc "$qsx=if ($old_soil <= 0, ($R * 0.001 * $C * $sflowtopo), ($R * $K * $C * $sflowtopo))"
+#r.mapcalc "$qsy=if ($old_soil <= 0, ($R * 0.001 * $C * $sflowtopo), ($R * $K * $C * $sflowtopo))"
+ #r.mapcalc "$qsx=$sflowtopo * cos($aspect)"
+ #r.mapcalc "$qsy=$sflowtopo * sin($aspect)"
+
echo ""
echo "*************************"
echo "step 5 of 8: calculating partial derivatives for sediment transport"
@@ -629,14 +657,24 @@
echo ""
echo "*************************"
- echo "step 6 of 8: calculating net erosion and deposition"
+ echo "step 6 of 8: calculating net erosion and deposition in meters of vertical change per cell"
echo "*************************"
echo ""
+
+ #Now we multiply R K and C to get some local conditions for the final erosion/deposition rate. If the old soil is 0 or negative(could happen i suppose) then we make k-factor really really small to simulate erosion on bedrock
+ #Then we add the x and y flux's for total flux "T" in Mg/Ha. Then if we multiply T times 100, we get a rate in grams/squaremeter. This new rate times the resolution gives us grams per cell width. Then we divide this by the density to get cubic centimeters per cell width. This measure is then divided by the area of the cell in centimeters (resolution squared x 100 squared) to get vertical change per cell width in centimeters. We dived this by 100 to get that measure in meters. Several of these facctors cancel out to make a final equation of "T/(10,000*density*resolution)". This equation changes the orignal T from Mg/Ha to vertical change in meters over the length of one cell's width. In order to convert the output back to Mg/Ha (standard rate for USPED/RUSLE equations), you can multiply the netchange output map by "(10000 x resolution x soil density)" to create a map of soil erosion/deposition rates across the landscape. The rest of this mampcalc statement just checks the amount of erodable soil in a given cell against the amount of erosion calculated, and keeps the cell from eroding past this amount (if there is soil, then if the amount of erosion is more than the amount of soil, just remove all the soil and stop, else remove the amount of caclulated erosion. It also runs an error catch that checks to make sure that soil depth is not negative, and if it is, corrects it).
- #Here we add the x and y flux's for total flux "T" in Mg/Ha. Then if we multiply T times 100, we get a rate in grams/squaremeter. This new rate times the resolution gives us grams per cell width. Then we divide this by the density to get cubic centimeters per cell width. This measure is then divided by the area of the cell in centimeters (resolution squared x 100 squared) to get vertical change per cell width in centimeters. We dived this by 100 to get that measure in meters. Several of these facctors cancel out to make a final equation of "T/(10,000*density*resolution)". This equation changes the orignal T from Mg/Ha to vertical change in meters over the length of one cell's width. In order to convert the output back to Mg/Ha (standard rate for USPED/RUSLE equations), you can multiply the netchange output map by "(10000 x resolution x soil density)" to create a map of soil erosion/deposition rates across the landscape. The rest of this mampcalc statement just checks the amount of erodable soil in a given cell against the amount of erosion calculated, and keeps the cell from eroding past this amount (if there is soil, then if the amount of erosion is more than the amount of soil, just remove all the soil and stop, else remove the amount of caclulated erosion. It also runs an error catch that checks to make sure that soil depth is not negative, and if it is, corrects it).
- r.mapcalc "$erdep=if( $old_soil >= 0, if( (-1 * ( ($qsxdx + $qsydy) / ( 10000 * $sdensity * $res ) ) ) >= $old_soil, ( -1 * $old_soil ), ( ($qsxdx + $qsydy) / ( 10000 * $sdensity * $res ) ) ), ( ($qsxdx + $qsydy) / ( 10000 * $sdensity * $res ) ) )"
+# OLD
+# r.mapcalc "$erdep= if( $old_soil >= 0, if( (-1 * ( ( ($EDx + $EDy) * $res ) / ( 10000 * $sdensity ) ) ) >= $old_soil, ( -1 * $old_soil ), ( ( ($EDx + $EDy) * $res )/ ( 10000 * $sdensity ) ) ), ( ( ($EDx + $EDy) * $res ) / ( 10000 * $sdensity ) ) )"
+#no res mult
+ r.mapcalc "$erdep= if( $old_soil >= 0, if( (-1 * ( ( ($qsxdx + $qsydy) * $res ) / ( 10000 * $sdensity ) ) ) >= $old_soil, ( -1 * $old_soil ), ( ( ($qsxdx + $qsydy) * $res )/ ( 10000 * $sdensity ) ) ), ( ( ($qsxdx + $qsydy) * $res ) / ( 10000 * $sdensity ) ) )"
+
+# res mult
+# r.mapcalc "$erdep= if( $old_soil >= 0, if( (-1 * ( ( ($qsxdx + $qsydy) * $res ) / ( 10000 * $sdensity ) ) ) >= $old_soil, ( -1 * $old_soil ), $res * ( ( ($qsxdx + $qsydy) * $res )/ ( 10000 * $sdensity ) ) ), $res * ( ( ($qsxdx + $qsydy) * $res ) / ( 10000 * $sdensity ) ) )"
+
+
if [ "$GIS_FLAG_y" -eq 1 ]; then
echo "You have selected to smooth your files. Running bandpass despeckling filter now"
echo ""
@@ -905,7 +943,7 @@
echo "************************"
echo ""
-if [ "$GIS_FLAG_n" -ne 1 ]; then
+if [ "$GIS_FLAG_n" -eq 0 ]; then
g.mremove -f --q rast=$prefx"netchange*"
fi
@@ -924,7 +962,12 @@
echo "Cleaning up..."
g.mremove -f --quiet rast=$prefx"flowacc*"
+ if [ "$GIS_FLAG_p" -eq 0 ]; then
g.mremove -f --quiet rast=$prefx"slope*"
+ if [ "$num_iters" -eq 1 ]; then
+ g.rename --quiet rast=$prefx"slope*1",$prefx"slope"
+ fi
+ fi
g.mremove -f --quiet rast=$prefx"aspect*"
Modified: grass-addons/LandDyn/r.soildepth/r.soildepth
===================================================================
--- grass-addons/LandDyn/r.soildepth/r.soildepth 2009-03-18 22:02:21 UTC (rev 36409)
+++ grass-addons/LandDyn/r.soildepth/r.soildepth 2009-03-18 22:11:19 UTC (rev 36410)
@@ -53,13 +53,8 @@
#% answer: 100000
#% required : yes
#%END
-#%flag
-#% key: r
-#% description: -r Keep soil production rate map
-#%END
-
if [ -z "$GISBASE" ] ; then
echo "You must be in GRASS GIS to run this program." >&2
exit 1
@@ -88,7 +83,6 @@
pc="temp.pc.deletable"
tc="temp.tc.deletable"
rate="temp.rate.deletable"
-meancurv="temp.meancurv.deletable"
echo ""
echo "STEP 1, calculating profile and tangental curvatures"
@@ -97,55 +91,30 @@
r.slope.aspect --quiet $elev pcurv=$pc tcurv=$tc
echo ""
-echo "STEP 2, caclulating rate of soil production from mean hillslope curvatures"
+echo "STEP 2, caclulating differential rates of soil mantling from mean hillslope curvatures and kappa value"
echo ""
-r.mapcalc "$meancurv=(($pc+$tc)/2)"
+r.mapcalc "$rate=$kappa*(1-(($pc+$tc)/2))"
-
-max=`eval r.info -r map=$meancurv | grep "max=" | cut -d"=" -f2`
-min=`eval r.info -r map=$meancurv | grep "min=" | cut -d"=" -f2`
-
echo ""
-echo "The raw max ($max) and min ($min) will be rescaled from 2 to 0 "
+echo "STEP 3, calculating soil depth for given amount of time ($t years)"
echo ""
+r.mapcalc "$soildepth = ($rate*$t)"
-r.mapcalc "$rate=$kappa*(2-($meancurv*(2/($max)-($min))))"
+meandepth=`eval r.univar -g map=$soildepth percentile=90 | grep "mean=" | cut -d'=' -f2`
-rate1=`eval r.univar -g map=$rate percentile=90 | grep "mean=" | cut -d'=' -f2`
-
-echo "mean rate = $rate1"
-
echo ""
-echo "STEP 3, calculating soil depth for given amount of time"
+echo "STEP 4, calculating bedrock elevation for given soil depth"
echo ""
-echo "age of landscape = $t"
-echo "output = $soildepth"
+r.mapcalc "$bedrock=($elev - $soildepth)"
-r.mapcalc "$soildepth = ($rate*10000)"
-
-meandepth=`eval r.univar -g map=$soildepth percentile=90 | grep "mean=" | cut -d'=' -f2`
-
-echo "mean depth = $meandepth"
echo ""
-
-if [ "$GIS_FLAG_r" -eq 1 ]; then
-g.remove --q rast=$pc,$tc
-g.rename --o --q rast=$rate,$soildepth"_prod_rate"
-else
-g.remove --q rast=$pc,$tc,$rate,$meancurv
-fi
-
-
-
-
+echo "Cleaning up..."
echo ""
-echo "STEP 4, calculating bedrock elevation for given soil depth"
-echo ""
-r.mapcalc "$bedrock=($elev) - ($soildepth)"
+g.remove --q rast=$pc,$tc,$rate
if [ -n "$GIS_OPT_soildepth" ]; then
echo ""
@@ -153,9 +122,11 @@
echo ""
echo "########################################################"
echo "# #"
+echo "# Mean soil depth = $meandepth #"
echo "# Soil depth map = $soildepth #"
echo "# Bedrock elevation map = $bedrock #"
-
+echo "# #"
+echo "########################################################"
else
echo ""
echo " DONE!"
@@ -163,19 +134,14 @@
echo "########################################################"
echo "# #"
echo "# #"
+echo "# Mean soil depth = $meandepth #"
echo "# Bedrock elevation map = $bedrock #"
-
+echo "# #"
+echo "########################################################"
g.remove --quiet rast=$soildepth
fi
-if [ "$GIS_FLAG_r" -eq 1 ]; then
-echo "# Soil production rate map = $soildepth"_prod_rate" #"
-echo "# #"
-echo "########################################################"
-else
-echo "# #"
-echo "########################################################"
-fi
+
More information about the grass-commit
mailing list