[GRASS-SVN] r41992 - in grass-addons/LandDyn:
devs_landcover_scripts/r.cfactor.py
devs_landcover_scripts/r.landcover.update.py
devs_landcover_scripts/r.soil.fertility.py
devs_landcover_scripts/r.villages.py r.accumulate.erdep.py
r.catchment.py r.landscape.evol.py r.soildepth.py
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu Apr 22 16:00:17 EDT 2010
Author: isaacullah
Date: 2010-04-22 16:00:17 -0400 (Thu, 22 Apr 2010)
New Revision: 41992
Added:
grass-addons/LandDyn/devs_landcover_scripts/r.cfactor.py/r.cfactor.py
grass-addons/LandDyn/devs_landcover_scripts/r.landcover.update.py/r.landcover.update.py
grass-addons/LandDyn/devs_landcover_scripts/r.soil.fertility.py/r.soil.fertility.py
grass-addons/LandDyn/devs_landcover_scripts/r.villages.py/r.villages.py
grass-addons/LandDyn/r.accumulate.erdep.py/r.accumulate.erdep.py
grass-addons/LandDyn/r.catchment.py/r.catchment.py
grass-addons/LandDyn/r.landscape.evol.py/Makefile
grass-addons/LandDyn/r.landscape.evol.py/r.landscape.evol description.odf_html_585d862d.gif
grass-addons/LandDyn/r.landscape.evol.py/r.landscape.evol description.odf_html_m100fb7e.gif
grass-addons/LandDyn/r.landscape.evol.py/r.landscape.evol description.odf_html_m11de82c.gif
grass-addons/LandDyn/r.landscape.evol.py/r.landscape.evol description.odf_html_m2c6cce6a.gif
grass-addons/LandDyn/r.landscape.evol.py/r.landscape.evol description.odf_html_m2f9c13ec.gif
grass-addons/LandDyn/r.landscape.evol.py/r.landscape.evol description.odf_html_m8e0f3ca.gif
grass-addons/LandDyn/r.landscape.evol.py/r.landscape.evol.py
grass-addons/LandDyn/r.landscape.evol.py/r.landscape.evol.py.html
grass-addons/LandDyn/r.soildepth.py/r.soildepth.py
Log:
Actually adding the files!
Added: grass-addons/LandDyn/devs_landcover_scripts/r.cfactor.py/r.cfactor.py
===================================================================
--- grass-addons/LandDyn/devs_landcover_scripts/r.cfactor.py/r.cfactor.py (rev 0)
+++ grass-addons/LandDyn/devs_landcover_scripts/r.cfactor.py/r.cfactor.py 2010-04-22 20:00:17 UTC (rev 41992)
@@ -0,0 +1,94 @@
+#!/usr/bin/python
+#
+############################################################################
+#
+# MODULE: r.cfactor
+# AUTHOR(S): Isaac Ullah, Michael Barton, Arizona State University
+# PURPOSE: Converts a map of landcover values to a c-factor map based
+# on a set of reclass rules
+# ACKNOWLEDGEMENTS: National Science Foundation Grant #BCS0410269
+# COPYRIGHT: (C) 2009 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: Converts a map of landcover values to a c-factor map based on a set of reclass rules
+#%END
+
+#%option
+#% key: inmap
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: input landcover map (integer values aligned with reclass rules)
+#% required : yes
+#%END
+
+#%option
+#% key: outcfact
+#% type: string
+#% gisprompt: string
+#% description: c_factor output map name
+#% answer: year1_cfactor
+#% required : yes
+#%END
+
+#%option
+#% key: cfact_rules
+#% type: string
+#% gisprompt: string
+#% description: path to recode rules file for c-factor map
+#% answer: /usr/local/grass-6.5.svn/scripts/rules/cfactor_recode_rules.txt
+#% required : yes
+#%END
+
+#%option
+#% key: cfact_color
+#% type: string
+#% gisprompt: string
+#% description: path to color rules file for c-factor map
+#% answer: /usr/local/grass-6.5.svn/scripts/rules/cfactor_colors.txt
+#% required : yes
+#%END
+
+import sys
+import os
+import subprocess
+import tempfile
+# first define some useful custom methods
+
+# m is a message (as a string) one wishes to have printed in the output window
+def grass_print(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
+#main block of code starts here
+def main():
+ #setting up variables for use later on
+ inmap = os.getenv("GIS_OPT_inmap")
+ outcfact = os.getenv("GIS_OPT_outcfact")
+ cfact_rules = os.getenv("GIS_OPT_cfact_rules")
+ cfact_color = os.getenv("GIS_OPT_cfact_color")
+ #setting initial conditions of map area
+ grass_com('g.region rast=' + inmap)
+ grass_com('r.mask --quiet input=' + inmap + ' maskcats=*')
+ #creating c-factor map and setting colors
+ grass_com('r.recode --quiet input=' + inmap + ' output=' + outcfact + ' rules=' + cfact_rules)
+ grass_com('r.colors map=' + outcfact + ' rules=' + cfact_color)
+ grass_print('\nCleaning up\n')
+ grass_com('g.remove --quiet rast=MASK')
+ grass_print('\nDONE\n')
+ return
+
+# 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:
+ main()
Property changes on: grass-addons/LandDyn/devs_landcover_scripts/r.cfactor.py/r.cfactor.py
___________________________________________________________________
Added: svn:executable
+ *
Added: grass-addons/LandDyn/devs_landcover_scripts/r.landcover.update.py/r.landcover.update.py
===================================================================
--- grass-addons/LandDyn/devs_landcover_scripts/r.landcover.update.py/r.landcover.update.py (rev 0)
+++ grass-addons/LandDyn/devs_landcover_scripts/r.landcover.update.py/r.landcover.update.py 2010-04-22 20:00:17 UTC (rev 41992)
@@ -0,0 +1,186 @@
+#!/usr/bin/python
+#
+############################################################################
+#
+# MODULE: r.landcover.update
+# AUTHOR(S): Isaac Ullah, Michael Barton, Arizona State University
+# PURPOSE: Updates a map of landcover by the amounts specified in an impacts
+# map
+# ACKNOWLEDGEMENTS: National Science Foundation Grant #BCS0410269
+# COPYRIGHT: (C) 2009 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 a map of landcover by the amounts specified in an impacts map
+#%END
+
+#%option
+#% key: inmap
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: input landcover map (values coded 0-max)
+#% required : yes
+#%END
+
+#%option
+#% key: impacts
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: map of impacts on landcover values
+#% required : yes
+#%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
+#% description: maximum value for landcover maps (number for climax veg)
+#% answer: 50
+#% required : yes
+#%END
+
+#%option
+#% key: outmap
+#% type: string
+#% gisprompt: string
+#% description: land cover output map name (no prefix)
+#% answer: landcover
+#% required : yes
+#%END
+
+#%option
+#% key: lc_rules
+#% type: string
+#% gisprompt: string
+#% description: Path to reclass rules file for making a "labels" map. If no rules specified, no labels map will be made.
+#% answer: /usr/local/grass-6.5.svn/scripts/rules/luse_reclass_rules.txt
+#% required : no
+#%END
+
+#%option
+#% key: lc_color
+#% type: string
+#% gisprompt: string
+#% description: Path to color rules file for landcover map
+#% answer: /usr/local/grass-6.5.svn/scripts/rules/luse_colors.txt
+#% required : no
+#%END
+
+#%flag
+#% key: s
+#% description: -s Output text file of land-use stats from the simulation (will be named "prefix"_luse_stats.txt, and will be overwritten if you run the simulation again with the same prefix)
+#%END
+
+import sys
+import os
+import subprocess
+import tempfile
+# first define some useful custom methods
+
+# m is a message (as a string) one wishes to have printed in the output window
+def grass_print(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_nw(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 charcater seperated list of info to stdout, n is the character that seperates individual pieces of information, and o is a defined blank list to write results to
+def out2list2(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.append(y0)
+ o.append(y1.strip('\n'))
+
+#main block of code starts here
+def main():
+ #setting up variables for use later on
+ inmap = os.getenv('GIS_OPT_inmap')
+ impacts = os.getenv('GIS_OPT_impacts')
+ sfertil = os.getenv('GIS_OPT_sfertil')
+ sdepth = os.getenv('GIS_OPT_sdepth')
+ outmap = os.getenv('GIS_OPT_outmap')
+ max = os.getenv('GIS_OPT_max')
+ lc_rules = os.getenv('GIS_OPT_lc_rules')
+ lc_color = os.getenv('GIS_OPT_lc_color')
+ txtout = outmap + '_landcover_stats.txt'
+ temp_rate = 'temp_rate'
+ reclass_out = outmap + '_labels'
+ #setting initial conditions of map area
+ grass_com('g.region --quiet rast=' + inmap)
+ grass_com('r.mask --quiet input=' + inmap + ' maskcats=*')
+ # calculating rate of regrowth based on current 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
+ grass_com('r.mapcalc "' + temp_rate + '=if(' + sdepth + ' <= 1, ( ( ( (-0.000118528 * (exp(' + sfertil + ',2))) + (0.0215056 * ' + sfertil + ') + 0.0237987 ) + ( ( -0.000118528 * (exp((100*' + sdepth + '),2))) + (0.0215056 * (100*' + sdepth + ')) + 0.0237987 ) ) / 2 ), ( ( ( (-0.000118528 * (exp(' + sfertil + ',2))) + (0.0215056 * ' + sfertil + ') + 0.0237987 ) + 1) / 2 ) )"')
+ #updating raw landscape category numbers based on agent impacts and newly calculated regrowth rate
+ grass_com('r.mapcalc "' + outmap + '=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 + ') ) ) ) )"')
+ grass_com('r.colors --quiet map=' + outmap + ' rules=' + lc_color)
+ try:
+ lc_rules
+ except NameError:
+ lc_rules = None
+ if lc_rules is None:
+ grass_print( "No Labels reclass rules specified, so no Labels map will be made")
+ else:
+ grass_print( 'Creating reclassed Lables map (' + reclass_out +') of text descriptions to raw landscape categories')
+ grass_com('r.reclass --quiet input=' + outmap + ' output=' + reclass_out + ' rules=' + lc_rules)
+ grass_com('r.colors --quiet map=' + reclass_out + ' rules=' + lc_color)
+ #checking total area of updated cells
+ statlist = []
+ out2list2('r.stats -a -n input=' + impacts + ' fs=, nv=* nsteps=1', ',', statlist)
+ grass_print('Total area of impacted zones = ' + statlist[1] + ' square meters\n')
+ #creating optional output text file of stats
+ if os.getenv('GIS_FLAG_s') == '1':
+ f = file(txtout, 'wt')
+ f.write('Stats for ' + prfx + '_landcover\n\nTotal area of impacted zones = ' + statlist[1] + ' square meters\n\n\nLandcover class #, Landcover description, Area (sq. m)\n')
+ p1 = grass_com_nw('r.stats --quiet -a -l -n input=' + prfx +'_landuse1 fs=, nv=* nsteps=' + max )
+ p2 = p1.stdout.readlines()
+ for y in p2:
+ f.write(y)
+ f.close()
+
+ grass_print('\nCleaning up\n')
+ grass_com('g.remove --quiet rast=MASK,' + temp_rate)
+ grass_print("\nDONE!\n")
+ return
+
+# 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:
+ main()
Property changes on: grass-addons/LandDyn/devs_landcover_scripts/r.landcover.update.py/r.landcover.update.py
___________________________________________________________________
Added: svn:executable
+ *
Added: grass-addons/LandDyn/devs_landcover_scripts/r.soil.fertility.py/r.soil.fertility.py
===================================================================
--- grass-addons/LandDyn/devs_landcover_scripts/r.soil.fertility.py/r.soil.fertility.py (rev 0)
+++ grass-addons/LandDyn/devs_landcover_scripts/r.soil.fertility.py/r.soil.fertility.py 2010-04-22 20:00:17 UTC (rev 41992)
@@ -0,0 +1,138 @@
+#!/usr/bin/python
+#
+############################################################################
+#
+# 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: s_fertility
+#% 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.5.svn/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
+
+import sys
+import os
+import subprocess
+import tempfile
+# first define some useful custom methods
+
+# m is a message (as a string) one wishes to have printed in the output window
+def grass_print(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_nw(m):
+ subprocess.Popen('%s' % m, shell='bash')
+ return
+# 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()
+# m is a grass/bash command that will generate some list of info to stdout, n is the character that seperates individual pieces of information, and o is a defined blank list to write results to
+def out2list2(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.append(y0)
+ o.append(y1.strip('\n'))
+
+#main block of code starts here
+def main():
+ #setting up variables for use later on
+ inmap = os.getenv('GIS_OPT_inmap')
+ impacts = os.getenv('GIS_OPT_impacts')
+ outmap = os.getenv('GIS_OPT_outmap')
+ max = os.getenv('GIS_OPT_max')
+ sf_color = os.getenv('GIS_OPT_sf_color')
+ txtout = outmap + '_sfertil_stats.txt'
+ #setting initial conditions of map area
+ grass_com('g.region rast=' + inmap)
+ grass_com('r.mask --quiet -o input=' + inmap + ' maskcats=*')
+ #updating raw soil fertility category numbers
+ grass_com('r.mapcalc "' + outmap + '=if(' + inmap + ' == ' + max + ' && isnull(' + impacts + '), ' + max + ', (if (' + inmap + ' < ' + max + ' && isnull(' + impacts + '), (' + inmap + ' + 1), if (' + inmap + ' > ' + max + ', (' + max + ' - ' + impacts + '), if (' + inmap + ' < 0, 0, (' + inmap + ' - ' + impacts + '))))))"')
+ grass_com('r.colors --quiet map=' + outmap + ' rules=' + sf_color)
+ #checking total area of updated cells
+ temparea = []
+ out2list2('r.stats -a -n input=' + impacts + ' fs=, nv=* nsteps=1', ',', temparea)
+ grass_print('\n\nTotal area of impacted zones = %s square meters\n\n' % temparea[1])
+ #creating optional output text file of stats
+ if os.getenv('GIS_FLAG_s') == '1':
+ f = file(txtout, 'wt')
+ f.write('Stats for ' + outmap+ '\n\nTotal area of impacted zones = ' + temparea[1] + ' square meters\n\nSoil Fertility Value,Area (sq. m)\n\n')
+ p1 = subprocess.Popen('r.stats -A -a -n input=' + outmap + ' fs=, nv=* nsteps=' + max, stdout=subprocess.PIPE, shell='bash')
+ p2 = p1.stdout.readlines()
+ for y in p2:
+ f.write(y)
+ f.close()
+ grass_print('\nCleaning up...\n\n')
+ grass_com('g.remove --quiet rast=MASK')
+ grass_print('DONE!\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.
+if __name__ == "__main__":
+ if ( len(sys.argv) <= 1 or sys.argv[1] != "@ARGS_PARSED@" ):
+ os.execvp("g.parser", [sys.argv[0]] + sys.argv)
+ else:
+ main()
Property changes on: grass-addons/LandDyn/devs_landcover_scripts/r.soil.fertility.py/r.soil.fertility.py
___________________________________________________________________
Added: svn:executable
+ *
Added: grass-addons/LandDyn/devs_landcover_scripts/r.villages.py/r.villages.py
===================================================================
--- grass-addons/LandDyn/devs_landcover_scripts/r.villages.py/r.villages.py (rev 0)
+++ grass-addons/LandDyn/devs_landcover_scripts/r.villages.py/r.villages.py 2010-04-22 20:00:17 UTC (rev 41992)
@@ -0,0 +1,112 @@
+#!/usr/bin/python
+#
+############################################################################
+#
+# MODULE: r.villages
+# AUTHOR(S): Isaac Ullah, Michael Barton, Arizona State University
+# PURPOSE: patches a map of village locations on a landcover landcover map
+# ACKNOWLEDGEMENTS: National Science Foundation Grant #BCS0410269
+# COPYRIGHT: (C) 2009 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: patches a map of village locations on a landcover landcover map
+#%END
+
+#%option
+#% key: inmap
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: input landcover map (integer values aligned with reclass rules)
+#% required : yes
+#%END
+
+#%option
+#% key: villages
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: input map of village locations (coded as all one integer value)
+#% required : yes
+#%END
+
+#%option
+#% key: val
+#% type: integer
+#% answer: 40
+#% description: value for villages to be coded onto output map
+#% required : yes
+#%END
+
+#%option
+#% key: output
+#% type: string
+#% gisprompt: string
+#% description: name for output map
+#% answer: year1
+#% required : yes
+#%END
+
+import sys
+import os
+import subprocess
+import tempfile
+# first define some useful custom methods
+
+# m is a message (as a string) one wishes to have printed in the output window
+def grass_print(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 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()
+
+#main block of code starts here
+def main():
+ #setting up variables for use later on
+ inmap = os.getenv('GIS_OPT_inmap')
+ villages = os.getenv('GIS_OPT_villages')
+ val = os.getenv('GIS_OPT_val')
+ outmap = os.getenv('GIS_OPT_output')
+ temp_villages = outmap + '_temporary_village_reclass'
+ #setting initial conditions of map area
+ grass_com('g.region --quiet rast=' + inmap)
+ grass_com('r.mask --quiet -o input=' + inmap +' maskcats=*')
+ #discovering the value of the input villages map
+ inval = out2var('r.stats -n input=%s fs=space nv=* nsteps=1' % villages)
+ #setting up color and reclass rules
+ colors = tempfile.NamedTemporaryFile()
+ colors.write('%s red' % val)
+ colors.flush()
+ reclass = tempfile.NamedTemporaryFile()
+ reclass.write('%s = %s Village\n' % (inval.strip('\n'), val))
+ reclass.flush()
+ #doing reclass and recolor
+ grass_print('%s = %s Village\n' % (inval.strip('\n'), val))
+ grass_print('r.reclass --quiet input=' + villages + ' output=' + temp_villages + ' rules=%s' % reclass.name)
+ grass_com('r.reclass --quiet input=' + villages + ' output=' + temp_villages + ' rules="%s"' % reclass.name)
+ grass_com('r.colors --quiet map=' + temp_villages + ' rules="%s"' % colors.name)
+ #patching maps together
+ grass_com('r.patch --quiet input=' + temp_villages + ',' + inmap + ' output=' + outmap)
+ grass_print('\nCleaning up...\n')
+ grass_com('g.remove --quiet rast=MASK,' + temp_villages)
+ colors.close()
+ reclass.close()
+ grass_print('\n\nDONE!\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.
+if __name__ == "__main__":
+ if ( len(sys.argv) <= 1 or sys.argv[1] != "@ARGS_PARSED@" ):
+ os.execvp("g.parser", [sys.argv[0]] + sys.argv)
+ else:
+ main()
+
Property changes on: grass-addons/LandDyn/devs_landcover_scripts/r.villages.py/r.villages.py
___________________________________________________________________
Added: svn:executable
+ *
Added: grass-addons/LandDyn/r.accumulate.erdep.py/r.accumulate.erdep.py
===================================================================
--- grass-addons/LandDyn/r.accumulate.erdep.py/r.accumulate.erdep.py (rev 0)
+++ grass-addons/LandDyn/r.accumulate.erdep.py/r.accumulate.erdep.py 2010-04-22 20:00:17 UTC (rev 41992)
@@ -0,0 +1,181 @@
+#!/usr/bin/env python
+
+
+############################################################################
+#
+# MODULE: r.accumulate.erdep.py
+# AUTHOR(S): Isaac Ullah, Arizona State University
+# PURPOSE: Takes the output "netchange" maps and adds them chronologically
+# to make a series showing the "to-date" net change for each year. Good for animations.
+# 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: Takes the output "netchange" maps and adds them chronologically to make a series showing the "to-date" net change for each year. Good for animations.
+#%END
+
+#%option
+#% key: pattern
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: Pattern of first part of file names (prefixes) for map series (leave off #'s, but include any "_"'s or "."'s between the prefix and #)
+#% required : yes
+#%END
+#%option
+#% key: startnum
+#% type: integer
+#% description: Smallest number of the input map series (ie. # of the first map you'd like to include in the cumulative series).
+#% required : yes
+#%END
+#%option
+#% key: endnum
+#% type: integer
+#% description: Largest number of the input map series (ie. # of the last map you'd like to include in the cumulative series).
+#% required : yes
+#%END
+#%option
+#% key: digits
+#% type: integer
+#% description: Total number of digits for input (and output) numbers. (for padding with leading zeros. if zero, numbers are given no leading zeroes)
+#% required : yes
+#% answer: 0
+#%END
+#%option
+#% key: infix
+#% type: string
+#% description: Infix inserted between the prefix and number of the output maps
+#% answer: _cumseries_
+#% required : yes
+#%END
+#%option
+#% key: suffix
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: Pattern of last part of file names (suffixes) for map series with infixed numbers (leave blank if numbers are on the end of the file name!!!)
+#% required : no
+#%END
+#%option
+#% key: statsout
+#% type: string
+#% description: Name of text file to write yearly cumulative erosion/deposition stats to (if blank, no stats file will be written)
+#% answer: erdep_stats.csv
+#% required : no
+#%END
+#%Flag
+#% key: e
+#% description: -e export all maps to PNG files in home directory (good for animation in other programs)
+#%End
+
+import sys
+import os
+import subprocess
+
+def grass_print(m):
+ return subprocess.Popen("g.message message='"'%s'"'" % m, shell='bash').wait()
+
+def main():
+ pattern = os.getenv("GIS_OPT_pattern")
+ startnum = int(os.getenv("GIS_OPT_startnum"))
+ endnum = int(os.getenv("GIS_OPT_endnum"))
+ infix = os.getenv("GIS_OPT_infix")
+ digits = int(os.getenv("GIS_OPT_digits"))
+ #create temp file for color rules
+ f1 = open('temp_color_rules.txt', 'w')
+ f1.write("100% 0 0 100\n5 purple\n0.5 blue\n0.000075 108 166 205\n0 230 230 230\n-0.000075 205 092 092\n-0.5 red\n-5 magenta\n0% 150 0 50")
+ f1.close()
+ #test to see if we make a stats file, and make it if true
+ if bool(os.getenv("GIS_OPT_statsout")) == True:
+ statsfile = file(os.getenv("GIS_OPT_statsout"), 'w')
+ statsfile.write('Erosion Stats,,,,,,Deposition Stats,,,,,\nMax,Min,Mean,Standard Deviation,99th percentile,,Max,Min,Mean,Standard Deviation,99th percentile\n')
+ #if clause tests if numbers are inixed, and then runs the loop accordingly
+ if bool(os.getenv("GIS_OPT_suffix")) == False:
+ grass_print("Numbers are suffixes to prefix: %s" % pattern)
+ for x in range((startnum - 1), endnum):
+ if (x + 1) == startnum:
+ outmap = "%s%s%s" % (pattern, infix, str(startnum).zfill(digits))
+ subprocess.Popen('g.copy --quiet rast=%s%s,%s' % (pattern, str(startnum).zfill(digits), outmap), shell='bash').wait()
+ subprocess.Popen('r.colors --quiet map=%s rules=$HOME/temp_color_rules.txt' % (outmap), shell='bash').wait()
+ else:
+ mapone = "%s%s%s" % (pattern, infix, str(x).zfill(digits))
+ maptwo = "%s%s" % (pattern, str(x + 1).zfill(digits))
+ outmap = "%s%s%s" % (pattern, infix, str(x + 1).zfill(digits))
+ grass_print('doing mapcalc statement for cum netchange map of year %s' % (str(x + 1).zfill(digits)))
+ subprocess.Popen('r.mapcalc "%s = (%s + %s)"' % (outmap, mapone, maptwo), shell='bash').wait()
+ grass_print('setting colors for statement for map %s' % outmap)
+ subprocess.Popen('r.colors --quiet map=%s rules=$HOME/temp_color_rules.txt' % (outmap), shell='bash').wait()
+ if ( os.getenv("GIS_FLAG_e") == "1" ):
+ grass_print('creating png image of map %s' % outmap)
+ subprocess.Popen('r.out.png --quiet input=%s output=%s.png' % (outmap, outmap), shell='bash').wait()
+ if bool(os.getenv("GIS_OPT_statsout")) == True:
+ grass_print('calculating erosion/deposition statistics for map %s' % outmap)
+ subprocess.Popen('r.mapcalc "temperosion=if(%s < 0, %s, null())"' % (outmap, outmap), shell='bash').wait()
+ subprocess.Popen('r.mapcalc "tempdep=if(%s > 0, %s, null())"' % (outmap, outmap), shell='bash').wait()
+ p1 = subprocess.Popen('r.univar -g -e map=temperosion percentile=1', stdout=subprocess.PIPE, shell='bash')
+ erstats = p1.stdout.readlines()
+ dict1 = {}
+ for y in erstats:
+ y0,y1 = y.split('=')
+ dict1[y0] = y1.strip()
+ p2 = subprocess.Popen('r.univar -g -e map=tempdep percentile=99', stdout=subprocess.PIPE, shell='bash')
+ depstats = p2.stdout.readlines()
+ dict2 = {}
+ for z in depstats:
+ z0,z1 = z.split('=')
+ dict2[z0] = z1.strip()
+ statsfile.write('%s,%s,%s,%s,%s,,%s,%s,%s,%s,%s\n' % (dict1['max'], dict1['min'], dict1['mean'], dict1['stddev'], dict1['percentile_1'], dict2['max'], dict2['min'], dict2['mean'], dict2['stddev'], dict2['percentile_99']))
+ else:
+ suffix = os.getenv("GIS_OPT_suffix")
+ grass_print("Numbers are infixes between prefix: %s and suffix: %s" % (pattern, suffix))
+ for x in range((startnum - 1), endnum):
+ if (x + 1) == startnum:
+ outmap = "%s%s%s%s" % (pattern, str(startnum).zfill(digits), infix, suffix)
+ subprocess.Popen('g.copy --quiet rast=%s%s%s,%s' % (pattern, str(startnum).zfill(digits), suffix, outmap), shell='bash').wait()
+ subprocess.Popen('r.colors --quiet map=%s rules=$HOME/temp_color_rules.txt' % (outmap), shell='bash').wait()
+ else:
+ mapone = "%s%s%s%s" % (pattern, str(x).zfill(digits), infix, suffix)
+ maptwo = "%s%s%s" % (pattern, str(x + 1).zfill(digits), suffix)
+ outmap = "%s%s%s%s" % (pattern, str(x + 1).zfill(digits), infix, suffix)
+ grass_print('doing mapcalc statement for cum netchange map of year %s' % (str(x + 1).zfill(digits)))
+ subprocess.Popen('r.mapcalc "%s = (%s + %s)"' % (outmap, mapone, maptwo), shell='bash').wait()
+ grass_print('setting colors for statement for map %s' % outmap)
+ subprocess.Popen('r.colors --quiet map=%s rules=$HOME/temp_color_rules.txt' % (outmap), shell='bash').wait()
+ if ( os.getenv("GIS_FLAG_e") == "1" ):
+ grass_print('creating png image of map %s' % outmap)
+ subprocess.Popen('r.out.png --quiet input=%s output=%s.png' % (outmap, outmap), shell='bash').wait()
+ if bool(os.getenv("GIS_OPT_statsout")) == True:
+ grass_print('calculating erosion/deposition statistics for map %s' % outmap)
+ subprocess.Popen('r.mapcalc "temperosion=if(%s < 0, %s, null())"' % (outmap, outmap), shell='bash').wait()
+ subprocess.Popen('r.mapcalc "tempdep=if(%s > 0, %s, null())"' % (outmap, outmap), shell='bash').wait()
+ p1 = subprocess.Popen('r.univar -g -e map=temperosion percentile=1', stdout=subprocess.PIPE, shell='bash')
+ erstats = p1.stdout.readlines()
+ dict1 = {}
+ for y in erstats:
+ y0,y1 = y.split('=')
+ dict1[y0] = y1.strip('\n')
+ p2 = subprocess.Popen('r.univar -g -e map=tempdep percentile=99', stdout=subprocess.PIPE, shell='bash')
+ depstats = p2.stdout.readlines()
+ dict2 = {}
+ for z in depstats:
+ z0,z1 = z.split('=')
+ dict2[z0] = z1.strip('\n')
+ statsfile.write('%s,%s,%s,%s,%s,,%s,%s,%s,%s,%s\n' % (dict1['max'], dict1['min'], dict1['mean'], dict1['stddev'], dict1['percentile_1'], dict2['max'], dict2['min'], dict2['mean'], dict2['stddev'], dict2['percentile_99']))
+ if bool(os.getenv("GIS_OPT_statsout")) == True:
+ statsfile.close()
+ subprocess.Popen('rm -f $HOME"/temp_color_rules.txt"', shell="bash")
+ subprocess.Popen('g.remove -f rast=temperosion,tempdep', shell="bash")
+ return
+
+if __name__ == "__main__":
+ if ( len(sys.argv) <= 1 or sys.argv[1] != "@ARGS_PARSED@" ):
+ os.execvp("g.parser", [sys.argv[0]] + sys.argv)
+ else:
+ grass_print(" Starting the process--hold on!")
+ grass_print("It is not done until you see DONE WITH EVERYTHING!")
+ main();
+ grass_print("DONE WITH EVERYTHING!")
Property changes on: grass-addons/LandDyn/r.accumulate.erdep.py/r.accumulate.erdep.py
___________________________________________________________________
Added: svn:executable
+ *
Added: grass-addons/LandDyn/r.catchment.py/r.catchment.py
===================================================================
--- grass-addons/LandDyn/r.catchment.py/r.catchment.py (rev 0)
+++ grass-addons/LandDyn/r.catchment.py/r.catchment.py 2010-04-22 20:00:17 UTC (rev 41992)
@@ -0,0 +1,296 @@
+#!/usr/bin/python
+#
+############################################################################
+#
+# MODULE: r.catchment
+# AUTHOR(S): Isaac Ullah, Arizona State University
+# PURPOSE: Creates a raster buffer of specified area around vector points
+# using cost distances. Module requires r.walk.
+# ACKNOWLEDGEMENTS: National Science Foundation Grant #BCS0410269
+# COPYRIGHT: (C) 2009 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: Creates a raster buffer of specified area around vector points using cost distances. Requires r.walk.
+#%END
+
+
+#%option
+#% key: elev
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: Input elevation map (DEM)
+#% required : yes
+#%END
+#%option
+#% key: incost
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: Input cost map (This will override the input elevation map, if none specified, one will be created from input elevation map with r.walk)
+#% required : no
+#%END
+#%option
+#% key: vect
+#% type: string
+#% gisprompt: old,vector,vector
+#% description: Name of input vector site points map
+#% required : yes
+#%END
+#%option
+#% key: frict
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: Optional map of friction costs. If no map selected, default friction=1 making output reflect time costs only
+#% answer:
+#% required : no
+#%END
+#%option
+#% key: a
+#% type: double
+#% description: Coefficients for walking energy formula parameters a,b,c,d
+#% answer: 0.72
+#% required : no
+#%END
+#%option
+#% key: b
+#% type: double
+#% description:
+#% answer: 6.0
+#% required : no
+#%END
+#%option
+#% key: c
+#% type: double
+#% description:
+#% answer: 1.9998
+#% required : no
+#%END
+#%option
+#% key: d
+#% type: double
+#% description:
+#% answer: -1.9998
+#% required : no
+#%END
+#%option
+#% key: lambda
+#% type: double
+#% description: Lambda value for cost distance calculation
+#% answer: 0
+#% required : no
+#%END
+#%option
+#% key: slope_factor
+#% type: double
+#% description: Slope factor determines travel energy cost per height step
+#% answer: -0.2125
+#% required : no
+#%END
+#%option
+#% key: buffer
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: Output buffer map
+#% required : yes
+#%END
+#%option
+#% key: sigma
+#% type: double
+#% description: Slope threshold for mask
+#% answer:
+#% required : no
+#%END
+#%option
+#% key: area
+#% type: integer
+#% description: Area of buffer (Integer value to nearest 100 square map units)
+#% answer: 5000000
+#% required : yes
+#%END
+#%option
+#% key: deviation
+#% type: integer
+#% description: Percentage to which output buffer can differ from desired buffer size (large values decrease run-time, but results will be less precise)
+#% answer: 10
+#% required : yes
+#%END
+#%option
+#% key: mapval
+#% type: integer
+#% description: Integer vlaue for out put catchment area (all other areas will be Null)
+#% answer: 1
+#% required : yes
+#%END
+#%flag
+#% key: k
+#% description: -k Use knight's move for calculating cost surface (slower but more accurate)
+#%END
+#%flag
+#% key: c
+#% description: -c Keep cost surface used to calculate buffers
+#%END
+
+
+import sys
+import os
+import subprocess
+import tempfile
+# first define some useful custom methods
+
+# m is a message (as a string) one wishes to have printed in the output window
+def grass_print(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_nw(m):
+ subprocess.Popen('%s' % m, shell='bash')
+ return
+# 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()
+
+# 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 charcater seperated list of info to stdout, n is the character that seperates individual pieces of information, and o is a defined blank list to write results to
+def out2list2(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.append(y0)
+ o.append(y1.strip('\n'))
+
+#main block of code starts here
+def main():
+ #setting up variables for use later on
+ elev = os.getenv('GIS_OPT_elev')
+ vect = os.getenv('GIS_OPT_vect')
+ lmbda = os.getenv('GIS_OPT_lambda')
+ slpfct = os.getenv('GIS_OPT_slope_factor')
+ a = os.getenv('GIS_OPT_a')
+ b = os.getenv('GIS_OPT_b')
+ c = os.getenv('GIS_OPT_c')
+ d = os.getenv('GIS_OPT_d')
+ sigma = os.getenv('GIS_OPT_sigma')
+ area = float(os.getenv('GIS_OPT_area'))
+ deviation = area*(float(os.getenv('GIS_OPT_deviation'))/100)
+ buffer = os.getenv('GIS_OPT_buffer')
+ mapval = os.getenv('GIS_OPT_mapval')
+
+ grass_print("Wanted buffer area=%s\n" % int(area))
+ grass_com('g.region --quiet rast=' + elev)
+####################################################
+ if bool(os.getenv('GIS_OPT_incost')) is True:
+ grass_print('\n\nstep 1 of 4: Using input cost surface\n')
+ cost = os.getenv('GIS_OPT_incost')
+ else:
+ grass_print('\n\nstep 1 of 4: Calculating cost surface\n')
+ cost = 'temporary.cost'
+ if bool(os.getenv('GIS_OPT_frict')) is True:
+ grass_print('Calculating costs using input friction map\n')
+ frict = os.getenv('GIS_OPT_frict')
+ else:
+ grass_print('Calculating for time costs only')
+ frict = "temporary.friction"
+ grass_com('r.mapcalc "' + frict + '=if(isnull(' + elev + '), null(), 1)"')
+ if os.getenv('GIS_FLAG_k') == '1':
+ grass_print('Using Knight\'s move\n')
+ grass_com('r.walk --quiet -k elevation=' + elev + ' friction=' + frict + ' output=' + cost +' start_points=' + vect + ' lambda=' + lmbda + ' percent_memory=100 nseg=4 walk_coeff=' + a + ',' + b + ',' + c + ',' + d + ' slope_factor=' + slpfct)
+ else:
+ grass_com('r.walk --quiet elevation=' + elev + ' friction=' + frict + ' output=' + cost +' start_points=' + vect + ' lambda=' + lmbda + ' percent_memory=100 nseg=4 walk_coeff=' + a + ',' + b + ',' + c + ',' + d + ' slope_factor=' + slpfct)
+ if bool(os.getenv('GIS_OPT_frict')) is False:
+ grass_com('g.remove --quiet rast=' + frict)
+#################################################
+ grass_print('\n\nstep 2 of 4: Creating optional slope mask\n')
+ if bool(os.getenv('GIS_OPT_sigma')) is True:
+ slope = "temporary.slope"
+ grass_com('r.slope.aspect --quiet --overwrite elevation=' + elev + ' slope=' + slope)
+ grass_com('r.mapcalc "MASK=if(' + slope + ' <= ' + sigma + ', 1, null())"')
+ else:
+ grass_print('no slope mask created')
+##################################################
+ grass_print('\n\nstep 3 of 4: Calculating buffer\n')
+ #set up for loop by getting initial variables from map
+ costdict = {}
+ out2dict('r.univar -g map=%s' % cost,'=' , costdict)
+ maxcost = float(costdict['max'])
+ grass_print("\nMaximum cost distance value= %s\n" % maxcost)
+ arealist = []
+ out2list2('r.stats -a -n input=' + cost + ' fs=, nv=* nsteps=1', ',', arealist)
+ maxarea = int(float(arealist[1]))
+ grass_print('Total map area = %s' % maxarea)
+ cutoff = maxcost
+ mincost = 0.0
+ testarea = maxarea
+ # set up error trap to prevent infinite loops
+ lastcut = 0
+ #start the loop, and home in on the target range
+ while testarea not in range((int(area)-int(deviation)),(int(area)+(int(deviation)+1))) and round(lastcut) != round(cutoff):
+ lastcut = cutoff
+ if testarea > area:
+ maxcost = cutoff
+ cutoff = (maxcost+mincost)/2.0
+ else:
+ mincost = cutoff
+ cutoff = (maxcost+mincost)/2.0
+ temp = tempfile.NamedTemporaryFile()
+ temp.write('0 thru %s = %s\n' % (int(cutoff), mapval))
+ temp.flush()
+ grass_print('0 thru %s = %s\n' % (int(cutoff), mapval))
+ grass_com('r.reclass --overwrite --quiet input=%s output=cost.reclass rules=%s' % (cost, temp.name))
+ temp.close()
+ arealist = []
+ out2list2('r.stats -a -n input=cost.reclass fs=, nv=* nsteps=1', ',', arealist)
+ testarea = int(float(arealist[1]))
+ grass_print('\nTesting with cost distance value = %s......\ncurrent area = %s\n' % (cutoff, testarea))
+ #loop is done, so print the final stats
+ difference = testarea-area
+ grass_print('\n\n*************************\n\n Final cost distance cutoff of %s produces a catchment of %s square meters.\nThe difference from the requested catchment size is %s square meters.' % (cutoff, testarea, difference))
+####################################################
+ grass_print('\n\nstep 4 of 4: Creating output map\n')
+ grass_com('r.mapcalc "%s=if(isnull(cost.reclass), null(), cost.reclass)"' % buffer)
+ grass_print(buffer)
+ grass_com('r.colors --quiet map=%s color=ryb' % buffer)
+ if os.getenv('GIS_FLAG_c') == '1':
+ if bool(os.getenv('GIS_OPT_incost')) is False:
+ grass_com('g.rename temporary.cost,' + buffer + '_cost_surface')
+ grass_print('Cleaning up...(keeping cost map)')
+ grass_com('g.remove --quiet rast=cost.reclass')
+ else:
+ grass_print('Cleaning up...1')
+ grass_com('g.remove --quiet rast=cost.reclass')
+ else:
+ if bool(os.getenv('GIS_OPT_incost')) is False:
+ grass_print('Cleaning up...2')
+ grass_com('g.remove --quiet rast=cost.reclass,temporary.cost' )
+ else:
+ grass_print('Cleaning up...3')
+ grass_com('g.remove --quiet rast=cost.reclass')
+ if bool(os.getenv('GIS_OPT_sigma')) is True:
+ grass_com('g.remove --quiet rast=' + slope)
+ grass_com('g.remove --quiet rast=MASK')
+ grass_print(' DONE!')
+ return
+
+# 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:
+ main()
Property changes on: grass-addons/LandDyn/r.catchment.py/r.catchment.py
___________________________________________________________________
Added: svn:executable
+ *
Added: grass-addons/LandDyn/r.landscape.evol.py/Makefile
===================================================================
--- grass-addons/LandDyn/r.landscape.evol.py/Makefile (rev 0)
+++ grass-addons/LandDyn/r.landscape.evol.py/Makefile 2010-04-22 20:00:17 UTC (rev 41992)
@@ -0,0 +1,7 @@
+MODULE_TOPDIR = ../..
+
+PGM=r.landscape.evol.py
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script
Added: grass-addons/LandDyn/r.landscape.evol.py/r.landscape.evol description.odf_html_585d862d.gif
===================================================================
(Binary files differ)
Property changes on: grass-addons/LandDyn/r.landscape.evol.py/r.landscape.evol description.odf_html_585d862d.gif
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: grass-addons/LandDyn/r.landscape.evol.py/r.landscape.evol description.odf_html_m100fb7e.gif
===================================================================
(Binary files differ)
Property changes on: grass-addons/LandDyn/r.landscape.evol.py/r.landscape.evol description.odf_html_m100fb7e.gif
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: grass-addons/LandDyn/r.landscape.evol.py/r.landscape.evol description.odf_html_m11de82c.gif
===================================================================
(Binary files differ)
Property changes on: grass-addons/LandDyn/r.landscape.evol.py/r.landscape.evol description.odf_html_m11de82c.gif
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: grass-addons/LandDyn/r.landscape.evol.py/r.landscape.evol description.odf_html_m2c6cce6a.gif
===================================================================
(Binary files differ)
Property changes on: grass-addons/LandDyn/r.landscape.evol.py/r.landscape.evol description.odf_html_m2c6cce6a.gif
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: grass-addons/LandDyn/r.landscape.evol.py/r.landscape.evol description.odf_html_m2f9c13ec.gif
===================================================================
(Binary files differ)
Property changes on: grass-addons/LandDyn/r.landscape.evol.py/r.landscape.evol description.odf_html_m2f9c13ec.gif
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: grass-addons/LandDyn/r.landscape.evol.py/r.landscape.evol description.odf_html_m8e0f3ca.gif
===================================================================
(Binary files differ)
Property changes on: grass-addons/LandDyn/r.landscape.evol.py/r.landscape.evol description.odf_html_m8e0f3ca.gif
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: grass-addons/LandDyn/r.landscape.evol.py/r.landscape.evol.py
===================================================================
--- grass-addons/LandDyn/r.landscape.evol.py/r.landscape.evol.py (rev 0)
+++ grass-addons/LandDyn/r.landscape.evol.py/r.landscape.evol.py 2010-04-22 20:00:17 UTC (rev 41992)
@@ -0,0 +1,603 @@
+#!/usr/bin/python
+
+############################################################################
+#
+# MODULE: r.landscape.evol.py
+# AUTHOR(S): iullah
+# COPYRIGHT: (C) 2009 GRASS Development Team/iullah
+#
+# 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! Note that all flags are disabled! Currently flags -z -b -n -f are hard coded to run by default.
+
+# 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
+# the Free Software Foundation; either version 2 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+#############################################################################/
+#%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! Note that all flags are disabled! Currently flags -z -b -n -f are hard coded to run by default.
+#%End
+#%option
+#% key: elev
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: Input elevation map (DEM)
+#% required : yes
+#% guisection: Input
+#%end
+#%option
+#% key: initbdrk
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: Initial bedrock elevations map (for first iteration only)
+#% answer:
+#% required : yes
+#% guisection: Input
+#%end
+#%option
+#% key: prefx
+#% type: string
+#% description: Prefix for all output maps
+#% answer: levol_
+#% required : yes
+#% guisection: Input
+#%end
+#%option
+#% key: outdem
+#% 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
+#%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
+#%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
+#%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: z
+#% description: -z Do not stay zoomed into output maps (return to default region settings)
+#% 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
+#% required : yes
+#% guisection: Variables
+#%end
+#%option
+#% key: K
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: Soil erodability index (K factor) map or constant
+#% answer: 0.42
+#% required : yes
+#% guisection: Variables
+#%end
+#%option
+#% key: sdensity
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: Soil density constant (for conversion from mass to volume)
+#% answer: 1.2184
+#% required : yes
+#% guisection: Variables
+#%end
+#%option
+#% key: C
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: Landcover index (C factor) map or constant
+#% answer: 0.001
+#% 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
+#% required : yes
+#% guisection: Variables
+#%end
+#%option
+#% key: raindays
+#% type: string
+#% description: Number of days of rain per year (integer)
+#% answer: 100
+#% 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
+#% 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
+#% 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
+#% required : yes
+#% guisection: Variables
+#%end
+#%option
+#% key: kappa
+#% type: string
+#% description: Hillslope diffusion (Kappa) rate map or constant (meters per kiloyear)
+#% answer: 1
+#% required : yes
+#% 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 (sq. m uplsope area)
+#% answer: 900
+#% required : yes
+#% guisection: Variables
+#%end
+#%option
+#% key: cutoff2
+#% 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
+#%end
+#%option
+#% key: cutoff3
+#% 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
+#%end
+#%option
+#% key: number
+#% type: integer
+#% 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
+
+#%flag
+#% key: w
+#% description: -w Calcuate for only sheetwash across entire map
+#% guisection: Flow_type
+#%end
+#%flag
+#% key: r
+#% description: -r Calcuate for only channelized flow across entire map
+#% guisection: Flow_type
+#%end
+#%flag
+#% key: d
+#% description: -d Calcuate for only diffusive flow (soil creep) across entire map
+#% guisection: Flow_type
+#%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
+#%end
+
+
+import sys
+import os
+import subprocess
+import tempfile
+
+# first define some useful custom methods
+
+# m is a message (as a string) one wishes to have printed in the output window
+def grass_print(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'))
+
+# 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)
+
+# 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
+ 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")
+ K = os.getenv("GIS_OPT_K")
+ 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)))
+ rain = os.getenv("GIS_OPT_rain")
+ raindays = os.getenv("GIS_OPT_raindays")
+ infilt = os.getenv("GIS_OPT_infilt")
+ 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
+ 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
+ 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 ( 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_com('r.mapcalc "%s=%s - %s"' % (old_soil, old_dem, 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)
+ if ( years == '1' ):
+ slope = '%sslope' % p
+ netchange = '%snetchange' % 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_print('\n##################################################\n\n*************************\n Year %s ' % o + 'step 1 of 7: calculating slope and aspect\n*************************\n')
+ grass_com('r.slope.aspect --quiet elevation=%s slope=%s aspect=%s pcurv=%s tcurv=%s' % (old_dem, slope, aspect, pc, tc))
+ grass_print('\n*************************\n Year %s ' % o + 'step 2 of 7: calculating upslope accumulated flow\n*************************\n')
+ if ( os.getenv("GIS_FLAG_f") == "1" ):
+ grass_print('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_print('Amount of free RAM being allocated for this step: %s Megabites' % mem)
+ grass_com('r.terraflow --q elev=' + old_dem + ' filled=' + p + '.filled direction=' + p + '.direct swatershed=' + p + '.sink accumulation=' + flowacc + ' tci=' + p + '.tci d8cut=infinity memory=' + mem+ ' STREAM_DIR=/var/tmp ')
+ grass_com('g.remove --quiet rast=%s.filled,%s.direct,%s.sink,%s.tci,%s' % (p, p, p, p, tmpdirection))
+ os.remove(os.sep + 'var' + os.sep +'tmp' + os.sep + 'STREAM*')
+ else:
+ grass_print('Using r.watershed to calculate overland flow accumulation per cell (number of cells uplsope from each cell)')
+ grass_com('r.watershed --quiet -fa elevation=' + old_dem + ' accumulation=' + flowacc + ' convergence=5')
+
+ grass_print('\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_print('calculating for sheetwash across the entire map')
+ grass_com('r.mapcalc "%s=%s * sin(%s)"' % (sflowtopo, flowacc, slope))
+ elif ( os.getenv("GIS_FLAG_w") == "0" and os.getenv("GIS_FLAG_d") == "0" and os.getenv("GIS_FLAG_r") == "1" ):
+ grass_print('calculating for channelzed flow across the entire map')
+ grass_com('r.mapcalc "%s=exp((%s),1.6) * exp(sin(%s),1.3)"' % (sflowtopo, flowacc, slope))
+ elif ( os.getenv("GIS_FLAG_w") == "0" and os.getenv("GIS_FLAG_d") == "1" and os.getenv("GIS_FLAG_r") == "0" ):
+ grass_print('calculating for diffusive flow across the entire map')
+ grass_com('r.mapcalc "%s=%s * exp((%s),2) * sin(%s)"' % (sflowtopo, kappa, r, 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)
+ 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 suspended 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_com('r.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 + ') ))))"')
+
+ grass_print('\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_com('r.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 + '))))"')
+
+ grass_com('r.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 + '))))"')
+
+ grass_print('\n*************************\n Year %s ' % o + 'step 5 of 7: calculating partial derivatives for sediment transport\n*************************\n\n')
+ grass_com('r.slope.aspect --q %s dx=%s' % (qsx ,qsxdx))
+ grass_com('r.slope.aspect --q %s dy=%s' % (qsy ,qsydy))
+
+ grass_print('\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 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).
+ grass_com('r.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 + ' ) ) )"')
+
+ grass_print('\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_com('r.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)))"')
+ temp_dem = '%stemp_dem' % p
+ grass_com('r.mapcalc "' + temp_dem + '=' + old_dem + ' + ' + netchange + '"')
+ #do patch-job to catch the shrinking edge problem
+ grass_com('r.patch --quiet input=%s,%s output=%s' % (temp_dem, old_dem, new_dem))
+ #set colors for elevation map to match other dems
+ grass_com('r.colors --q map=' + new_dem + ' rast=' + os.getenv("GIS_OPT_elev"))
+ grass_com('g.remove --quiet rast=%s' % 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_print('\nstep 8.5: Experimental bedrock weathering\n')
+ grass_com('r.mapcalc "' + meancurv + '=((' + pc + ' + ' + tc + ') / 2)"')
+ # create dictionary to record max and min curvature
+ statdict2 = {}
+ out2dict('r.info -r map=%s' % meancurv, '=', statdict2)
+ grass_print('\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_com('r.mapcalc "' + rate + '=' + kappa + '*(2-(' + meancurv + '*(2/(' + statdict2['max'] + ')-(' + 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_com('r.mapcalc "' + new_bdrk + '=' + old_bdrk + ' - ' + rate + '"')
+ grass_print('Calculating new soil depths using new bedrock map')
+ grass_com('r.mapcalc "' + new_soil + '=if ((' + new_dem + ' - ' + new_bdrk + ') < 0, 0, (' + new_dem + ' - ' + 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_com('g.remove --quiet rast=' + rate + ',' + meancurv )
+ else:
+ grass_print('\nCalculating new soil depths keeping bedrock elevations static\n')
+ grass_com('r.mapcalc "' + new_soil + '=if ((' + new_dem + ' - ' + initbdrk + ') < 0, 0, (' + new_dem + ' - ' + initbdrk + '))"')
+ #setting colors for new soil depth map
+ grass_print('reading color rules from %s' % sdcolors.name)
+ grass_print (sdcolors.readlines())
+ grass_com('r.colors --quiet map=%s rules=%s' % (new_soil , 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_com('r.mapcalc "' + tmperosion + '=if(' + netchange + ' < 0, ' + netchange + ', null())"')
+ grass_com('r.mapcalc "' + tmpdep + '=if(' + netchange + ' > 0, ' + netchange + ', null())"')
+ #grab the stats fromt hese temp files and save them to dictionaries
+ soilstats = {}
+ out2dict('r.univar -g -e map=' + new_soil + ' percentile=99', '=', soilstats)
+ erosstats = {}
+ out2dict('r.univar -g -e map=' + tmperosion + ' percentile=1', '=', erosstats)
+ depostats = {}
+ out2dict('r.univar -g -e map=' + tmpdep + ' percentile=99', '=', depostats)
+ #clena up temp maps
+ grass_com('g.remove --quiet rast=' + tmperosion + ',' + tmpdep)
+ #write stats to a new line in the stats file
+ grass_print('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 foroutput netchange map
+ if ( os.getenv("GIS_FLAG_n") == "1" ):
+ grass_print('Not keeping a netchange map')
+ grass_com('g.remove --quiet rast=' + netchange)
+ else:
+ grass_com('r.colors --quiet map=%s rules=%s' % (netchange, nccolors.name))
+ sdcolors.close()
+ nccolors.close()
+ grass_print('\n*************************\nDone with Year %s ' % o + '\n*************************\n')
+ if ( m == '0' and os.getenv("GIS_FlAG_e") == "1" ):
+ grass_print('\nkeeping initial soil depths map ' + old_soil + '\n')
+ elif m == '0':
+ grass_com('g.remove --quiet rast=' + old_soil)
+ grass_print('\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.
+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
+ years = os.getenv("GIS_OPT_number")
+ prefx = os.getenv("GIS_OPT_prefx")
+ #make the stats out file with correct column headers
+ if os.getenv("GIS_OPT_statsout") == "":
+ proc = subprocess.Popen("g.gisenv get=MAPSET", stdout=subprocess.PIPE, shell='bash')
+ out = proc.stdout.read()
+ mapset = out.strip()
+ 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')
+ #let's grab the current resolution
+ reg1 = {}
+ out2dict('g.region -p -m -g', '=', reg1)
+ # we must set the region to the map being used. otherwise r.flow will not work
+ grass_print('Setting region to %s' % os.getenv("GIS_OPT_elev"))
+ reg2 = {}
+ out2dict('g.region -a -p -m -g rast=%s' % os.getenv("GIS_OPT_elev"), '=', reg2)
+ #print resolution before and after
+ grass_print('Old resolution was %s' % reg1['nsres'])
+ grass_print('New resolution is %s' % reg2['nsres'])
+ grass_print('\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_print("Total number of iterations to be run is %s years" % years)
+ # This is the loop!
+ for x in range(int(years)):
+ grass_print("Iteration = %s" % (x + 1))
+ main(x, (x + 1), prefx, statsout, reg2['nsres']);
+ #since we are now done with the loop, close the stats file.
+ f.close()
+ #reset the region if we were asked to
+ if os.getenv("GIS_FLAG_z") == "1":
+ grass_print('\nIterations complete, restoring default region settings\n')
+ grass_com('g.region -d -g')
+ else:
+ grass_print('\nIterations complete, keeping region set to output maps\n')
+ if os.getenv("GIS_FLAG_k") == "1":
+ grass_print('\nTemporary maps have NOT been deleted\n')
+ else:
+ grass_print('\nCleaning up temporary maps...\n\n')
+ if os.getenv("GIS_FLAG_p") == "0":
+ grass_com('g.remove -f --quiet rast=%sslope*' % prefx)
+ grass_com('g.mremove -f --quiet rast=%sflowacc*' % prefx)
+ grass_com('g.mremove -f --quiet rast=%saspect*' % prefx)
+ grass_com('g.mremove -f --quiet rast=%ssflowtopo*' % prefx)
+ grass_com('g.mremove -f --quiet rast=%sqsx*' % prefx)
+ grass_com('g.mremove -f --quiet rast=%sqsy*' % prefx)
+ grass_com('g.mremove -f --quiet rast=%serosdep*' % prefx)
+ grass_com('g.mremove -f --quiet rast=%spc*' % prefx)
+ grass_com('g.mremove -f --quiet rast=%stc*' % prefx)
+ grass_print('\nDone\n\n')
+ grass_print("\n\nDone with everything")
+
+
+
Property changes on: grass-addons/LandDyn/r.landscape.evol.py/r.landscape.evol.py
___________________________________________________________________
Added: svn:executable
+ *
Added: grass-addons/LandDyn/r.landscape.evol.py/r.landscape.evol.py.html
===================================================================
--- grass-addons/LandDyn/r.landscape.evol.py/r.landscape.evol.py.html (rev 0)
+++ grass-addons/LandDyn/r.landscape.evol.py/r.landscape.evol.py.html 2010-04-22 20:00:17 UTC (rev 41992)
@@ -0,0 +1,341 @@
+<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.0 Transitional//EN">
+<HTML>
+<HEAD>
+ <META HTTP-EQUIV="CONTENT-TYPE" CONTENT="text/html; charset=utf-8">
+ <TITLE></TITLE>
+ <META NAME="GENERATOR" CONTENT="OpenOffice.org 3.1 (Unix)">
+ <META NAME="CREATED" CONTENT="0;0">
+ <META NAME="CHANGEDBY" CONTENT="Isaac Ullah">
+ <META NAME="CHANGED" CONTENT="20100422;11051600">
+</HEAD>
+<BODY LANG="en-US" DIR="LTR">
+<H2>DESCRIPTION</H2>
+<P><EM>r.landscape.evol</EM> takes as input a raster digital
+elevation model (DEM) of surface topography and an input raster DEM
+of bedrock elevations, as well as several environmental variables,
+and computes the net change in elevation due to erosion and
+deposition on the hill-slopes using the USPED equation, and in the
+stream channels using a process equation based on the excess stream
+power or shear stress. The module has the ability to run recursively,
+looping over several iterations. The time interval represented by
+each iteration is determined by the scale of the input environmental
+variables, and as such, all input variables should be on the same
+time scale. The script creates a new map where each raster cell
+carries a numerical value, which represents the simulated meters of
+erosion or deposition (ED) estimated for that cell, under the
+specified conditions of rainfall intensity, soil erodability, water
+flow, and vegetation cover. This map of net ED is then added to (for
+deposition) or subtracted from (for erosion) the topography map of
+the previous time step, to create a new topography map (i.e., as a
+DEM) after a cycle of landuse and landscape change.</P>
+<P><B>R</B>, <B>K</B>, and <B>C</B> are environmental factors in the
+USPED equation that relate to the intensity of yearly rainfall, the
+erodability of soil, and the degree to which vegetation cover
+prevents erosion (See below for a detailed description of these
+factors). These factors largely determine the amount of erosion or
+deposition that occur on the hill-slopes. <B>cutoff1</B>, <B>cutoff2,
+</B><SPAN STYLE="font-weight: normal">and </SPAN><B>cutoff3</B> are
+values of flow accumulation (amount of upslope area in square meters)
+that determine where surface processes change from soil-creep to
+laminar overland flow (sheetwash), from laminar overland flow to
+channelized overland flow (rills/gullies), and from channelized
+overland flow to full stream flow respectively. <B>kappa</B> is the
+rate of diffusion for soil-creep in meters per 1000 years. <B>sdensity</B>
+is the density of the soil in grams per cubic centimeters. <B>rain</B><SPAN STYLE="font-weight: normal">
+is the total annual precipitation measured in meters (or the average
+annual rainfall in meters per year). </SPAN><B>raindays</B><SPAN STYLE="font-weight: normal">
+is the total number of days on which it rained in one year (or an
+average value of days per year). </SPAN><B>infilt</B><SPAN STYLE="font-weight: normal">
+is the proportion of rainfall that infiltrates into the soil and thus
+does not contribute to runoff (values are between 0 and 1). </SPAN><B>Kt</B><SPAN STYLE="font-weight: normal">
+is the stream transport efficiency variable that describes the
+cohesivness of the stream channel beds (0.001 for normal
+gravel/sandy/silt channel bed to 0.000001 for a bedrock channel bed).
+</SPAN><B>loadexp</B><SPAN STYLE="font-weight: normal"> is the stream
+transport type variable that determines the type of stream transport
+modeled (1.5 for bedload transport, or 2.5 for suspended load
+transport). </SPAN><B>alpha</B><SPAN STYLE="font-weight: normal"> is
+the critical slope threshold above which the model will simulate the
+cumulative effects of mass wasting (landsliding). These</SPAN>
+measures all need to be determined empirically for a given landscape
+under a given climatic condition, but the defaults are average values
+for the Circum-Mediterranean Basin.
+</P>
+<P>By default, <EM>r.watershed</EM> is used to calculate flow
+accumulation modeling using the MFD alglrithm included in GRASS 6.4
+and higher. This can be made backwards compatable by checking the -f
+flag, which will use <I>r.terraflow </I><SPAN STYLE="font-style: normal">to
+compute a flow accumulation model using the SFD algorithm. This will,
+however, produce much less accurate results, and users are therefore
+encouraged to used GRASS 6.4 or higher.</SPAN></P>
+<P> The user may use the <B>statsout</B> option to define the name of
+the file that contains the statistics of erosion, deposition, and
+soil depths over all iterations. The default name is
+<TT>"mapset"_"prefix"_lsevol_stats.txt</TT> (in
+the users home directory).
+</P>
+<H2>CALCULATING SURFACE EROSION AND DEPOSITION</H2>
+<P>Because physical laws that govern the flow of water across
+landscapes and its ability to erode, entrain, transport, and deposit
+sediments can be expressed in mathematical form, they can be
+translated into a scripting algorithm that modifies raster landscapes
+(i.e., in the GIS) in ways analogous to the ways in which real-world
+landscapes change. There are various mathematical expressions of the
+relevant surface processes in the geomorphological literature
+depending for example on the processes selected to be represented,
+the simplicity of representation desired, and the degree of
+resolution desired (Clevis, et al. 2006; Degani, et al. 1979; Mitas
+and Mitasova 1998; Mitasova, Hofierka, et al. 1996; Mitasova and
+Mitas 2001a, b; Peeters, et al. 2006; Singh and Phadke 2006; Warren,
+et al. 2005; Wischmeier, et al. 1971; Wischmeier and Smith 1978). We
+use the Unit Stream Power Erosion-Deposition (USPED) equation,
+derived in part from the widely-used Revised Universal Soil Loss
+Equation (RUSLE) (American Society of Agricultural Engineers 2003;
+Degani, et al. 1979; Mitasova, et al. 2001; Mitasova, Mitas, et al.
+1996; Mitasova, et al. 2004; Singh and Phadke 2006; Warren, et al.
+2005; Wischmeier 1976; Wischmeier, et al. 1971; Wischmeier and Smith
+1978), to calculate net erosion and deposiiton across each landscape
+cell above the flow accumualtion breakpoint <B>cutoff3</B>. USPED was
+developed for hillslopes, small watersheds, and small channels (i.e.,
+rills and gullies) (Warren, et al. 2005), and is less applicable to
+larger streams and rivers. Therefore we use a different process
+equation to model erosion and deposition in stream channels (see
+below).
+</P>
+<P>Net erosion and deposition rates on hillslopes are computed from
+the change in sediment flow across cells of a DEM that have flow
+accumulation values less than <B>cutoff3</B>. We approximate sediment
+flow rate from sediment transport capacity, assuming that water
+flowing over landscapes normally carries sediment at capacity.
+Transport capacity is calculated by combining a rainfall coefficient
+(R, MJ mm/ha h yr), soil erodability coefficient (K, Mg ha h/ha MJ
+mm), and coefficient for the ability of vegetation to prevent erosion
+(C, unitless) from RUSLE with with an estimate of topographically
+driven stream power as shown in equation (1)</P>
+<P><IMG SRC="r.landscape.evol_html_m11de82c.gif" NAME="Object2" ALIGN=ABSMIDDLE HSPACE=8 WIDTH=157 HEIGHT=21></P>
+<P>where <I>A</I> is the upslope contributing area (a measure of
+water flowing through a cell) and <EM>B</EM> is the slope of the
+cell. The exponents <EM>m</EM> and <EM>n</EM> are empirically derived
+and vary for water flowing over nearly level ground, on hillslopes,
+in water catchments at the heads of gullies, or in small channels.
+The sediment flow rate is largely determined by the amount of water
+flowing (contributing area), its velocity (a function of slope), the
+erodability of the substrate (K factor), and the ability of the
+vegetation cover to prevent erosion (C factor).</P>
+<P>Implementing the USPED algorithm in a GRASS script combines GIS
+modules for calculating slope, aspect, and flow accumulation (the
+amount of water that flows across each cell) using map algebra. Data
+used by the script includes a map of initial surface topography (a
+raster DEM), soil erodability (a constant for uniform soil or a
+raster map for variable soil), vegetation cover (a constant or raster
+map), and rainfall intensity (a constant only). We also create an
+underlying bedrock topography map (a raster DEM) to limit the total
+depth of unconsolidated sediment that can be eroded. Soil
+erodability, vegetation cover, and rainfall are expressed as the
+K-factor <I>(K),</I> C-factor (<I>C</I><SPAN STYLE="font-style: normal">)</SPAN>,
+and R-factor (<I>R</I>)<SPAN STYLE="font-style: normal"> components</SPAN>
+of the RUSLE and have been calculated empirically for a variety of
+setting (Boellstorff and Benito 2005; MartÃnez-Casasnovas, 2000;
+Essa 2004; Hammad, et al. 2004; Renard, et al. 1997; Renard and
+Freimund 1994).
+</P>
+<P>For areas of the DEM that have flow accumulation values greater
+than <B>cutoff3 </B><SPAN STYLE="font-weight: normal">(ie. areas
+that are proper streams), we use a case of the transport limited
+process law that is formulated for water flowing in stream channels
+(Howard 1980; Tucker and Hancock 2010). This is done by first
+calculating the reach average shear stress (</SPAN><FONT FACE="Times New Roman, serif"><SPAN STYLE="font-weight: normal">τ</SPAN></FONT><SPAN STYLE="font-weight: normal">),
+here estimated for a cellular landscape simply as:</SPAN></P>
+<P><IMG SRC="r.landscape.evol_html_m2f9c13ec.gif" NAME="Object1" ALIGN=ABSMIDDLE HSPACE=8 WIDTH=119 HEIGHT=22></P>
+<P> <SPAN STYLE="font-weight: normal">Where: </SPAN><I><SPAN STYLE="font-weight: normal">9806.65</SPAN></I><SPAN STYLE="font-weight: normal">
+is a constant related to the gravitational acceleration of water, </SPAN><I><SPAN STYLE="font-weight: normal">B</SPAN></I><SPAN STYLE="font-weight: normal">
+is the slope of the cell in degrees, and </SPAN><I><SPAN STYLE="font-weight: normal">D</SPAN></I><SPAN STYLE="font-weight: normal">
+is the instantaneous depth of flowing water in the cell. </SPAN><I><SPAN STYLE="font-weight: normal">D
+</SPAN></I><SPAN STYLE="font-style: normal"><SPAN STYLE="font-weight: normal">is</SPAN></SPAN><SPAN STYLE="font-weight: normal">
+here assumed to be roughly equivalent to the depth of flow during the
+average minute of rainfall, calculated by:</SPAN></P>
+<P STYLE="font-weight: normal"><IMG SRC="r.landscape.evol_html_m2c6cce6a.gif" NAME="Object3" ALIGN=ABSMIDDLE HSPACE=8 WIDTH=137 HEIGHT=42></P>
+<P><SPAN STYLE="font-weight: normal">Where: </SPAN><I><SPAN STYLE="font-weight: normal">R</SPAN></I><SUB><I><SPAN STYLE="font-weight: normal">m</SPAN></I></SUB><SPAN STYLE="font-weight: normal">
+is the total annual precipitation in meters, </SPAN><I><SPAN STYLE="font-weight: normal">i</SPAN></I><SPAN STYLE="font-weight: normal">
+is the proportion of rainfall that infiltrates rather than </SPAN><SPAN STYLE="font-weight: normal">runs
+off, </SPAN><I><SPAN STYLE="font-weight: normal">A</SPAN></I><SPAN STYLE="font-style: normal"><SPAN STYLE="font-weight: normal">
+is the uplsope accumulated area per unit contour width at the cell,
+</SPAN></SPAN><I><SPAN STYLE="font-weight: normal">R</SPAN></I><SUB><I><SPAN STYLE="font-weight: normal">d</SPAN></I></SUB><SPAN STYLE="font-style: normal"><SPAN STYLE="font-weight: normal">
+is the number of days on which it rained in a one year period, and
+</SPAN></SPAN><I><SPAN STYLE="font-weight: normal">1440</SPAN></I><SPAN STYLE="font-style: normal"><SPAN STYLE="font-weight: normal">
+is a constant relating to the number of minutes in a day.</SPAN></SPAN></P>
+<P STYLE="font-style: normal; font-weight: normal">Then the transport
+capacity is calculated by:</P>
+<P STYLE="font-style: normal; font-weight: normal"><IMG SRC="r.landscape.evol_html_m100fb7e.gif" NAME="Object4" ALIGN=ABSMIDDLE HSPACE=8 WIDTH=76 HEIGHT=28></P>
+<P><SPAN STYLE="font-weight: normal">Where: </SPAN><I><SPAN STYLE="font-weight: normal">K</SPAN></I><SUB><I><SPAN STYLE="font-weight: normal">t</SPAN></I></SUB><SPAN STYLE="font-style: normal"><SPAN STYLE="font-weight: normal">
+is the transport efficiency factor related to the character of the
+stream bed (0.001 for normal sediment to 0.000001 for bedrock), and </SPAN></SPAN><I><SPAN STYLE="font-weight: normal">n</SPAN></I><SPAN STYLE="font-style: normal"><SPAN STYLE="font-weight: normal">
+is an empirically determined exponent related to the dominant type of
+transport in the stream system (1.5 for bedload transport or 2.5
+suspended load transport).</SPAN></SPAN></P>
+<P>Net erosion and deposition rates are then computed across the
+entire DEM as change in sediment flow in the x and y directions
+across a cell as follows”</P>
+<P><IMG SRC="r.landscape.evol_html_m8e0f3ca.gif" NAME="Object6" ALIGN=ABSMIDDLE HSPACE=8 WIDTH=204 HEIGHT=38></P>
+<P><SPAN STYLE="font-weight: normal">where ED is net erosion or
+deposition rate for sediment and </SPAN><EM><FONT FACE="Times New Roman, serif"><SPAN STYLE="font-weight: normal">α</SPAN></FONT></EM><SPAN STYLE="font-weight: normal">
+is the topographic aspect (i.e., direction of slope) for a cell.
+Whether flowing water will erode or deposit sediment in a particular
+cell is determined by the </SPAN><EM><SPAN STYLE="font-style: normal"><U><SPAN STYLE="font-weight: normal">change</SPAN></U></SPAN></EM><SPAN STYLE="font-weight: normal">
+in sediment flow (transport capacity) from one cell to the next. If
+the transport capacity increases (for example, due to an increase in
+the steepness of the slope or amount of flowing water), more sediment
+will be entrained and erosion will occur; if the transport capacity
+decreases (for example, due to a decrease in slope or water flow)
+sediment will be deposited.</SPAN></P>
+<P>The output of this GRASS implementation of these transport
+equations must be modified in several ways in order to make it
+appropriate for landscape evolution simulation. First, because of the
+way slope is calculated in <EM>r.slope.aspect</EM>, the flux <I>T</I>
+is actually calculated one cell downslope from where is really
+occurs. This causes problems when USPED is iterated over many cycles,
+and creates oscillating "spikes" in positive and negative
+flux values resulting in the calculation of alternating deep pits and
+high mounds at sensitive areas on the landscape. To overcome this,
+<EM>r.landscape.evol</EM> uses a nieghborhood algorithm in <EM>r.mapcalc</EM>
+to put the calculated value of <I>T</I> back into the cell that is
+most uplsope from where it is originally calculated.
+</P>
+<P>Additionally, control must be kept for the amount of erodible
+sediment available to moved. <EM>r.landscape.evol</EM> explicitly
+tracks this by taking the difference between the input bedrcok
+elevation DEM, and the current surface topography DEM, and creating a
+map of "soil" depth. This map tracks the amount of material
+assumed to be available for entrainment and transport by surface
+processes. A simple logical algorithm is used to prevent unduly large
+amounts of erosion from being calculated in areas devoid of erodible
+materials (ie. at bedrock outcrops). Where this condition occurs, <I>K</I>
+or <I>K</I><SUB><I>t </I></SUB>is made to be very small, resulting in
+only extremely small amounts of erosion.
+</P>
+<P>Another major issue is that the total flux <I>T </I>is in units of
+Tons/Ha, which means it must be converted in order to calculate the
+change in elevation at each cell (<I>m</I><SUB><I>vert</I></SUB>).
+This is done via a simple algorithm that uses the density of the soil
+and the cell resolution:</P>
+<P><IMG SRC="r.landscape.evol_html_585d862d.gif" NAME="Object5" ALIGN=ABSMIDDLE HSPACE=8 WIDTH=174 HEIGHT=20></P>
+<P>Where: <I>10000</I> is the number of meters per hectare, <I>Sd </I>is
+the density of the soil, and <I>Res </I>is the cell resolution
+(width). In order to convert the output back to Tons/Ha (standard
+rate for USPED/RUSLE equations), you can multiply the <B>netchange</B>
+output map by "(10000 x resolution x soil density)" to
+create a map of soil erosion/deposition rates across the landscape.
+</P>
+<H2>SEE ALSO</H2>
+<UL>
+ <LI><P STYLE="margin-bottom: 0in">The <A HREF="http://medland.asu.edu/">MEDLAND</A>
+ project at Arizona State University
+ </P>
+ <LI><P><A HREF="r.watershed.html">r.watershed</A>, <A HREF="r.terraflow.html">r.terraflow</A>,
+ <A HREF="r.mapcalc.html">r.mapcalc</A>
+ </P>
+</UL>
+<H2>REFERENCES</H2>
+<P>American Society of Agricultural Engineers 2003 Honoring the
+Universal Soil Loss Equation: Historic Landmark Dedication Pamphlet.
+Purdue University Department of Agricultural and Biological
+Engineering.
+</P>
+<P>Clevis, Q., G. E. Tucker, G. Lock, S. T. Lancaster, N. Gasparini,
+A. Desitter and R. L. Bras 2006 Geoarchaeological simulation of
+meandering river deposits and settlement distributions: A
+three-dimensional approach. Geoarchaeology 21(8):843-874.
+</P>
+<P>Degani, A., L. A. Lewis and B. B. Downing 1979 Interactive
+Computer Simulation of the Spatial Process of Soil Erosion.
+Professional Geographer 31(2):184-190.
+</P>
+<P STYLE="margin-left: 0.5in; text-indent: -0.5in">Howard, A. D.
+1980. Thresholds in river regimes. <SPAN STYLE="font-style: normal">Thresholds
+in geomorphology</SPAN>, 227–258.
+</P>
+<P>Mitas, L. and H. Mitasova 1998 Distributed soil erosion simulation
+for effective erosion prevention. Water Resources Research
+34(3):505-516.
+</P>
+<P>Mitasova, H., J. Hofierka, M. Zlocha and L. R. Iverson 1996
+Modelling topographic potential for erosion and deposition using GIS.
+International Journal of Geographical Information Systems
+10(5):629-641.
+</P>
+<P>Mitasova, H. and L. Mitas 2001a Modeling Physical Systems. In
+Geographic Information Systems and Environmental Modeling, edited by
+B. O. Parks, M. Crane and K. C. Clarke, pp. 189-210. Prentice Hall,
+New York. 2001b Multiscale soil erosion simulations for land use
+management. In Landscape erosion and landscape evolution modeling,
+edited by R. Harmon and W. Doe, pp. 321-347. Kluwer Academic/Plenum
+Publishers, New York.
+</P>
+<P>Mitasova, H., L. Mitas and W. M. Brown 2001 Multiscale simulation
+of land use impact on soil erosion and deposition patterns. In
+Sustaining the Global Farm. Selected Papers from the 10th
+International Soil Conservation Organization Meeting, May 1999,
+Purdue University, edited by D. E. Stott, R. H. Mohtar and G. C.
+Steinhardt, pp. 1163-1169. USDA-ARS National Soil Erosion Research
+Laboratory, Purdue.
+</P>
+<P>Mitasova, H., L. Mitas, W. M. Brown and D. Johnston 1996
+Multidimensional Soil Erosion/Deposition Modeling Part III: Process
+based erosion simulation. Geographic Modeling and Systems Laboratory,
+University of Illinois at Urban-Champaign.
+</P>
+<P>Mitasova, H., C. Thaxton, J. Hofierka, R. McLaughlin, A. Moore and
+M. L 2004 Path sampling method for modeling overland water flow,
+sediment transport and short term terrain evolution in Open Source
+GIS. In Proceedings of the XVth International Conference on
+Computational Methods in Water Resources (CMWR XV), edited by C. T.
+Miller, M. W. Farthing, V. G. Gray and G. F. Pinder, pp. 1479-1490.
+Elsevier, Chapel Hill, NC, USA.
+</P>
+<P>Peeters, I., T. Rommens, G. Verstraeten, G. Govers, A. Van
+Rompaey, J. Poesen and K. Van Oost 2006 Reconstructing ancient
+topography through erosion modelling. Geomorphology 78(3-4):250-264.
+</P>
+<P>Rawls, W. J. 1983 Estimating soil bulk denisty from particle size
+analysis and organic matter content. Soil Science 135(2):123.
+</P>
+<P>Renard, K. G., G. R. Foster, G. A. Weesies, D. K. McCool and D. C.
+Yoder 1997 Predicting soil erosion by water: a guide to conservation
+planning with the Revised Universal Soil Loss Equation (RUSLE). In
+Agriculture Handbook, pp. 1–251. vol. 703. US Department of
+Agriculture, Washington, DC.
+</P>
+<P>Renard, K. G. and J. R. Freimund 1994 Using monthly precipitation
+data to estimate the R-factor in the revised USLE. Journal of
+Hydrology 157(1-4):287-306.
+</P>
+<P>Singh, R. and V. S. Phadke 2006 Assessing soil loss by water
+erosion in Jamni River Basin, Bundelkhand region, India, adopting
+universal soil loss equation using GIS. Current Science
+90(10):1431-1435.
+</P>
+<P>Tucker, G. E. and G. R Hancock 2010 Modelling landscape
+evolution. <SPAN STYLE="font-style: normal">Earth Surface Processes
+and Landforms</SPAN> 35(1): 28–50.
+</P>
+<P>Warren, S. D., H. Mitasova, M. G. Hohmann, S. Landsberger, F. Y.
+Iskander, T. S. Ruzycki and G. M. Senseman 2005 Validation of a 3-D
+enhancement of the Universal Soil Loss Equation for prediction of
+soil erosion and sediment deposition. Catena 64:281-296.
+</P>
+<P>Wischmeier, W. H. 1976 Use and Misuse of the Universal Soil Loss
+Equation. Journal of Soil and Water Conservation 31:5-9.
+</P>
+<P>Wischmeier, W. H., C. B. Johnson and B. V. Cross 1971 A Soil
+Erodibility Nomograph for Farmland and Construction Sites. Journal of
+Soil and Water Conservation 26:189-92.
+</P>
+<P>Wischmeier, W. H. and D. D. Smith 1978 Predicting Rainfall-Erosion
+Losses - A Guide to Conservation Planning. USDA Agriculture Handbook
+282.
+</P>
+<P><BR><BR>
+</P>
+<P><I>Last changed: $Date: 2009-23-1 (Fri, 21 Jan 2009) $</I></P>
+</BODY>
+</HTML>
\ No newline at end of file
Added: grass-addons/LandDyn/r.soildepth.py/r.soildepth.py
===================================================================
--- grass-addons/LandDyn/r.soildepth.py/r.soildepth.py (rev 0)
+++ grass-addons/LandDyn/r.soildepth.py/r.soildepth.py 2010-04-22 20:00:17 UTC (rev 41992)
@@ -0,0 +1,159 @@
+#!/usr/bin/python
+#
+############################################################################
+#
+# MODULE: r.soildepth
+# AUTHOR(S): Isaac Ullah, Arizona State University
+# PURPOSE: Create soil depth map based on hillslope curvature.
+# ACKNOWLEDGEMENTS: National Science Foundation Grant #BCS0410269
+# Based on equations from Heimsath et al, 1997
+# 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: Create soil depth map based on hillslope curvature
+#%End
+#%option
+#% key: elev
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: Input elevation map (DEM)
+#% required : yes
+#%END
+#%option
+#% key: bedrock
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: Output bedrock elevation map
+#% required : yes
+#%END
+#%option
+#% key: soildepth
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: Output soil depths map
+#% required : no
+#%END
+#%option
+#% key: min
+#% type: double
+#% description: What are the actual (empircal) minimum soil depths in this area (in meters)?
+#% answer: 0.0001
+#% required : yes
+#%END
+#%option
+#% key: max
+#% type: double
+#% description: What are the actual (empircal) maximum soil depths in this area (in meters)?
+#% answer: 5.0
+#% required : yes
+#%END
+#%flag
+#% key: k
+#% description: -k Keep the soil rate map (map name will be the same as specified in input option 'bedrock" with suffix "_rate" appended to it)
+#%END
+#%flag
+#% key: s
+#% description: -s Print some useful statistics about the output soil depths (written to stdout)
+#%END
+
+import sys
+import os
+import subprocess
+import tempfile
+
+# m is a message (as a string) one wishes to have printed in the output window
+def grass_print(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 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 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()
+
+def main():
+ if os.getenv('GIS_OPT_soildepth') == None:
+ soildepth = "temp.soildepth"
+ else:
+ soildepth = os.getenv('GIS_OPT_soildepth')
+ elev = os.getenv('GIS_OPT_elev')
+ bedrock = os.getenv('GIS_OPT_bedrock')
+ min = os.getenv('GIS_OPT_min')
+ max = os.getenv('GIS_OPT_max')
+ slope = "temp_slope_deletable"
+ pc = "temp_pc_deletable"
+ tc = "temp_tc_deletable"
+ temprate = "temp_rate_deletable"
+ #let's grab the current resolution
+ reg1 = {}
+ out2dict('g.region -p -m -g', '=', reg1)
+ res = int(float(reg1['nsres']))
+ # 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()
+ grass_print('STEP 1, calculating profile and tangental curvatures\n')
+ grass_com('r.slope.aspect --quiet --overwrite format=percent elevation=%s slope=%s pcurv=%s tcurv=%s' % (elev, slope, pc, tc))
+ grass_print('STEP 2, Calculating soil depth ratios across the landscape\n')
+ grass_com('r.mapcalc "' + temprate + '=eval(x = (-1*(' + pc + '+' + tc + ')), y = (if(' + slope + ' < 1, 1, ' + slope + ')), ((50*x)+50)/y)"')
+ # create dictionary to record max and min rate so we can rescale it according the user supplied max and min desired soil depths
+ ratedict = {}
+ out2dict('r.info -r map=%s' % temprate, '=', ratedict)
+ grass_print('STEP 3, Calculating actual soil depths across the landscape (based on user input min and max soil depth values)\n')
+ #creating and running a linear regression to scale the calculated landform soil depth ratios into real soil dpeths using user specified min and max soil depth values
+ #grass_com('r.mapcalc "' + soildepth + '=eval(x = ((' + ratedict['max'] + '-' + max +')/(' + ratedict['min'] + '-' + min + ')), y = (' + ratedict['max'] + '-(y*' + ratedict['min'] + ')), (y*' + temprate + ')+x)"') #####this is the discreet way using mapcalc, I thik the r.recode way (below) may be faster
+ recode = tempfile.NamedTemporaryFile()
+ recode.write('%s:%s:%s:%s' % (ratedict['min'] , ratedict['max'] , min, max))
+ recode.flush()
+ grass_com('r.recode --quiet input=%s output=%s rules=%s' % (temprate, soildepth, recode.name))
+ recode.close()
+ #grab some stats if asked to
+ if os.getenv('GIS_FLAG_s') == '1':
+ depthdict = {}
+ out2dict('r.univar -ge map=' + soildepth + ' percentile=90', '=', depthdict)
+ grass_print('STEP 4, calculating bedrock elevation map\n')
+ grass_com('r.mapcalc "' + bedrock + '=(' + elev + ' - ' + soildepth + ')"')
+ grass_print('Cleaning up...')
+ if os.getenv('GIS_FLAG_k') == '1':
+ grass_com('g.remove --quiet rast=' + pc + ',' + tc)
+ grass_com('g.rename --quiet rast=' + temprate + ',' + bedrock + '_rate')
+ else:
+ grass_com('g.remove --quiet rast=' + pc + ',' + tc + ',' + temprate)
+ if os.getenv('GIS_OPT_soildepth') is None:
+ grass_com('g.remove --quiet rast=' + soildepth)
+ else:
+ grass_com('r.colors --quiet map=%s rules=%s' % (soildepth , sdcolors.name))
+ grass_print('\nDONE!\n')
+ if os.getenv('GIS_FLAG_s') == '1':
+ for key in depthdict.keys():
+ grass_print('%s=%s' % (key, depthdict[key]))
+ grass_print('Total volume of soil is %s cubic meters' % (float(depthdict['sum'])*res))
+ return
+
+# 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:
+ main()
+
+
+
+
+
Property changes on: grass-addons/LandDyn/r.soildepth.py/r.soildepth.py
___________________________________________________________________
Added: svn:executable
+ *
More information about the grass-commit
mailing list