[GRASS-SVN] r49099 - grass-addons/raster/LandDyn/r.landscape.evol.py

svn_grass at osgeo.org svn_grass at osgeo.org
Fri Nov 4 19:01:20 EDT 2011


Author: isaacullah
Date: 2011-11-04 16:01:20 -0700 (Fri, 04 Nov 2011)
New Revision: 49099

Modified:
   grass-addons/raster/LandDyn/r.landscape.evol.py/r.landscape.evol.py
Log:
Minor update to flow depth calculation to estimate flow depth more realistically. No change to functionality.

Modified: grass-addons/raster/LandDyn/r.landscape.evol.py/r.landscape.evol.py
===================================================================
--- grass-addons/raster/LandDyn/r.landscape.evol.py/r.landscape.evol.py	2011-11-04 23:00:03 UTC (rev 49098)
+++ grass-addons/raster/LandDyn/r.landscape.evol.py/r.landscape.evol.py	2011-11-04 23:01:20 UTC (rev 49099)
@@ -3,10 +3,10 @@
 ############################################################################
 #
 # MODULE:       r.landscape.evol.py
-# AUTHOR(S):    iullah
-# COPYRIGHT:    (C) 2009 GRASS Development Team/iullah
+# AUTHOR(S):    Isaac Ullah and Michael Barton
+# COPYRIGHT:    (C) 2010 GRASS Development Team/Isaac Ullah
 #
-#  description: Create raster maps of net erosion/depostion, the modified terrain surface (DEM) after net erosion/deposition using the USPED equation, bedrock elevations after soil production, and soil depth maps. This module uses appropriate flow on different landforms by default; however, singular flow regimes can be chosen instead. THIS SCRIPT WILL PRODUCE MANY TEMPORARY MAPS AND REQUIRES A LOT OF FREE FILE SPACE! 
+#  description: Simulates the cumulative effect of erosion and deposition on a landscape over time. This module uses appropriate flow on different landforms by default; however, singular flow regimes can be chosen by manipulating the cutoff points. This module requires GRASS 6.4 or greater. THIS SCRIPT WILL PRODUCE MANY TEMPORARY MAPS AND REQUIRES A LOT OF FREE FILE SPACE! 
 
 #  This program is free software; you can redistribute it and/or modify
 #  it under the terms of the GNU General Public License as published by
@@ -20,217 +20,155 @@
 #
 #############################################################################/
 #%Module
-#% description: Create raster maps of net erosion/depostion, the modified terrain surface (DEM) after net erosion/deposition using the USPED equation, bedrock elevations after soil production, and soil depth maps. This module uses appropriate flow on different landforms by default; however, singular flow regimes can be chosen instead. THIS SCRIPT WILL PRODUCE MANY TEMPORARY MAPS AND REQUIRES A LOT OF FREE FILE SPACE! 
+#% description: Simulates the cumulative effect of erosion and deposition on a landscape over time. This module uses appropriate flow on different landforms by default; however, singular flow regimes can be chosen by manipulating the cutoff points. This module requires GRASS 6.4 or greater. THIS SCRIPT WILL PRODUCE MANY TEMPORARY MAPS AND REQUIRES A LOT OF FREE FILE SPACE! 
 #%End
 #%option
 #% key: elev
 #% type: string
 #% gisprompt: old,cell,raster
-#% description: Input elevation map (DEM)
+#% description: Input elevation map (DEM of surface)
 #% required : yes
-#% guisection: Input
 #%end
 #%option
 #% key: initbdrk
 #% type: string
 #% gisprompt: old,cell,raster
-#% description: Initial bedrock elevations map (for first iteration only)
+#% description: Bedrock elevations map (DEM of bedrock)
 #% answer: 
 #% required : yes
-#% guisection: Input
 #%end
 #%option
-#% key: prefx
+#% key: K
 #% type: string
-#% description: Prefix for all output maps
-#% answer: levol_
+#% gisprompt: old,cell,raster
+#% description: Soil erodability index (K factor) map or constant
+#% answer: 0.42
 #% required : yes
-#% guisection: Input
 #%end
 #%option
-#% key: outdem
+#% key: sdensity
 #% type: string
-#% description: Name stem for output elevation map(s) (preceded by prefix and followed by numerical suffix if more than one iteration)
-#% answer: elevation
-#% required: yes
-#% guisection: Input
+#% gisprompt: old,cell,raster
+#% description: Soil density map or constant [T/m3] for conversion from mass to volume
+#% answer: 1.2184
+#% required : yes
 #%end
 #%option
-#% key: outsoil
-#% type: string
-#% description: Name stem for the output soil depth map(s) (preceded by prefix and followed by numerical suffix if more than one iteration)
-#% answer: soildepth
-#% required: yes
-#% guisection: Input
+#% key: Kt
+#% type: double
+#% description: Stream transport efficiency variable (0.001 for a soft substrate, 0.0001 for a normal substrate, 0.00001 for a hard substrate, 0.000001 for a very hard substrate)
+#% answer: 0.0001
+#% options: 0.001,0.0001,0.00001,0.000001
+#% required : yes
 #%end
 #%option
-#% key: outbdrk
-#% type: string
-#% description: Name stem for the output bedrock map(s) (required if the -b option is NOT checked; preceded by prefix and followed by numerical suffix if more than one iteration)
-#% answer: bedrock
-#% required: no
-#% guisection: Input
+#% key: loadexp
+#% type: double
+#% description: Stream transport type variable (1.5 for mainly bedload transport, 2.5 for mainly suspended load transport)
+#% answer: 1.5
+#% options: 1.5,2.5
+#% required : yes
 #%end
-#%Option
-#% key: statsout
-#% type: string
-#% description: Name for the statsout text file (optional, if none provided, a default name will be used)
-#% required: no
-#% guisection: Input
-#%end
-#%flag
-#% key: g
-#% description: -g Don't put header on statsout text file and always append data, even if stats file already exists (useful if script is being run by an outside program)
-#% guisection: Input
-#%end
-#%flag
-#% key: l
-#% description: -l Don't output yearly soil depth maps
-#% guisection: Input
-#%end
-#%flag
-#% key: k
-#% description: -k Keep all temporary maps (do not clean up)
-#% guisection: Input
-#%end
-#%flag
-#% key: e
-#% description: -e Keep initial soil depths map (for year 0)
-#% guisection: Input
-#%end
-#%flag
-#% key: b
-#% description: -b Enable experimental bedrock weathering (soil production)
-#% guisection: Input
-#%end
-#%flag
-#% key: p
-#% description: -p Keep all slope maps 
-#% guisection: Input
-#%end
-#%flag
-#% key: n
-#% description: -n Don't output yearly net elevation change (netchange) maps
-#% guisection: Input
-#%end
-
 #%option
-#% key: R
-#% type: string
-#% description: Rainfall (R factor) constant (AVERAGE FOR WHOLE MAP AREA)
-#% answer: 5.66
+#% key: kappa
+#% type: double
+#% description: Hillslope diffusion (Kappa) rate map or constant [m/kyr]
+#% answer: 1
 #% required : yes
-#% guisection: Variables
 #%end
 #%option
-#% key: K
+#% key: C
 #% type: string
 #% gisprompt: old,cell,raster
-#% description: Soil erodability index (K factor) map or constant
-#% answer: 0.42
+#% description: Landcover index (C factor) map or constant
+#% answer: 0.005
 #% required : yes
-#% guisection: Variables
 #%end
 #%option
-#% key: sdensity
+#% key: rain
 #% type: string
 #% gisprompt: old,cell,raster
-#% description: Soil density constant (for conversion from mass to volume)
-#% answer: 1.2184
+#% description: Precip totals for the average storm [mm]
+#% answer: 10
 #% required : yes
-#% guisection: Variables
 #%end
 #%option
-#% key: C
+#% key: R
 #% type: string
-#% gisprompt: old,cell,raster
-#% description: Landcover index (C factor) map or constant
-#% answer: 0.001
+#% description: Rainfall (R factor) constant (AVERAGE FOR WHOLE MAP AREA)
+#% answer: 5.66
 #% required : yes
-#% guisection: Variables
 #%end
 #%option
-#% key: rain
-#% type: string
-#% gisprompt: old,cell,raster
-#% description: Anuall precip totals map or constant (in meters per year)
-#% answer: 0.2
+#% key: storms
+#% type: integer
+#% description: Number of storms per year (integer)
+#% answer: 75
 #% required : yes
-#% guisection: Variables
 #%end
 #%option
-#% key: raindays
-#% type: string
-#% description: Number of days of rain per year (integer)
-#% answer: 100
+#% key: stormlength
+#% type: double
+#% description: Average length of the storm [h]
+#% answer: 6.0
 #% required : yes
-#% guisection: Variables
 #%end
 #%option
-#% key: infilt
-#% type: string
-#% description: Proportion of rain that infiltrates into the ground rather than runs off (0.0 to 1.0)
-#% answer: 0.1
+#% key: speed
+#% type: double
+#% description: Average velocity of flowing water in the drainage [m/s]
+#% answer: 1.4
 #% required : yes
-#% guisection: Variables
 #%end
 #%option
-#% key: Kt
-#% type: string
-#% description: Stream transport efficiency variable (0.001 for a normal substrate (ie. sediment) to 0.000001 for a very hard substrate like bedrock)
-#% answer: 0.001
-#% options: 0.001,0.0001,0.00001,0.00001,0.000001
+#% key: cutoff1
+#% type: double
+#% description: Flow depth [m] breakpoint value for shift from diffusion to overland flow
+#% answer: 10
 #% required : yes
-#% guisection: Variables
 #%end
 #%option
-#% key: loadexp
-#% type: string
-#% description: Stream transport type varaible (1.5 for mainly bedload transport, 2.5 for mainly suspended load transport)
-#% answer: 1.5
-#% options: 1.5,2.5
+#% key: cutoff2
+#% type: double
+#% description: Flow depth [m] breakpoint value for shift from overland flow to rill/gully flow (if value is the same as cutoff1, no sheetwash processes will be modeled)
+#% answer: 80
 #% required : yes
-#% guisection: Variables
 #%end
 #%option
-#% key: kappa
-#% type: string
-#% description: Hillslope diffusion (Kappa) rate map or constant (meters per kiloyear)
-#% answer: 1
+#% key: cutoff3
+#% type: double
+#% description: Flow depth [m] breakpoint value for shift from rill/gully flow to stream flow (if value is the same as cutoff2, no rill processes will be modeled)
+#% answer: 160
 #% required : yes
-#% guisection: Variables
 #%end
 #%option
-#% key: alpha
+#% key: smoothing
 #% type: string
-#% description: Critical slope threshold for mass movement of sediment (in degrees above horizontal)
-#% answer: 40
+#% description: Amount of modal-smoothing (low smoothing is recommended)
+#% answer: low
+#% options: no,low,high
 #% required : yes
-#% guisection: Variables
 #%end
 #%option
-#% key: cutoff1
+#% key: prefx
 #% type: string
-#% description: Flow accumultion breakpoint value for shift from diffusion to overland flow (sq. m uplsope area)
-#% answer: 900
+#% description: Prefix for all output maps
+#% answer: levol_
 #% required : yes
-#% guisection: Variables
 #%end
 #%option
-#% key: cutoff2
+#% key: outdem
 #% type: string
-#% description: Flow accumultion breakpoint value for shift from overland flow to rill/gully flow (sq. m uplsope area)
-#% answer: 11250
-#% required : yes
-#% guisection: Variables
+#% description: Name stem for output elevation map(s) (preceded by prefix and followed by numerical suffix if more than one iteration)
+#% answer: elevation
+#% required: yes
 #%end
 #%option
-#% key: cutoff3
+#% key: outsoil
 #% type: string
-#% description: Flow accumultion breakpoint value for shift from rill/gully flow to stream flow (sq. m uplsope area)
-#% answer: 225000
-#% required : yes
-#% guisection: Variables
+#% description: Name stem for the output soil depth map(s) (preceded by prefix and followed by numerical suffix if more than one iteration)
+#% answer: soildepth
+#% required: yes
 #%end
 #%option
 #% key: number
@@ -238,117 +176,87 @@
 #% description: Number of iterations (cycles) to run
 #% answer: 1
 #% required : yes
-#% guisection: Variables
 #%end
 
 
-#these are commented out as we currently utilize the profile curvature method described by Heimsath et al...
-# 	#%option
-# 	#% key: Ba
-# 	#% type: string
-# 	#% description: Rate of average soil production (Ba)
-# 	#% answer: 0.00008
-# 	#% required : yes
-# 	#% guisection: Variables
-# 	#%end
-# 	#%option
-# 	#% key: Bb
-# 	#% type: string
-# 	#% description: Relationship between soil depth and production rate (Bb)
-# 	#% answer: 0.1
-# 	#% required : yes
-# 	#% guisection: Variables
-# 	#%end
+# #%option
+# #% key: alpha                 ###This may be added back in if I can find a good equation for mass movement
+# #% type: integer
+# #% description: Critical slope threshold for mass movement of sediment (in degrees above horizontal)
+# #% answer: 40
+# #% required : yes
+# #%end
 
 #%flag
-#% key: w
-#% description: -w Calcuate for only sheetwash across entire map
-#% guisection: Flow_type
+#% key: p
+#% description: -p Output a vector points map with sampled values of flow accumulation and curvatures suitable for determining cutoff values. NOTE: Overrides all other output options, and exits after completion. The output vector points map will be named  "PREFIX_#_randomly_sampled_points".
+#% guisection: Optional
 #%end
 #%flag
-#% key: r
-#% description: -r Calcuate for only channelized flow across entire map
-#% guisection: Flow_type
+#% key: 1
+#% description: -1 Calculate streams as 1D difference instead of 2D divergence (forces heavy incision in streams)
+#% guisection: Optional
 #%end
 #%flag
+#% key: k
+#% description: -k Keep ALL temporary maps (overides flags -drst). This will make A LOT of maps!
+#% guisection: Optional
+#%end
+#%flag
 #% key: d
-#% description: -d Calcuate for only diffusive flow (soil creep) across entire map
-#% guisection: Flow_type
+#% description: -d Don't output yearly soil depth maps
+#% guisection: Optional
 #%end
 #%flag
-#% key: f
-#% description: -f Use r.terrflow instead of r.watershed to calculate flow accumulation ( GRASS 6.3.x users MUST use this flag!)
-#% answer: 0
-#% guisection: Flow_type
+#% key: r
+#% description: -r Don't output yearly maps of the erosion/deposition rates ("ED_rate" map, in vertical meters)
+#% guisection: Optional
 #%end
+#%flag
+#% key: s
+#% description: -s Keep all slope maps 
+#% guisection: Optional
+#%end
+#%flag
+#% key: t
+#% description: -t Keep yearly maps of the Transport Capacity at each cell ("Qs" maps)
+#% guisection: Optional
+#%end
+#%flag
+#% key: e
+#% description: -e Keep yearly maps of the Excess Transport Capacity (divergence) at each cell ("DeltaQs" maps)
+#% guisection: Optional
+#%end
+#%Option
+#% key: statsout
+#% type: string
+#% description: Name for the statsout text file (optional, if none provided, a default name will be used)
+#% required: no
+#% guisection: Optional
+#%end
 
-
 import sys
 import os
 import subprocess
+import math
+import random
 import tempfile
 grass_install_tree = os.getenv('GISBASE')
 sys.path.append(grass_install_tree + os.sep + 'etc' + os.sep + 'python')
 import grass.script as grass
 
-# first define some useful custom methods
-
-# m is a message (as a string) one wishes to have printed in the output window
-def grass_message(m):
-	subprocess.Popen('g.message message="%s"' % m, shell='bash').wait()
-	return
-
-# m is grass (or bash) command to execute (written as a string). script will wait for grass command to finish
-def grass_com(m):
-	subprocess.Popen('%s' % m, shell='bash').wait()
-	return
-
-# m is grass (or bash) command to execute (written as a string). script will not wait for grass command to finish
-def grass_com_nowait(m):
-	subprocess.Popen('%s' % m, shell='bash')
-	return
-
-# m is a grass/bash command that will generate some list of keyed info to stdout, n is the character that seperates the key from the data, o is a defined blank dictionary to write results to
-def out2dict(m, n, o):
-    p1 = subprocess.Popen('%s' % m, stdout=subprocess.PIPE, shell='bash')
-    p2 = p1.stdout.readlines()
-    for y in p2:
-        y0,y1 = y.split('%s' % n)
-        o[y0] = y1.strip('\n')
-
-# m is a grass/bash command that will generate some list of info to stdout seperated by newlines, n is a defined blank list to write results to
-def out2list(m, n):
-        p1 = subprocess.Popen('%s' % m, stdout=subprocess.PIPE, shell='bash')
-        p2 = p1.stdout.readlines()
-        for y in p2:
-            n.append(y.strip('\n'))
-#SAVE!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+# First define two useful custom methods that we will need later on
 # m is a grass/bash command that will generate some info to stdout. You must invoke this command in the form of "variable to be made" = out2var('command')
 def out2var(m):
         pn = subprocess.Popen('%s' % m, stdout=subprocess.PIPE, shell='bash')
         return pn.stdout.read()
-
-# returns amoutn of free ramused like: var = getram()
-def getram():
-    if sys.platform == "win32":
-        #mem1 = os.popen('mem | find "total"').readlines()
-        mem = subprocess.Popen('mem | find "total"',  stdout=subprocess.PIPE)
-        mem1 = mem.stdout.readlines()
-        FreeMemory = int(mem1[0].split()[0])
-    else:
-        mem = subprocess.Popen('free -mo  | grep "Mem" | tr -s /[:blank:] /[:] | cut -d":" -f4')
-        mem1 = mem.stdout.read()
-        FreeMemory = int(mem1.strip('\n'))
-    return FreeMemory
     
-# now define  "main",  our main block of code, here defined because of the way g.parser needs to be called with python codes for grass (see below)        
-
+# Now define  "main",  our main block of code, here defined because of the way g.parser needs to be called with python codes for grass (see below)        
 # m = last iteration number, o = iteration number, p = prefx, q = statsout, r = resolution of input elev map
 def main(m, o, p, q, r):
-    # get variables from user input
+    # Get variables from user input
     years = os.getenv("GIS_OPT_number")
     initbdrk = os.getenv("GIS_OPT_initbdrk")
-    outbdrk = os.getenv("GIS_OPT_outbdrk")
     outdem = os.getenv("GIS_OPT_outdem")
     outsoil = os.getenv("GIS_OPT_outsoil")
     R = os.getenv("GIS_OPT_R")
@@ -356,238 +264,359 @@
     sdensity = os.getenv("GIS_OPT_sdensity")
     C = os.getenv("GIS_OPT_C")
     kappa = os.getenv("GIS_OPT_kappa")
-    aplpha = os.getenv("GIS_OPT_method")
-    #Ba = os.getenv("GIS_OPT_Ba")
-    #Bb = os.getenv("GIS_OPT_Bb")
-    sigma = os.getenv("GIS_OPT_sigma")
-    nbhood = os.getenv("GIS_OPT_nbhood")
-    cutoff1 = str(int(os.getenv("GIS_OPT_cutoff1")) / int(float(r) * float(r)))
-    cutoff2 = str(int(os.getenv("GIS_OPT_cutoff2")) / int(float(r) * float(r)))
-    cutoff3 = str(int(os.getenv("GIS_OPT_cutoff3")) / int(float(r) * float(r)))
+    cutoff1 = os.getenv("GIS_OPT_cutoff1")
+    cutoff2 = os.getenv("GIS_OPT_cutoff2")
+    cutoff3 = os.getenv("GIS_OPT_cutoff3")
     rain = os.getenv("GIS_OPT_rain")
-    raindays = os.getenv("GIS_OPT_raindays")
-    infilt = os.getenv("GIS_OPT_infilt")
+    storms = os.getenv("GIS_OPT_storms")
+    stormlengthsecs = float(os.getenv("GIS_OPT_stormlength"))*3600.00                  # number of seconds in the storm
+    stormtimet = stormlengthsecs / (float(os.getenv("GIS_OPT_speed")) * float(r))       # number of hydrologic instants in the storm
+    timet = stormlengthsecs/stormtimet                                                                      # length of a single hydrologic instant in seconds
     Kt = os.getenv("GIS_OPT_Kt")
     loadexp = os.getenv("GIS_OPT_loadexp")
-    method = os.getenv("GIS_OPT_method")
-    # make some variables for temporary map names
-    flowacc = '%sflowacc%s' % (p, o)
-    aspect = '%saspect%s' % (p, o)
-    pc = '%spc%s' % (p, o)
-    tc = '%stc%s' % (p, o)
-    meancurv = '%smeancurv%s' % (p, o)
-    rate = '%srate%s' % (p, o)
-    sflowtopo = '%ssflowtopo%s' % (p, o)
-    qsx = '%sqsx%s' % (p, o)
-    qsy = '%sqsy%s' % (p, o)
-    qsxdx = '%sqsx_dx%s' % (p, o)
-    qsydy = '%sqsy_dy%s' % (p, o)
-    erdep = '%serosdep%s' % (p, o)
-    # make color rules for netchange maps
+    # Make some variables for temporary map names
+    aspect = '%saspect%04d' % (p, o)
+    flowacc = '%sflowacc%04d' % (p, o)
+    flacclargenums = flowacc + '_largenums'
+    flowdir = '%sflowdir%04d' % (p, o)
+    pc = '%spc%04d' % (p, o)
+    tc = '%stc%04d' % (p, o)
+    meancurv = '%smeancurv%04d' % (p, o)
+    rate = '%srate%04d' % (p, o)
+    # Make color rules for netchange maps
     nccolors = tempfile.NamedTemporaryFile()
     nccolors.write('100% 0 0 100\n1 blue\n0.5 indigo\n0.01 green\n0 white\n-0.01 yellow\n-0.5 orange\n-1 red\n0% 150 0 50')
     nccolors.flush()
-    # make color rules for soil depth maps
+    # Make color rules for soil depth maps
     sdcolors = tempfile.NamedTemporaryFile()
     sdcolors.write('100% 0:249:47\n20% 78:151:211\n6% 194:84:171\n0% 227:174:217')
     sdcolors.flush()
-    # if first iteration, use input maps. Otherwise, use maps generated from previous iterations
+    # If first iteration, use input maps. Otherwise, use maps generated from previous iterations
     if ( o == 1 ):
         old_dem = '%s' % os.getenv("GIS_OPT_elev")
         old_bdrk = '%s' % os.getenv("GIS_OPT_initbdrk")
         old_soil = "%s%s_init" % (prefx, os.getenv("GIS_OPT_outsoil"))
-        grass.mapcalc('${out}=${m1}-${m2}', out = old_soil, m1 = old_dem, m2 = old_bdrk)
+        grass.mapcalc('${old_soil}=${old_dem}-${old_bdrk}', old_soil = old_soil, old_dem = old_dem, old_bdrk = old_bdrk)
     else :
-        old_dem = '%s%s%s' % (p, os.getenv("GIS_OPT_outdem"), m)
-        old_bdrk = '%s%s%s' % (p, os.getenv("GIS_OPT_outbdrk"), m)
-        old_soil = '%s%s%s' % (p, os.getenv("GIS_OPT_outsoil"), m)
-    #checking for special condition of there being only one run, and setting variables accordingly (one year runs have no numbers suffixed to the output map names)
+        old_dem = '%s%s%04d' % (p, os.getenv("GIS_OPT_outdem"), m)
+        old_bdrk = '%s%s%04d' % (p, os.getenv("GIS_OPT_outbdrk"), m)
+        old_soil = '%s%s%04d' % (p, os.getenv("GIS_OPT_outsoil"), m)
+    #Checking for special condition of there being only one run, and setting variables accordingly (one year runs have no numbers suffixed to the output map names)
     if ( years == '1' ):
         slope = '%sslope' % p
-        netchange = '%snetchange' % p
+        qs = '%sQs' % (p)
+        netchange = '%sED_rate' % p
         new_dem ='%s%s' % (p, outdem)
-        new_bdrk = '%s%s' % (p, outbdrk)
         new_soil = '%s%s' % (p, outsoil)
     else:
-        slope = '%sslope%s' % (p, o)
-        netchange = '%snetchange%s' % (p, o)
-        new_dem = '%s%s%s' % (p, outdem, o)
-        new_bdrk = '%s%s%s' % (p, outbdrk, o)
-        new_soil = '%s%s%s' % (p, outsoil, o)
-    grass.message('\n##################################################\n\n*************************\n Year %s ' % o + 'step 1 of 7: calculating slope and aspect\n*************************\n')
-    grass.run_command('r.slope.aspect',  quiet = True, elevation = old_dem, slope = slope, aspect = aspect,  pcurv = pc, tcurv = tc)
-    grass.message('\n*************************\n Year %s ' % o + 'step 2 of 7: calculating upslope accumulated flow\n*************************\n')       
-    if ( os.getenv("GIS_FLAG_f") == "1" ):
-        grass.message('using r.terraflow to calculate overland flow accumulation per cell (number of cells uplsope from each cell)\n\n')   
-        #First we need to grab the amount of free RAM for r.terraflow
-        mem = getram()
-        #r.terraflow can't handle it if you tell it to use more than 2 Gigs of RAM, so if you have more than that, we have to tell r.terraflow to only use up to 2 Gigs of the free RAM... 
-        if int(mem) > 2000:
-            mem = '2000'    
-        grass.message('Amount of free RAM being allocated for this step: %s Megabites' % mem)
-        grass.run_command('r.terraflow',  quiet = True,  elev = old_dem, filled = p + '.filled', direction = p + '.direct', swatershed = p + '.sink', accumulation = flowacc, tci = p + '.tci',  d8cut = 'infinity', memory = mem, STREAM_DIR = os.path.expanduser('~'))
-        grass.run_command('g.remove', quiet =True, rast = p + '.filled, ' + p + '.direct, ' + p + '.sink, ' + p + '.tci')
-        os.remove(os.path.expanduser('~') + os.sep + 'STREAM*')
+        slope = '%sslope%04d' % (p, o)
+        netchange = '%sED_rate%04d' % (p, o)
+        new_dem = '%s%s%04d' % (p, outdem, o)
+        new_soil = '%s%s%04d' % (p, outsoil, o)
+    if ( os.getenv("GIS_FLAG_p") == "1" ):
+        grass.message('1) Calculating slope and curvatures')
+        grass.run_command('r.slope.aspect',  quiet = True, elevation = old_dem, slope = slope, pcurv = pc, tcurv = tc)
     else:
-        grass.message('Using r.watershed to calculate overland flow accumulation per cell (number of cells uplsope from each cell)')
-        if '<flag name="f">' in out2var('r.watershed --interface-description'):
-            grass.run_command('r.watershed',  quiet = True,  flags = 'fa',  elevation = old_dem, accumulation = flowacc,  convergence = '5')
+        grass.message('\n##################################################\n\n*************************\n Year %s -- ' % o + 'step 1: calculating slope\n*************************\n')
+        grass.run_command('r.slope.aspect',  quiet = True, elevation = old_dem, aspect = aspect, slope = slope)
+    if ( os.getenv("GIS_FLAG_p") == "1" ):
+        grass.message('2) Calculating map of rainfall excess')
+    else:
+        grass.message('\n*************************\n Year %s -- ' % o + 'step 2: calculating accumulated flow depths\n*************************\n')       
+        grass.message('Calculating rainfall excess rates from the input rainfall and landcover data. (units are in 100,000ths of a meter')
+    rainexcess = "%s_rainfall_excess_map_%04d"% (p, o)
+    grass.mapcalc('${rainexcess}=int((${rain}) * (252.288*exp(${C}, 0.401896)))', rainexcess = rainexcess, rain = rain, C = C)
+    if ( os.getenv("GIS_FLAG_p") == "1" ):
+        grass.message('3) Calculating accumulated flow (in meters of accumulated depth')
+    grass.message('Calculating overland flow accumulation per cell (total vertical depth of water that passes over each cell in one storm)')
+    try:
+        grass.run_command('r.watershed',  quiet = True,  flags = 'fa',  elevation = old_dem, flow = rainexcess, accumulation = flacclargenums, drainage = flowdir, convergence = '5')
+    except:
+        grass.run_command('r.watershed',  quiet = True,  flags = 'a',  elevation = old_dem, flow = rainexcess, accumulation = flacclargenums, drainage = flowdir, convergence = '5')
+#    if '<flag name="f">' in out2var('r.watershed --interface-description'):
+#        grass.run_command('r.watershed',  quiet = True,  flags = 'fa',  elevation = old_dem, flow = runoffp, accumulation = flacclargenums, drainage = flowdir, convergence = '5')
+#    else:
+#        grass.run_command('r.watershed',  quiet = True,  flags = 'a',  elevation = old_dem, flow = runoffp, accumulation = flacclargenums, drainage = flowdir, convergence = '5')
+    grass.mapcalc('${flowacc}=${flacclargenums}/100000', flowacc = flowacc, flacclargenums = flacclargenums)
+    if ( os.getenv("GIS_FLAG_k") == "1" ):
+        pass
+    else:
+        grass.run_command('g.remove', quiet = True, rast = rainexcess + "," + flacclargenums)
+    if ( os.getenv("GIS_FLAG_p") == "1" ):
+        grass.message('4) Determining number of sampling points using formula: "ln(#cells_in_input_map)*100"')
+        flaccstats = grass.parse_command('r.univar', flags = 'g', map = flowacc)
+        numpts = int(math.log(int(flaccstats['n']))*100)
+        grass.message('5) Creating random points and sampling values of flow accumulation, curvatures, and slope.')
+        vout = '%s%s_randomly_sampled_points' % (p, numpts)
+        grass.run_command('r.random', quiet = True, input = flowacc, cover = pc, n = numpts, vector_output = vout)
+        grass.run_command('v.db.renamecol', quiet = True, map = vout, column = 'value,Flow_acc')
+        grass.run_command('v.db.renamecol', quiet = True, map = vout, column = 'covervalue,Princ_curv')
+        grass.run_command('v.db.addcol', quiet = True, map = vout, columns = 'Tang_curv double precision, Slope double precision')
+        grass.run_command('v.what.rast', quiet = True, vector = vout, raster = tc, column = "Tang_curv")
+        grass.run_command('v.what.rast', quiet = True, vector = vout, raster = slope, column = "Slope")
+        if ( os.getenv("GIS_FLAG_k") == "1" ):
+            grass.message('--Keeping the created maps (Flow Accumulation, Slope, Principle Curvature, Tangential Curvature)')
         else:
-            grass.run_command('r.watershed',  quiet = True,  flags = 'a',  elevation = old_dem, accumulation = flowacc,  convergence = '5')
+            grass.message('6) Cleaning up...')
+            grass.run_command('g.remove', quiet = True,  rast = slope + "," + pc + "," + tc + "," + flowacc)
+        grass.run_command('r.mask', flags = 'r')
+        if ismask == 1:
+            grass.run_command('g.rename', quiet = "True", rast = tempmask +',MASK')
+        grass.message('FINISHED. \nRandom sample points map "%s" created successfully.\n' % vout)
+        sys.exit(0)
 
-    grass.message('\n*************************\n Year %s ' % o + 'step 3 of 7: calculating basic sediment transport rates\n*************************\n')  
-    error_message = '\n############################################################\n !!!!!!!!!YOU MUST SELECT ONLY ONE TYPE OF EROSION!!!!!!!!!!!\n ############################################################ \n \n'
-    if ( os.getenv("GIS_FLAG_w") == "1" and  os.getenv("GIS_FLAG_d") == "0" and os.getenv("GIS_FLAG_r") == "0" ):
-        grass.message('calculating for sheetwash across the entire map')
-        grass.mapcalc('${sflowtopo} = ${flowacc} * sin(${slope})',  sflowtopo = sflowtopo,  flowacc = flowacc,  slope = slope)
-    elif ( os.getenv("GIS_FLAG_w") == "0" and  os.getenv("GIS_FLAG_d") == "0" and os.getenv("GIS_FLAG_r") == "1" ):
-        grass.message('calculating for channelzed flow across the entire map')
-        grass.mapcalc('${sflowtopo} = exp(${flowacc}, 1.6) * exp(sin(${slope}), 1.3)',  sflowtopo = sflowtopo,  flowacc = flowacc,  slope = slope)
-    elif ( os.getenv("GIS_FLAG_w") == "0" and  os.getenv("GIS_FLAG_d") == "1" and os.getenv("GIS_FLAG_r") == "0" ):
-        grass.message('calculating for diffusive flow across the entire map')
-        grass.mapcalc('${sflowtopo} = ${kappa} * exp((${r},2) * sin(${slope})',  out = sflowtopo,  kappa = kappa,  r = r,  slope = slope)
-    elif ( os.getenv("GIS_FLAG_w") == "0" and  os.getenv("GIS_FLAG_d") == "1" and os.getenv("GIS_FLAG_r") == "1" ):
-        sys.exit(error_message)
-    elif ( os.getenv("GIS_FLAG_w") == "1" and  os.getenv("GIS_FLAG_d") == "0" and os.getenv("GIS_FLAG_r") == "1" ):
-        sys.exit(error_message)
-    elif ( os.getenv("GIS_FLAG_w") == "1" and  os.getenv("GIS_FLAG_d") == "1" and os.getenv("GIS_FLAG_r") == "0" ):
-        sys.exit(error_message)
-    elif ( os.getenv("GIS_FLAG_w") == "1" and  os.getenv("GIS_FLAG_d") == "1" and os.getenv("GIS_FLAG_r") == "1" ):
-        sys.exit(error_message)
+    grass.message('\n*************************\n Year %s -- ' % o + 'step 3: calculating sediment transport rates (units variable depending upon process) \n*************************\n')  
+    # This step calculates the force of the flowing water at every cell on the landscape using the proper transport process law for the specific point in the flow regime. For upper hillslopes (below cutoff point 1) this done by multiplying the diffusion coeficient by the accumulated flow/cell res width. For midslopes (between cutoff 1 and 2) this is done by multiplying slope by accumulated flow with the m and n exponents set to 1. For channel catchment heads (between cutoff 2 and 3), this is done by multiplying slope by accumulated flow with the m and n exponents set to 1.6 and 1.3 respectively. For Channelized flow in streams (above cutoff 3), this is done by calculating the reach average shear stress (hydraulic radius [here estimated for a cellular landscape simply as the depth of flow]  times  slope times accumulated flow [cells] times gravitatiopnal acceleration of water [9806.65 newtons], all raised to the appropriate exponant for the type of transport (bedload or suspe
 nded load), and then divided by the resolution. Depth of flow is calculated as a mean "instantaneous depth" during any given rain event, here estimated by the maximum depth of an idealized unit hydrograph with base equal to the duration of the storm, and area equal to the total accumulated excess rainfall during the storm. Then finally calculates the stream power or sediment carrying capacity (qs) of the water flowing at each part of the map by multiplying the reach average shear stress (channelized flow in streams) or the estimated flow force (overland flow) by the transport coeficient (estimated by R*K*C for hillslopes or kt for streams). This is a "Transport Limited" equation, however, we add some constraints on detachment by checking to see if the sediment supply has been exhausted: if the current soil depth is 0 or negative (checking for a negative value is kind of an error trap) then we make the transport coefficient small (0.000001) to simulate erosion on bedrock. Bec
 ause diffusion and USPED require 2D divergence later on, we calculate these as vectors in the X and Y directions. Stream flow only needs 1D difference, so it's calulated in the direction of flow.
+    
+    qsx4 = '%sQsx_diffusion%04d' % (p, o)
+    qsy4 = '%sQsy_diffusion%04d' % (p, o)
+    grass.mapcalc('${qsx4}=(${kappa} * sin(${slope}) * cos(${aspect}))', qsx4 = qsx4, kappa = kappa, slope = slope, aspect = aspect)
+    grass.mapcalc('${qsy4}=(${kappa} * sin(${slope}) * sin(${aspect}))', qsy4 = qsy4, kappa = kappa, slope = slope, aspect = aspect)
+    if cutoff1 == cutoff2:
+        qsx3 = qsx4
+        qsy3 = qsy4
     else:
-        # This step calculates the force of the flowing water at every cell on the landscape using the proper transport process law for the specific point in the flow regime. For upper hillslopes (below cutoff point 1) this done by multiplying the diffusion coeficient by the accumulated flow/cell res width. For midslopes (between cutoff 1 and 2) this is done by multiplying slope by accumulated flow with the m and n exponents set to 1. For channel catchment heads (between cutoff 2 and 3), this is done by multiplying slope by accumulated flow with the m and n exponents set to 1.6 and 1.3 respectively. For Channelized flow in streams (above cutoff 3), this is done by calculating the reach average shear stress (hydraulic radius [here estimated for a cellular landscape simply as the depth of flow]  times  slope times accumulated flow [cells] times gravitatiopnal acceleration of water [9806.65 newtons], all raised to the appropriate exponant for the type of transport (bedload or s
 uspended load]). Depth of flow is claculated as a mean "instantaneous depth" during any given rain event, here assumed to be roughly equivelent tothe average depth of flow during 1 minute. 
-        grass.mapcalc('${sflowtopo}= if(${flowacc} >= ${cutoff3}, exp((9806.65 * ((${rain} - (${rain} * ${infilt})) / (${raindays} * 1440)) * ${flowacc} * ${slope}), ${loadexp}), if(${flowacc} >= ${cutoff2} && ${flowacc} < ${cutoff3}, (exp((${flowacc}*${r}),1.6000000) * exp(sin(${slope}),1.3000000)), if(${flowacc} >= ${cutoff1} && ${flowacc} < ${cutoff2}, ((${flowacc}*${r}) * sin(${slope})), (${kappa} * sin(${slope}) ))))', sflowtopo = sflowtopo, rain = rain, infilt = infilt, raindays = raindays, flowacc = flowacc, cutoff1 = cutoff1,  cutoff2 = cutoff2,  cutoff3 = cutoff3,  slope = slope, kappa = kappa,  loadexp = loadexp,  r = r)
-    
-    grass.message('\n*************************\n Year %s ' % o + 'step 4 of 7: calculating sediment transport capacity in x and y directions\n*************************\n\n')
-        #This step calculates the stream power or sediment carrying capacity (qs) of the water flowing at each part of the map by multiplying the reach average shear stress (channelized flow in streams) or the estimated flow force (overland flow) by the transport coeficient (estimated by R*K*C for hillslopes or kt for streams). This is a "Transport Limited" equation, however, we add some constraints on detachment by checking to see if the sediment supply has been exhausted: if the current soil depth is 0 or negative (checking for a negative value is kind of an error trap) then we make the transport coefficient small (0.000001) to simulate erosion on bedrock.
-    grass.mapcalc('${qsx}=if(${old_soil} <= 0, (0.000001 * ${sflowtopo} * cos(${aspect})), if(${old_soil} > 0 && ${flowacc} >= ${cutoff3}, (${Kt} * ${sflowtopo} * cos(${aspect})), (${R} * ${K} * ${C} * ${sflowtopo} * cos(${aspect}))))', qsx = qsx, old_soil = old_soil, flowacc = flowacc, sflowtopo = sflowtopo, aspect = aspect, cutoff3 = cutoff3, Kt = Kt, R = R, K = K, C = C)
-    grass.mapcalc('${qsy}=if(${old_soil} <= 0, (0.000001 * ${sflowtopo} * sin(${aspect})), if(${old_soil} > 0 && ${flowacc} >= ${cutoff3}, (${Kt} * ${sflowtopo} * sin(${aspect})), (${R} * ${K} * ${C} * ${sflowtopo} * sin(${aspect}))))', qsy = qsy, old_soil = old_soil, flowacc = flowacc, sflowtopo = sflowtopo, aspect = aspect, cutoff3 = cutoff3, Kt = Kt, R = R, K = K, C = C)
+        qsx3 = '%sQsx_sheetwash%04d' % (p, o)
+        qsy3 = '%sQsy_sheetwash%04d' % (p, o)
+        grass.mapcalc('${qsx3}=( (${R}*${K}*${C}*${flowacc}*${res}*sin(${slope})) * cos(${aspect}) )', qsx3 = qsx3, R = R, K = K, C =C, slope = slope, res = r, flowacc = flowacc, aspect = aspect,  sdensity = sdensity)
+        grass.mapcalc('${qsy3}=( (${R}*${K}*${C}*${flowacc}*${res}*sin(${slope})) * sin(${aspect}) )', qsy3 = qsy3, R = R, K = K, C =C, slope = slope, res = r, flowacc = flowacc, aspect = aspect,  sdensity = sdensity)
+    if cutoff2 == cutoff3:
+        qsx2 = qsx3
+        qsy2 = qsy3
+    else:
+        qsx2 = '%sQsx_rills%04d' % (p, o)
+        qsy2 = '%sQsy_rills%04d' % (p, o)
+        grass.mapcalc('${qsx2}=( (${R}*${K}*${C}*exp((${flowacc}*${res}),1.6000000)*exp(sin(${slope}),1.3000000)) * cos(${aspect}))', qsx2 = qsx2, R = R, K = K, C =C,  slope = slope, res = r, flowacc = flowacc, aspect = aspect,  sdensity = sdensity)
+        grass.mapcalc('${qsy2}=( (${R}*${K}*${C}*exp((${flowacc}*${res}),1.6000000)*exp(sin(${slope}),1.3000000)) * sin(${aspect}))', qsy2 = qsy2, R = R, K = K, C =C,  slope = slope, res = r, flowacc = flowacc, aspect = aspect,  sdensity = sdensity)
+    if ( os.getenv("GIS_FLAG_1") == "1" ):
+        #This is the 1D version
+        qs1 = '%sQs_1D_streams%04d' % (p, o)
+        grass.mapcalc('${qs1}=(${Kt} * exp(9806.65*((2*${flowacc})/(${stormtimet}+(0.25*${stormtimet})))*sin(${slope}), ${loadexp}) )',  qs1 = qs1, flowacc = flowacc, stormtimet = stormtimet, slope = slope, loadexp = loadexp,  Kt = Kt,  sdensity = sdensity) 
+    else:
+        
+        #This is the 2D version
+        qsx1 = '%sQsx_streams%04d' % (p, o)
+        qsy1 = '%sQsy_streams%04d' % (p, o)
+        grass.mapcalc('${qsx1}=(${Kt} * exp(9806.65*(${flowacc}/(0.595*${stormtimet}))*sin(${slope}), ${loadexp}) ) * cos(${aspect})',  qsx1 = qsx1, flowacc = flowacc, stormtimet = stormtimet, slope = slope, loadexp = loadexp,  Kt = Kt,  sdensity = sdensity, aspect = aspect) 
+        grass.mapcalc('${qsy1}=(${Kt} * exp(9806.65*(${flowacc}/(0.595*${stormtimet}))*sin(${slope}), ${loadexp}) ) * sin(${aspect})',  qsy1 = qsy1, flowacc = flowacc, stormtimet = stormtimet, slope = slope, loadexp = loadexp,  Kt = Kt,  sdensity = sdensity, aspect = aspect)
 
-
-    grass.message('\n*************************\n Year %s ' % o + 'step 5 of 7: calculating partial derivatives for sediment transport\n*************************\n\n')
-    grass.run_command('r.slope.aspect',  quiet = True,  elevation = qsx, dx = qsxdx)
-    grass.run_command('r.slope.aspect',  quiet = True,  elevation = qsy, dy = qsydy)
-
-    grass.message('\n*************************\n Year %s ' % o + 'step 6 of 7: calculating net erosion and deposition in meters of vertical change per cell\n*************************\n\n')
-    #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 th
 e 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).
-    grass.mapcalc('${erdep}= if( ${old_soil} >= 0, if( (-1 * ( ( (${qsxdx} + ${qsydy}) ) / ( 10000 * ${sdensity} ) ) ) >= ${old_soil}, ( -1 * ${old_soil} ), ( ( (${qsxdx} + ${qsydy}) )/ ( 10000 * ${sdensity} ) ) ), ( ( (${qsxdx} + ${qsydy}) ) / ( 10000 * ${sdensity} ) ) )', erdep = erdep, old_soil = old_soil, qsxdx = qsxdx, sdensity = sdensity, qsydy = qsydy)
+    grass.message('\n*************************\n Year %s -- ' % o + 'step 4: calculating divergence/difference of sediment transport for each process and the actual amount of erosion or deposition in vertical meters/cell/year\n*************************\n\n')
+    #Here is where we figure out the change in transport capacity, and thus the actual amount of erosion an deposition that would occur. There are two ways of doing this. On planar and convex surfaces (i.e., ridgetops, flats, hillslopes), it is better to take the 2D divergence of sediment flux (we use r.slope.aspect to calculate this), but on highly convex surfaces (i.e., in channels) it is better to take the 1D difference between one cell, and the cell that is immediately downstream from it. This all assumes that the system is always operating at Transport Capacity, or if it is not, then is still behaves as if it were (ie., that the actual differences in transported sediment between the cells would be proportional to the system operating at capacity). Thus, under this assumption, the divergence of capacity is equals to actual amount of sediment eroded/deposited.  
     
-    grass.message('\n*************************\n Year %s ' % o + 'step 7 of 7: calculating terrain evolution, new bedrock elevations, and new soil depths\n *************************\n\n')
-    #put the net dz back where it is supposed to go (correct for the fact that, in grass, derivitaves of slope are calculated and stored one cell downslope from the cell where they actually belong) must be calulatedand then subtract it from dem to make new dem
-    grass.mapcalc('${netchange} =eval(x=if((${aspect}  < 22.5 ||  ${aspect}  >= 337.5) && ${aspect}  != 0, (${erdep} [1,0]), if (${aspect}  >= 22.5 && ${aspect}  < 67.5, (${erdep} [1,-1]), if (${aspect}  >= 67.5 && ${aspect}  < 112.5, (${erdep} [0,-1]), if (${aspect}  >= 112.5 && ${aspect}  < 157.5, (${erdep} [-1,-1]), if (${aspect}  >= 157.5 && ${aspect}  < 202.5, (${erdep} [-1,0]), if (${aspect}  >= 202.5 && ${aspect}  < 247.5, (${erdep} [-1,1]), if (${aspect}  >= 247.5 && ${aspect}  < 292.5, (${erdep} [0,1]), if (${aspect}  >= 292.5 && ${aspect}  < 337.5, (${erdep} [1,1]), (${erdep} ))))))))), (if(isnull(x), ${erdep} , x)))', netchange = netchange, aspect = aspect, erdep = erdep)
-
-    temp_dem = '%stemp_dem' % p
-    grass.mapcalc('' + temp_dem + '=' + old_dem + ' + ' + netchange + '', temp_dem = temp_dem, old_dem = old_dem, netchange = netchange)
-    #do patch-job to catch the shrinking edge problem
+    #This is the way we implemnt this: First calculate, we calculate the divergence/differnce for EACH of the different flow processes on the ENTIRE map (i.e., make one map per process, difference for streams, divergence for USPED and diffusion). Then, we cut out the pieces of each of these maps that correspond to the correct landforms from each specific process (based on the user-input cutoffs in flow accumulation), and patch them together into a single map (NOTE: see output unit conversions section below to see how we get all the units to line up during this process). This counters the "boundary effect" that happens when running the differential equations for divergence across the boundary of two different flow processes.  Then we may still have to run a median smoother on the patched map to get rid of any latent spikes.
+    
+    qsxdx4 = '%sDelta_Qsx_diffusion%04d' % (p, o)
+    qsydy4 = '%sDelta_Qsy_diffusion%04d' % (p, o)
+    grass.run_command('r.slope.aspect',  quiet = True,  elevation = qsx4, dx = qsxdx4)
+    grass.run_command('r.slope.aspect',  quiet = True,  elevation = qsy4, dy = qsydy4)
+    #this commented out section is the 1D difference for this process
+    #grass.mapcalc('${nc4}=if(${flowdir} == 7, (${qs4}[-1,-1]-${qs4}), if (${flowdir} == 6, (${qs4}[-1,0]-${qs4}), if (${flowdir} == 5, (${qs4}[-1,1]-${qs4}), if (${flowdir} == 4, (${qs4}[0,1]-${qs4}), if (${flowdir} == 3, (${qs4}[1,1]-${qs4}), if (${flowdir} == 2, (${qs4}[1,0]-${qs4}), if (${flowdir} == 1, (${qs4}[1,-1]-${qs4}), if (${flowdir} == 8, (${qs4}[0,-1]-${qs4}), ${qs4}))))))))', nc4 = nc4, qsxdx4 = qsydy4)
+    if cutoff1 == cutoff2:
+        qsxdx3 = qsxdx4
+        qsydy3 = qsydy4
+    else:
+        qsxdx3 = '%sDelta_Qsx_sheetwash%04d' % (p, o)
+        qsydy3 = '%sDelta_Qsy_sheetwash%04d' % (p, o)
+        grass.run_command('r.slope.aspect',  quiet = True,  elevation = qsx3, dx = qsxdx3)
+        grass.run_command('r.slope.aspect',  quiet = True,  elevation = qsy3, dy = qsydy3)
+        #this commented out section is the 1D difference for this process        
+        #grass.mapcalc('${nc3}=if(${flowdir} == 7, (${qs3}[-1,-1]-${qs3}), if (${flowdir} == 6, (${qs3}[-1,0]-${qs3}), if (${flowdir} == 5, (${qs3}[-1,1]-${qs3}), if (${flowdir} == 4, (${qs3}[0,1]-${qs3}), if (${flowdir} == 3, (${qs3}[1,1]-${qs3}), if (${flowdir} == 2, (${qs3}[1,0]-${qs3}), if (${flowdir} == 1, (${qs3}[1,-1]-${qs3}), if (${flowdir} == 8, (${qs3}[0,-1]-${qs3}), ${qs3}))))))))', nc3 = nc3, flowdir = flowdir, qs3 = qs3)
+    if cutoff2 == cutoff3:
+        qsxdx2 = qsxdx3
+        qsydy2 = qsydy3
+    else:
+        qsxdx2 = '%sDelta_Qsx_rills%04d' % (p, o)
+        qsydy2 = '%sDelta_Qsy_rills%04d' % (p, o)
+        grass.run_command('r.slope.aspect',  quiet = True,  elevation = qsx2, dx = qsxdx2)
+        grass.run_command('r.slope.aspect',  quiet = True,  elevation = qsy2, dy = qsydy2)
+        #this commented out section is the 1D difference for this process
+        #grass.mapcalc('${nc2}=if(${flowdir} == 7, (${qs2}[-1,-1]-${qs2}), if (${flowdir} == 6, (${qs2}[-1,0]-${qs2}), if (${flowdir} == 5, (${qs2}[-1,1]-${qs2}), if (${flowdir} == 4, (${qs2}[0,1]-${qs2}), if (${flowdir} == 3, (${qs2}[1,1]-${qs2}), if (${flowdir} == 2, (${qs2}[1,0]-${qs2}), if (${flowdir} == 1, (${qs2}[1,-1]-${qs2}), if (${flowdir} == 8, (${qs2}[0,-1]-${qs2}), ${qs2}))))))))', nc2 = nc2, flowdir = flowdir, qs2 = qs2)
+    
+    if ( os.getenv("GIS_FLAG_1") == "1" ):
+        #this is the 1D difference for this process
+        qsd1 = '%sDelta_Qs_1D_streams%04d' % (p, o)
+        grass.mapcalc('${qsd1}=if(${flowdir} == 7, (${qs1}[-1,-1]-${qs1}), if (${flowdir} == 6, (${qs1}[-1,0]-${qs1}), if (${flowdir} == 5, (${qs1}[-1,1]-${qs1}), if (${flowdir} == 4, (${qs1}[0,1]-${qs1}), if (${flowdir} == 3, (${qs1}[1,1]-${qs1}), if (${flowdir} == 2, (${qs1}[1,0]-${qs1}), if (${flowdir} == 1, (${qs1}[1,-1]-${qs1}), if (${flowdir} == 8, (${qs1}[0,-1]-${qs1}), ${qs1}))))))))', qsd1 = qsd1, flowdir = flowdir, qs1 = qs1)
+    else:
+        #this is the 2D version
+        qsxdx1 = '%sDelta_Qsx_streams%04d' % (p, o)
+        qsydy1 = '%sDelta_Qsy_streams%04d' % (p, o)
+        grass.run_command('r.slope.aspect',  quiet = True,  elevation = qsx1, dx = qsxdx1)
+        grass.run_command('r.slope.aspect',  quiet = True,  elevation = qsy1, dy = qsydy1)
+    
+    #This is the smoothing routine. First we calculate the rate of Erosion and Deposition by converting the Delta QS of the different processes to vertical meters by dividing by the soil denisity (with apropriate constants to get into the correct units, see UNIT CONVERSION note below), and for streams, also expand from the storm to the year level. All units of this initial (temporary) ED_rate map will be in m/cell/year. 
+    #OUTPUT UNIT CONVERSIONS: In the case of the diffusion equation, the output units are in verticle meters of sediment per cell per year, so these will be left alone. In the case of stream flow, the output units are kg/m2/storm, so need to multiply by 1000 to get T/m2/storm, and then divide by the soil density (T/m3) to get verticle meters of sediment/cell/storm, and will be multiplied by the number of storms/year in order to get vertical meters of sediment/cell/year. In the case of USPED, the output is in T/Ha/year, so first multiply by 0.1 to get T/m2/year and then divide by soil density (T/m3) to get verticle meters of sediment/cell/year. 
+    tempnetchange1 = '%sTEMPORARY_UNSMOOTHED_ED_rate%04d' % (p, o)
+    if ( os.getenv("GIS_FLAG_1") == "1" ):
+        grass.mapcalc('${tempnetchange1}=if(${flowacc} >= ${cutoff3}, ((${qsd1})/(${sdensity}*1000))*(${storms}*0.25*${stormtimet}), if(${flowacc} >= ${cutoff2} && ${flowacc} < ${cutoff3}, ((${qsxdx2}+${qsydy2})*0.1)/${sdensity}, if(${flowacc} >= ${cutoff1} && ${flowacc} < ${cutoff2}, ((${qsxdx3}+${qsydy3})*0.1)/${sdensity}, ${qsxdx4}+${qsydy4})))',  tempnetchange1 = tempnetchange1, qsd1 = qsd1, qsxdx2 = qsxdx2, qsydy2 = qsydy2, qsxdx3 = qsxdx3, qsydy3 = qsydy3, qsxdx4 = qsxdx4, qsydy4 = qsydy4, flowacc = flowacc, cutoff1 = cutoff1,  cutoff2 = cutoff2,  cutoff3 = cutoff3, sdensity = sdensity, storms = storms, stormtimet = stormtimet)
+    else:
+        grass.mapcalc('${tempnetchange1}=if(${flowacc} >= ${cutoff3}, ((${qsxdx1} + ${qsydy1})/(${sdensity}*1000))*(${storms}*0.25*${stormtimet}), if(${flowacc} >= ${cutoff2} && ${flowacc} < ${cutoff3}, ((${qsxdx2}+${qsydy2})*0.1)/${sdensity}, if(${flowacc} >= ${cutoff1} && ${flowacc} < ${cutoff2}, ((${qsxdx3}+${qsydy3})*0.1)/${sdensity}, ${qsxdx4}+${qsydy4})))',  tempnetchange1 = tempnetchange1, qsxdx1 = qsxdx1, qsydy1 = qsydy1, qsxdx2 = qsxdx2, qsydy2 = qsydy2, qsxdx3 = qsxdx3, qsydy3 = qsydy3, qsxdx4 = qsxdx4, qsydy4 = qsydy4, flowacc = flowacc, cutoff1 = cutoff1,  cutoff2 = cutoff2,  cutoff3 = cutoff3, sdensity = sdensity, storms = storms, stormtimet = stormtimet)
+    #Make some temp maps of just erosion rate and just deposition rate so we can grab some stats from them (If Smoothing is 'no', then we'll also use these stats later in the stats file)
+    grass.message('Running presmoothing filters...')
+    tmperosion = p + 'tmperosion%04d' % o
+    tmpdep = p + 'tmpdep%04d' % o
+    grass.mapcalc('${tmperosion}=if(${tempnetchange1} < -0, ${tempnetchange1}, null())', tmperosion = tmperosion, tempnetchange1 = tempnetchange1)
+    grass.mapcalc('${tmpdep}=if(${tempnetchange1} > 0, ${tempnetchange1}, null())', tmpdep = tmpdep, tempnetchange1 = tempnetchange1)
+    #Grab the stats from these temp files and save them to dictionaries
+    erosstats = grass.parse_command('r.univar', flags = 'g', map = tmperosion)
+    depostats = grass.parse_command('r.univar', flags = 'g', map = tmpdep)
+    maximum = depostats['max']
+    minimum = erosstats['min']
+    erosbreak =  float(erosstats['mean']) - float(erosstats['stddev'])
+    deposbreak = float(depostats['mean']) + float(depostats['stddev'])
+    scalemin = float(erosstats['mean']) - (2.0*float(erosstats['stddev']))
+    scalemax = float(depostats['mean']) + (2.0*float(depostats['stddev']))
+    #Use the stats we gathered to do some smoothing with a hi-cut and lo-cut filter (with soft-knee limiting) of the unsmoothed ED_rate map. Values from 1 standard dev below the mean erosion to 1 standard deviation above the mean deposition will remain the same. Values above that will be scaled linearly so that the cell(s) that had the original maximum value of erosion (some value much much greater than the mean) will be rescaled to two standad deviations below the mean erosion, and vice versa for maximum deposition. This brings any values that were really unreasonnable as originally calculated (spikes) into the range of what the maximum values should be on a normally distrubuted dataset, but does so with out a "brick wall" style of limiting, which would make all values above some cutoff equal to a theoretical maximum. Instead, this "soft-knee" style of limiting allows some scaling at the high ends, which allows for the smoothed value of very high cells to still be relativel
 y higher than other cells that are also above the cutoff, but not as high as those very high cells.
+    tempnetchange2 = '%sTEMPORARY_PRESMOOTHED_ED_rate%04d' % (p, o)
+    grass.mapcalc('${tempnetchange2}=graph(${tempnetchange1}, ${minimum},${scalemin}, ${erosbreak},${erosbreak}, ${deposbreak},${deposbreak}, ${maximum},${scalemax})', tempnetchange2 = tempnetchange2, tempnetchange1 =tempnetchange1, minimum = minimum, scalemin = scalemin, erosbreak = erosbreak, deposbreak = deposbreak, maximum = maximum, scalemax = scalemax)
+    
+    #Check if additional modal smoothing is requested (and if so the desired median smoothing neighborhood), and do what it asks (using mapcalc because it's a little faster than r.neighbors)
+    if len(os.getenv("GIS_OPT_smoothing")) is 2:
+        grass.message('Enacting "no" modal smoothing...')
+        grass.run_command('g.rename',  quiet = True, rast = tempnetchange2 + ',' + netchange)
+    elif len(os.getenv("GIS_OPT_smoothing")) is 3:
+        grass.message('Enacting "low" (3x3) modal smoothing.')
+#        grass.mapcalc('${tmperosion}=if(${tempnetchange2} < -0, ${tempnetchange2}, null())', tmperosion = tmperosion, tempnetchange2 = tempnetchange2)
+#        grass.mapcalc('${tmpdep}=if(${tempnetchange2} > 0, ${tempnetchange2}, null())', tmpdep = tmpdep, tempnetchange2 = tempnetchange2)
+#        erosstats2 = grass.parse_command('r.univar', flags = 'g', map = tmperosion)
+#        depostats2 = grass.parse_command('r.univar', flags = 'g', map = tmpdep)
+#        erosbreak1 =  float(erosstats2['mean']) - float(erosstats2['stddev'])
+#        deposbreak1 = float(depostats2['mean']) + float(depostats2['stddev'])
+        grass.run_command('r.neighbors', quiet = True, input = tempnetchange2, output = netchange, method = 'mode', size = '3')
+        #grass.mapcalc("${netchange}=if(${tempnetchange2} > ${deposbreak1} || ${tempnetchange2} < ${erosbreak1}, mode(${tempnetchange2}[-1,-1], ${tempnetchange2}[-1,0], ${tempnetchange2}[-1,1], ${tempnetchange2}[0,-1], ${tempnetchange2}[0,1], ${tempnetchange2}[1,-1], ${tempnetchange2}[1,0], ${tempnetchange2}[1,1]))", netchange = netchange, tempnetchange2 = tempnetchange2, deposbreak1 = deposbreak1, erosbreak1 = erosbreak1)
+        #Grab the stats from these new smoothed netchange maps and save them to dictionaries (Note that these are overwriting the stats variables previously gathered for the hi/lo cut filter)
+        grass.mapcalc('${tmperosion}=if(${netchange} < -0, ${netchange}, null())', tmperosion = tmperosion, netchange = netchange)
+        grass.mapcalc('${tmpdep}=if(${netchange} > 0, ${netchange}, null())', tmpdep = tmpdep, netchange = netchange)
+        erosstats1 = grass.parse_command('r.univar', flags = 'g', map = tmperosion)
+        depostats1 = grass.parse_command('r.univar', flags = 'g', map = tmpdep)
+    elif len(os.getenv("GIS_OPT_smoothing")) is 4:
+        grass.message('Enacting "high" (5x5) modal smoothing.')
+        grass.run_command('r.neighbors', quiet = True, input = tempnetchange2, output = netchange, method = 'mode', size = '5')
+        #grass.mapcalc("${netchange}= mode(${tempnetchange2}[-1,-1], ${tempnetchange2}[-1,0], ${tempnetchange2}[-1,1], ${tempnetchange2}[0,-1], ${tempnetchange2}[0,1], ${tempnetchange2}[1,-1], ${tempnetchange2}[1,0], ${tempnetchange2}[1,1], ${tempnetchange2}[-2,-2], ${tempnetchange2}[-2,-1], ${tempnetchange2}[-2,0], ${tempnetchange2}[-2,1], ${tempnetchange2}[-2,2], ${tempnetchange2}[-1,2], ${tempnetchange2}[0,2], ${tempnetchange2}[1,2], ${tempnetchange2}[2,2], ${tempnetchange2}[2,1], ${tempnetchange2}[2,0], ${tempnetchange2}[2,-1], ${tempnetchange2}[2,-2], ${tempnetchange2}[1,-2], ${tempnetchange2}[0,-2], ${tempnetchange2}[-1,-2])", netchange = netchange, tempnetchange2 = tempnetchange2)
+        #Grab the stats from these new smoothed netchange maps and save them to dictionaries (Note that these are overwriting the stats variables previously gathered for the hi/lo cut filter)
+        grass.mapcalc('${tmperosion}=if(${netchange} < -0, ${netchange}, null())', tmperosion = tmperosion, netchange = netchange)
+        grass.mapcalc('${tmpdep}=if(${netchange} > 0, ${netchange}, null())', tmpdep = tmpdep, netchange = netchange)
+        erosstats1 = grass.parse_command('r.univar', flags = 'g', map = tmperosion)
+        depostats1 = grass.parse_command('r.univar', flags = 'g', map = tmpdep)
+    else:
+        grass.message('There was a problem reading the median-smoothing variable, so maps will not be median-smoothed.')
+        grass.run_command('g.rename',  quiet = True, rast = tempnetchange2 + ',' + netchange)    
+    #Set the netchange map colors to the rules we've provided above
+    grass.run_command('r.colors', quiet = True, map = netchange, rules = nccolors.name)
+    #Clean up temp maps
+    grass.run_command('g.remove',  quiet = True, rast = tmperosion + ',' + tmpdep + ',' + tempnetchange1 + ',' + tempnetchange2)
+    
+    grass.message('\n*************************\n Year %s -- ' % o + 'step 5: calculating terrain evolution and new soil depths\n *************************\n\n')
+    #Set up a temp dem, and then do initial addition of ED change to old DEM. This mapcalc statement first 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 (could happen, I suppose), and if it is, corrects it). 
+    temp_dem = '%stemp_dem%04d' % (p, o)
+    grass.mapcalc('${temp_dem}=eval(x=if(${old_soil} > 0.0 && (-1*${netchange}) <= ${old_soil}, ${netchange}, if((-1*${netchange}) > ${old_soil},   (-1*${old_soil}), 0)), ${old_dem} + x)', temp_dem = temp_dem, old_soil = old_soil, old_dem = old_dem,  netchange = netchange)
+    #Do patch-job to catch the shrinking edge problem (the edge cells have no upstream cell, so get turned null in the calculations in step 4)
     grass.run_command('r.patch', quiet = True,  input = temp_dem + ',' + old_dem,  output= new_dem)
-    #set colors for elevation map to match other dems
+    #Set colors for elevation map to match other dems and then get rid of the temp dem
     grass.run_command('r.colors',  quiet = True, map = new_dem, rast = os.getenv("GIS_OPT_elev"))
     grass.run_command('g.remove', quiet = True, rast = temp_dem)
-    #if asked, calculate amount of bedrock weathered and update soildepths map with this info, else just make a new soildepths map based on amount of erosion deposition
-    if ( os.getenv("GIS_FLAG_b") == "1" ):
-        grass.message('\nstep 8.5: Experimental bedrock weathering\n')
-        grass.mapcalc('${meancurv}=((${pc} + ${tc}) / 2)', meancurv = meancurv, pc = pc,  tc = tc)
-        # create dictionary to record max and min curvature
-        statdict2 = grass.parse_command('r.info',  flags = 'r', map =meancurv)
-        grass.message('\nThe raw max (' + statdict2['max'] + ') and min (' + statdict2['min'] + ') curvature will be rescaled from 2 to 0\n')
-        # create map of bedrock weathering rates
-        grass.mapcalc('${rate}=${kappa}*(2-(${meancurv}*(2/(${max})-(${min}))))', rate = rate, kappa = kappa, meancurv = meancurv,  max = statdict2['max]'], min = statdict2['min'])
-        #rate is actually the net change in bedrock elevation due to soil production, so lets use it to find the new bedrock elev, and the new soil depth!
-        grass.mapcalc('${new_bdrk}=${old_bdrk} - ${rate}', new_bdrk = new_bdrk, old_bdrk = old_bdrk, rate = rate)
-        grass.message('Calculating new soil depths using new bedrock map')
-        grass.mapcalc('${new_soil}=if ((${new_dem} - ${new_bdrk}) < 0, 0, (${new_dem} - ${new_bdrk}))',  new_soil = new_soil, new_dem = new_dem, new_bdrk = new_bdrk)
-        #these are the old soil equations that I failed to be able to implement... I leave them in for documentation purposes
-        #r.mapcalc "$new_bdrk=$initbdrk - ($Ba * ($Bb*($erdep - $initbdrk)))"
-        #r.mapcalc "$new_soil=if (($erdep - $initbdrk) < 0, 0, ($erdep - $initbdrk))"
-        grass.run_command('g.remove', quiet = True, rast = rate + ',' + meancurv)
-    else:
-        grass.message('\nCalculating new soil depths keeping bedrock elevations static\n')
-        grass.mapcalc('${new_soil}=if ((${new_dem} - ${initbdrk}) < 0, 0, (${new_dem} - ${initbdrk}))', new_soil = new_soil, new_dem = new_dem, initbdrk = initbdrk)
-    #setting colors for new soil depth map
-    grass.message('reading color rules from %s' % sdcolors.name)
-    grass.message (sdcolors.readlines())
+    #calculate the new soildepth map and set colors
+    grass.mapcalc('${new_soil}=if ((${new_dem} - ${initbdrk}) < 0, 0, (${new_dem} - ${initbdrk}))', new_soil = new_soil, new_dem = new_dem, initbdrk = initbdrk)
     grass.run_command('r.colors', quiet = True, map = new_soil, rules = sdcolors.name)
-    #make some temp maps of just erosion and just deposition so we can grab some stats from them
-    tmperosion = p + 'tmperosion%s' % o
-    tmpdep = p + 'tmpdep%s' % o
-    grass.mapcalc('' + tmperosion + '=if(' + netchange + ' < 0, ' + netchange + ', null())', tmperosion = tmperosion, netchange = netchange)
-    grass.mapcalc('' + tmpdep + '=if(' + netchange + ' > 0, ' + netchange + ', null())', tmpdep = tmpdep, netchange = netchange)
-
-    #grab the stats fromt hese temp files and save them to dictionaries
+    grass.message('\n*************************\n Year %s -- ' % o + 'step 6: writing stats to output file\n *************************\n\n')
+    #Finish gathering stats (just need the soil depth stats now)
+    soilstats = {}
     soilstats = grass.parse_command('r.univar', flags = 'ge', map = new_soil, percentile = '99')
-    erosstats = grass.parse_command('r.univar', flags = 'ge', map = tmperosion, percentile = '1')
-    depostats = grass.parse_command('r.univar', flags = 'ge', map = tmpdep,  percentile = '99')
-    #clean up temp maps
-    grass.run_command('g.remove',  quiet = True, rast = tmperosion + ',' + tmpdep)
-    #write stats to a new line in the stats file
+    #Write stats to a new line in the stats file
+    #HEADER of the file should be: "Year,Mean Erosion,Standard Deviation Erosion,Smoothed Maximum Erosion,Original Maximum Erosion,,Mean Deposition,Standard Deviation Deposition,Smoothed Maximum Deposition,Original Maximum Deposition,,Mean Soil Depth, Standard Deviation Soil Depth,Minimum Soil Depth, Maximum Soil Depth"
     grass.message('Outputing stats to textfile: ' + q)
-    f.write('\n%s' % o + ',,' + erosstats['mean'] + ',' + erosstats['min'] + ',' + erosstats['max'] + ',' + erosstats['percentile_1'] + ',,' + depostats['mean'] + ',' + depostats['min'] + ',' + depostats['max'] + ',' + depostats['percentile_99'] + ',,' + soilstats['mean'] + ',' + soilstats['min'] + ',' + soilstats['max'] + ',' + soilstats['percentile_99'])
-    #delete netchange map, if asked to, otherwise set colors for output netchange map
-    if ( os.getenv("GIS_FLAG_n") == "1" ):
-        grass.message('Not keeping a netchange map')
-        grass.run_command('g.remove', quiet = True, rast = netchange)
+    if len(os.getenv("GIS_OPT_smoothing")) is 3 or len(os.getenv("GIS_OPT_smoothing")) is 4:
+        f.write('\n%s' % o + ',,' + erosstats1['mean'] + ',' + erosstats1['stddev'] + ',' +erosstats1['min'] + ',' + minimum + ',,' + depostats1['mean'] + ',' + depostats1['stddev'] + ',' + depostats1['max'] + ',' + maximum + ',,' + soilstats['mean'] + ',' + soilstats['stddev'] + ',' + soilstats['min'] + ',' + soilstats['max'])
     else:
-        grass.run_command('r.colors', quiet = True, map = netchange, rules = nccolors.name)
+        f.write('\n%s' % o + ',,' + erosstats['mean'] + ',' + erosstats['stddev'] + ',' + str(scalemin) + ',' + minimum + ',,' + depostats['mean'] + ',' + depostats['stddev'] + ',' + str(scalemax) + ',' + maximum + ',,' + soilstats['mean'] + ',' + soilstats['stddev'] + ',' + soilstats['min'] + ',' + soilstats['max'])
+
+    #Clean up temporary maps
+    if os.getenv("GIS_FLAG_k") == "1":
+        grass.message('\nTemporary maps will NOT be deleted!!!!\n')
+    else:
+        grass.message('\nCleaning up temporary maps...\n\n')
+        if os.getenv("GIS_FLAG_s") == "1":
+            grass.message('Keeping Slope map.')
+        else:
+            grass.run_command('g.remove', quiet =True, rast = slope)
+        if os.getenv("GIS_FLAG_d") == "1":
+            grass.message('Not keeping Soil Depth map.')
+            grass.run_command('g.remove', quiet =True, rast = old_soil)
+            #check if this is the last year and remove the "new-soil" map too
+            if ( o == int(os.getenv("GIS_OPT_number"))):
+                grass.run_command('g.remove', quiet =True, rast = new_soil)
+        else:
+            #check if this is the first year, and if so, remove the temporary "soildepths_init" map
+            if ( o == 1 ):
+                grass.run_command('g.remove', quiet = True, rast = old_soil)
+        if ( os.getenv("GIS_FLAG_e") == "1" ):
+            grass.message('Keeping Excess Transport Capacity (divergence) maps for all processes.')
+        else:
+            if ( os.getenv("GIS_FLAG_1") == "1" ):
+                grass.run_command('g.remove', quiet = True, rast = qsxdx4 + ',' + qsydy4 + ',' + qsxdx2 + ',' + qsydy2 + ',' + qsxdx3 + ',' + qsydy3 + ',' +  qsd1)
+            else:
+                grass.run_command('g.remove', quiet = True, rast = qsxdx4 + ',' + qsydy4 + ',' + qsxdx2 + ',' + qsydy2 + ',' + qsxdx3 + ',' + qsydy3 + ',' +  qsxdx1 + ',' +  qsydy1)
+        if ( os.getenv("GIS_FLAG_t") == "1" ):
+            grass.message('Keeping Transport Capacity maps for all processes.')
+        else:
+            if ( os.getenv("GIS_FLAG_1") == "1" ):
+                grass.run_command('g.remove', quiet = True, rast = qsx4 + ',' + qsy4 + ',' + qsx2 + ',' + qsy2 + ',' + qsx3 + ',' + qsy3 + ',' +  qs1)
+            else:
+                grass.run_command('g.remove', quiet = True, rast = qsx4 + ',' + qsy4 + ',' + qsx2 + ',' + qsy2 + ',' + qsx3 + ',' + qsy3 + ',' +  qsx1 + ',' +  qsy1)
+        if ( os.getenv("GIS_FLAG_r") == "1" ):
+            grass.message('Not keeping an Erosion and Deposition rate map.')
+            grass.run_command('g.remove', quiet = True, rast = netchange)
+        grass.run_command('g.remove', quiet =True, rast = flowdir + ',' + flowacc + ',' + aspect)
     sdcolors.close()
     nccolors.close()
     grass.message('\n*************************\nDone with Year %s ' % o + '\n*************************\n')
-    if ( m == '0' and os.getenv("GIS_FlAG_e") == "1" ):
-        grass.message('\nkeeping initial soil depths map ' + old_soil + '\n')
-    elif m == '0':
-        grass.run_command('g.remove', quiet = True, rast = old_soil)
-    grass.message('\nIf made, raster map ' + netchange + ' shows filtered net erosion/deposition\nIf made, raster map ' + new_soil + ' shows soildpeths\nRaster map ' + new_dem + ' shows new landscape (new elevations) after net erosion/depostion\n*************************\n')
 
-
-# here is where the code in "main" actually gets executed. This way of programming is neccessary for the way g.parser needs to run.
+#Here is where the code in "main" actually gets executed. This way of programming is neccessary for the way g.parser needs to run.
 if __name__ == "__main__":
     if ( len(sys.argv) <= 1 or sys.argv[1] != "@ARGS_PARSED@" ):
         os.execvp("g.parser", [sys.argv[0]] + sys.argv)
     else:
-        # set up some basic variables
+        # Set up some basic variables
         years = os.getenv("GIS_OPT_number")
         prefx = os.getenv("GIS_OPT_prefx")
-        #make the stats out file with correct column headers
+        #Make the statsout file with correct column headers
         if os.getenv("GIS_OPT_statsout") == "":
             env = grass.gisenv()
             mapset = env['MAPSET']
             statsout = '%s_%s_lsevol_stats.txt' % (mapset, prefx)
         else:
             statsout = os.getenv("GIS_OPT_statsout")
-        f = file(statsout, 'wt')
-        f.write('Year,,Mean Erosion,Max Erosion,Min Erosion,99th Percentile Erosion,,Mean Deposition,Min Deposition,Max Deposition,99th Percentile Deposition,,Mean Soil Depth,Min Soil Depth,Max Soil Depth,99th Percentile Soil Depth')
-        #region MUST be set to the map being used. otherwise r.flow will not work.
-        #let's grab the current region
+        if os.path.isfile(statsout):
+            f = file(statsout, 'a')
+        else:
+            f = file(statsout, 'wt')
+            f.write('These statistics are in units of vertical meters (depth) per cell\nYear,,Mean Erosion,Standard Deviation Erosion,Smoothed Maximum Erosion,Original Maximum Erosion,,Mean Deposition,Standard Deviation Deposition,Smoothed Maximum Deposition,Original Maximum Deposition,,Mean Soil Depth, Standard Deviation Soil Depth,Minimum Soil Depth, Maximum Soil Depth')
+        if ( os.getenv("GIS_FLAG_p") == "1" ):
+            grass.message('Making sample points map for determining cutoffs.')
+        else:
+            grass.message('\n##################################################\n##################################################\n\n STARTING SIMULATION\n\nBeginning iteration sequence. This may take some time.\nProcess is not finished until you see the message: \'Done with everything\'\n _____________________________________________________________\n_____________________________________________________________\n')
+            grass.message("Total number of iterations to be run is %s years" % years)
+        #Check to see if there is a MASK already, and if so, temporarily rename it
+        if "MASK" in grass.list_grouped('rast')[grass.gisenv()['MAPSET']]:
+            ismask = 1
+            tempmask = 'mask_%i' % random.randint(0,100)
+            grass.run_command('g.rename', quiet = "True", rast = 'MASK,' + tempmask)
+        else:
+            ismask = 0
+        #Set MASK to input DEM
+        grass.run_command('r.mask', input = os.getenv("GIS_OPT_elev"))
+        #Get the region settings
         region1 = grass.region()
-        #now get map extents and resolution of input maps
-        mapextent = grass.raster_info(os.getenv("GIS_OPT_elev"))
-        #Now test to see if region extents and res matches map extents and res
-        if (region1['e'] != mapextent['east'] or region1['w'] != mapextent['west'] or region1['n'] != mapextent['north'] or region1['s'] != mapextent['south'] or region1['nsres'] != mapextent['nsres'] or region1['ewres'] != mapextent['ewres'] ):
-            grass.fatal('Region MUST be set to match input map: ' + os.getenv("GIS_OPT_elev") + '/n Please run g.region and then try again!')
-        grass.message('\n##################################################\n##################################################\n\n STARTING SIMULATION\n\nBeginning iteration sequence. This may take some time.\nProcess is not finished until you see the message: \'Done with everything\'\n _____________________________________________________________\n_____________________________________________________________\n')
-        grass.message("Total number of iterations to be run is %s years" % years)
         # This is the loop!
         for x in range(int(years)):
             grass.message("Iteration = %s" % (x + 1))
             main(x, (x + 1), prefx, statsout,  region1['nsres']);
-        #since we are now done with the loop, close the stats file.
+        #Since we are now done with the loop, close the stats file and remove the MASK, and optionally replace the old mask
         f.close()
-        grass.message('\nIterations complete, keeping region set to output maps\n')
-        if os.getenv("GIS_FLAG_k") == "1":
-            grass.message('\nTemporary maps have NOT been deleted\n')
-        else:
-            grass.message('\nCleaning up temporary maps...\n\n')
-            if os.getenv("GIS_FLAG_p") == "0":
-                grass.run_command('g.mremove', quiet =True, flags = 'f', rast = prefx + 'slope*')
-            grass.run_command('g.mremove', quiet =True, flags = 'f', rast = prefx + 'flowacc*')
-            grass.run_command('g.mremove', quiet =True, flags = 'f', rast = prefx + 'aspect*')
-            grass.run_command('g.mremove', quiet =True, flags = 'f', rast = prefx + 'sflowtop*')
-            grass.run_command('g.mremove', quiet =True, flags = 'f', rast = prefx + 'qsx*')
-            grass.run_command('g.mremove', quiet =True, flags = 'f', rast = prefx + 'qsy*')
-            grass.run_command('g.mremove', quiet =True, flags = 'f', rast = prefx + 'erosdep*')
-            grass.run_command('g.mremove', quiet =True, flags = 'f', rast = prefx + 'pc*')
-            grass.run_command('g.mremove', quiet =True, flags = 'f', rast = prefx + 'tc*')
-            grass.message('\nDone\n\n')
-        grass.message("\n\nDone with everything")
+        grass.run_command('r.mask', flags = 'r')
+        if ismask == 1:
+            grass.run_command('g.rename', quiet = "True", rast = tempmask +',MASK')
+        grass.message('\nIterations complete!\n\nDone with everything')
 
 
 
+



More information about the grass-commit mailing list