[GRASS-SVN] r64601 - in grass-addons/grass7/raster: r.landscape.evol r.viewshed.cva
svn_grass at osgeo.org
svn_grass at osgeo.org
Fri Feb 13 01:49:03 PST 2015
Author: neteler
Date: 2015-02-13 01:49:02 -0800 (Fri, 13 Feb 2015)
New Revision: 64601
Added:
grass-addons/grass7/raster/r.landscape.evol/r.landscape.evol.html
grass-addons/grass7/raster/r.landscape.evol/r.landscape.evol.py
grass-addons/grass7/raster/r.landscape.evol/r.landscape.evol_description.odf.odt
grass-addons/grass7/raster/r.landscape.evol/r.landscape.evol_description.odt
grass-addons/grass7/raster/r.viewshed.cva/r.viewshed.cva.py
Removed:
grass-addons/grass7/raster/r.landscape.evol/description.html
grass-addons/grass7/raster/r.landscape.evol/r.landscape.evol
grass-addons/grass7/raster/r.landscape.evol/r.landscape.evol description.odf.odt
grass-addons/grass7/raster/r.landscape.evol/r.landscape.evol description.odt
grass-addons/grass7/raster/r.viewshed.cva/r.viewshed.cva
Log:
Addons: Python scripts need .py extension for G7; fix white space in file names
Deleted: grass-addons/grass7/raster/r.landscape.evol/description.html
===================================================================
--- grass-addons/grass7/raster/r.landscape.evol/description.html 2015-02-13 03:10:26 UTC (rev 64600)
+++ grass-addons/grass7/raster/r.landscape.evol/description.html 2015-02-13 09:49:02 UTC (rev 64601)
@@ -1,366 +0,0 @@
-<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 either 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><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>and <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. Note that some
-experimentation is required in order to find the best possible values
-for these cutoffs, and the -p flag will provide some output data that
-may be useful for this. <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>
-is the total annual precipitation measured in meters (or the average
-annual rainfall in meters per year). <b>raindays</b>
-is the total number of days on which it rained in one year (or an
-average value of days per year). <b>infilt</b>
-is the proportion of rainfall that infiltrates into the soil and thus
-does not contribute to runoff (values are between 0 and 1). <b>Kt</b>
-is the stream transport efficiency variable that describes the
-cohesiveness of the stream channel beds (0.001 for normal
-gravel/sandy/silt channel bed to 0.000001 for a bedrock channel bed).
-<b>loadexp</b> 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). <b>alpha</b> is
-the critical slope threshold above which the model will simulate the
-cumulative effects of mass wasting (landsliding). These
-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>By default, <em>r.watershed</em> is used to calculate flow
-accumulation modeling using the MFD algorithm included in GRASS 6.4
-and higher. This can be made backwards compatible by checking the -f
-flag, which will use <i>r.terraflow </i>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.
-<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).
-
-<h3>Calculating Erosion and Deposition</h3>
-<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>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)
-<center>
-<img src="r_landscape_evol_equation1.gif"><br>
-</center>
-<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>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>),
-and R-factor (<i>R</i>) components
-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>For areas of the DEM that have flow accumulation values greater
-than <b>cutoff3 </b>(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 (<FONT FACE="Times New Roman, serif">τ</FONT>),
-here estimated for a cellular landscape simply as:
-<center>
-<p><img src="r_landscape_evol_equation2.gif"><br>
-</center>
-<p> Where: <i>9806.65</i>
-is a constant related to the gravitational acceleration of water, <i>B</i>
-is the slope of the cell in degrees, and <i>D</i>
-is the instantaneous depth of flowing water in the cell. <i>D
-</i>is
-here assumed to be roughly equivalent to the depth of flow during the
-average minute of rainfall, calculated by:
-<center>
-<img src="r_landscape_evol_equation3.gif"><br>
-</center>
-<p>Where: <i>R</i><sub><i>m</i></sub>
-is the total annual precipitation in meters, <i>i</i>
-is the proportion of rainfall that infiltrates rather than runs
-off, <i>A</i>
-is the uplsope accumulated area per unit contour width at the cell,
-<i>R</i><sub><i>d</i></sub>
-is the number of days on which it rained in a one year period, and
-<i>1440</i>
-is a constant relating to the number of minutes in a day.
-<p>Then the transport capacity is calculated by:
-<center>
-<img src="r_landscape_evol_equation4.gif"><br>
-</center>
-<p>Where: <i>K</i><sub><i>t</i></sub>
-is the transport efficiency factor related to the character of the
-stream bed (0.001 for normal sediment to 0.000001 for bedrock), and <i>n</i>
-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).
-<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:
-<center>
-<img src="r_landscape_evol_equation5.gif"><br>
-</center>
-<p>where ED is net erosion or
-deposition rate for sediment and <em><FONT FACE="Times New Roman, serif">α</FONT></em>
-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 <em><U>change</U></em>
-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.
-<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 neighborhood 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>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 bedrock
-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>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:
-<center>
-<img src="r_landscape_evol_equation6.gif"><br>
-</center>
-<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 raster cell resolution x <b>sdensity</b>)" to
-create a map of soil erosion/deposition rates across the landscape.
-
-<h3>Determining Cutoff Points</h3>
-<p>
-To get started with <em>r.landscape.evol</em>, you need to determine the appropriate values for <b>cutoff1</b>, <b>cutoff2</b>, and <b>cutoff3</b>, which are transition points between different types of erosive processes. These are in units of flow accumulation scaled to actual surface flow as determined in r.watershed from the values of rainfall and flow hindrance from vegetation. To do this, you should parameterize the module as best as possible, EXCEPT for the three "cutoffs". Then, run the module with the <b>-p</b> flag, which will make a random points vector file with the values of scaled flow accumulation (scaled to actual rainfall and vegetation), profile curvature, and tangential curvature in the associated table. Plotting the log of the scaled flow accumulation against each of these two curvatures will help you to determine reasonable values for the cutoffs, as each transition should show a unique relationship between curvature and flow accumulations. See
the figures below for examples:
-
-<center>
-<img src="r_landscape_evol_Flow_acc_vs_curvature.png" width="1000" height="500" alt="Log Scaled Flow Accumulation versus Topographic Curvatures"><br>
-
-Log Scaled Flow Accumulation versus Topographic Curvatures.<br><br>
-
-<img src="r_landscape_evol_Map1.png" width="500" height="284" alt="Map showing the spatial patterns of the cutoffs determined from the previous figure"><br>
-
-Map showing the spatial patterns of the cutoffs determined from the previous figure.<br><br>
-</center>
-
-<p>
-<h3>Note About Climate Parameters</h3>
-<p>
-r.landscape.evol accepts an external "climate file", which should be a comma separated plain text file with four columns in the order of, "<b>rain</b>,<b>R</b>,<b>storms</b>,<b>stormlength</b>" (without headers). Each of these columns must exist, although there need not be values in every column (i.e., you can enter a single value for any of these parameters in the command line, and combine that with populated columns for the other values). Note that the climate file must have the same number of rows as there are iterations of the simulation (<b>years</b>).
-
-<h2>SEE ALSO</h2>
-<ul>
- <li><p>The <a href="http://medland.asu.edu/">MEDLAND</a>
- project at Arizona State University
-
- <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>
-
- <li><p>Mitasova, H., C. M. Barton, I. I. Ullah, J. Hofierka, and R. S. Harmon 2013 GIS-based soil erosion modeling. In Remote Sensing and GIScience in Geomorphology, edited by J. Shroder and M. P. Bishop. 3:228-258. San Diego: Academic Press.
-
-
-</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>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>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>Howard, A. D. 1980. Thresholds in river regimes. Thresholds
-in geomorphology, 227-258.
-
-<p>Mitas, L. and H. Mitasova 1998 Distributed soil erosion simulation
-for effective erosion prevention. Water Resources Research
-34(3):505-516.
-
-<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>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>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>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>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>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>Rawls, W. J. 1983 Estimating soil bulk denisty from particle size
-analysis and organic matter content. Soil Science 135(2):123.
-
-<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-51. vol. 703. US Department of
-Agriculture, Washington, DC.
-
-<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>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>Tucker, G. E. and G. R Hancock 2010 Modelling landscape
-evolution. Earth Surface Processes
-and Landforms 35(1): 28-50.
-
-<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>Wischmeier, W. H. 1976 Use and Misuse of the Universal Soil Loss
-Equation. Journal of Soil and Water Conservation 31:5-9.
-
-<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>Wischmeier, W. H. and D. D. Smith 1978 Predicting Rainfall-Erosion
-Losses - A Guide to Conservation Planning. USDA Agriculture Handbook
-282.
-
-
-<h2>AUTHORS</h2>
-Isaac I. Ullah, C. Michael Barton, and Helena Mitasova
-
-<p><i>Last changed: $Date: 2014-010-02 05:03:58 -0400 (Thur, 02 Oct 2014) $</i>
Deleted: grass-addons/grass7/raster/r.landscape.evol/r.landscape.evol
===================================================================
--- grass-addons/grass7/raster/r.landscape.evol/r.landscape.evol 2015-02-13 03:10:26 UTC (rev 64600)
+++ grass-addons/grass7/raster/r.landscape.evol/r.landscape.evol 2015-02-13 09:49:02 UTC (rev 64601)
@@ -1,687 +0,0 @@
-#!/usr/bin/python
-
-############################################################################
-#
-# MODULE: r.landscape.evol.py
-# AUTHOR(S): Isaac Ullah and Michael Barton
-# COPYRIGHT: (C) 2012 GRASS Development Team/Isaac Ullah
-#
-# description: Simulates the cumulative effect of erosion and deposition on a landscape over time. This module uses appropriate flow on different landforms by default; however, singular flow regimes can be chosen by manipulating the cutoff points. This module requires GRASS 6.4 or greater. THIS SCRIPT WILL PRODUCE MANY TEMPORARY MAPS AND REQUIRES A LOT OF FREE FILE SPACE!
-
-# This program is free software; you can redistribute it and/or modify
-# it under the terms of the GNU General Public License as published by
-# 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: Simulates the cumulative effect of erosion and deposition on a landscape over time. This module uses appropriate flow on different landforms by default; however, singular flow regimes can be chosen by manipulating the cutoff points. This module requires GRASS 6.4 or greater. THIS SCRIPT WILL PRODUCE MANY TEMPORARY MAPS AND REQUIRES A LOT OF FREE FILE SPACE!
-#%End
-#%option
-#% key: elev
-#% type: string
-#% gisprompt: old,cell,raster
-#% description: Input elevation map (DEM of surface)
-#% required : yes
-#%end
-#%option
-#% key: initbdrk
-#% type: string
-#% gisprompt: old,cell,raster
-#% description: Bedrock elevations map (DEM of bedrock)
-#% answer:
-#% required : yes
-#%end
-#%option
-#% key: k
-#% type: string
-#% gisprompt: old,cell,raster
-#% description: Soil erodability index (K factor) map or constant
-#% answer: 0.42
-#% required : no
-#% guisection: Landscape Evolution
-#%end
-#%option
-#% key: sdensity
-#% type: string
-#% gisprompt: old,cell,raster
-#% description: Soil density map or constant [T/m3] for conversion from mass to volume
-#% answer: 1.2184
-#% required : no
-#% guisection: Landscape Evolution
-#%end
-#%option
-#% key: kt
-#% type: double
-#% description: Stream transport efficiency variable (0.001 for a soft substrate, 0.0001 for a normal substrate, 0.00001 for a hard substrate, 0.000001 for a very hard substrate)
-#% answer: 0.0001
-#% required : no
-#% options : 0.001,0.0001,0.00001,0.000001
-#% guisection: Landscape Evolution
-#%end
-#%option
-#% key: loadexp
-#% type: double
-#% description: Stream transport type variable (1.5 for mainly bedload transport, 2.5 for mainly suspended load transport)
-#% answer: 1.5
-#% options: 1.5,2.5
-#% required : no
-#% guisection: Landscape Evolution
-#%end
-#%option
-#% key: kappa
-#% type: double
-#% description: Hillslope diffusion (Kappa) rate map or constant [m/kyr]
-#% answer: 1
-#% required : no
-#% guisection: Landscape Evolution
-#%end
-#%option
-#% key: c
-#% type: string
-#% gisprompt: old,cell,raster
-#% description: Landcover index (C factor) map or constant
-#% answer: 0.005
-#% required : no
-#% guisection: Landscape Evolution
-#%end
-#%option
-#% key: rain
-#% type: string
-#% gisprompt: old,cell,raster
-#% description: Precip totals for the average storm [mm] (or path to climate file of comma separated values of "rain,R,storms,stormlength", with a new line for each year of the simulation)
-#% answer: 20.61
-#% guisection: Climate
-#%end
-#%option
-#% key: r
-#% type: string
-#% description: Rainfall (R factor) constant (AVERAGE FOR WHOLE MAP AREA) (or path to climate file of comma separated values of "rain,R,storms,stormlength", with a new line for each year of the simulation)
-#% answer: 4.54
-#% guisection: Climate
-#%end
-#%option
-#% key: storms
-#% type: string
-#% description: Number of storms per year (integer) (or path to climate file of comma separated values of "rain,R,storms,stormlength", with a new line for each year of the simulation)
-#% answer: 25
-#% guisection: Climate
-#%end
-#%option
-#% key: stormlength
-#% type: string
-#% description: Average length of the storm [h] (or path to climate file of comma separated values of "rain,R,storms,stormlength", with a new line for each year of the simulation)
-#% answer: 24.0
-#% guisection: Climate
-#%end
-#%option
-#% key: speed
-#% type: double
-#% description: Average velocity of flowing water in the drainage [m/s]
-#% answer: 1.4
-#% required : no
-#% guisection: Hydrology
-#%end
-#%option
-#% key: manningn
-#% type: string
-#% gisprompt: old,cell,raster
-#% description: Map or constant of the value of Manning's "N" value for channelized flow.
-#% answer: 0.05
-#% required : no
-#% guisection: Hydrology
-#%end
-#%option
-#% key: flowcontrib
-#% type: string
-#% gisprompt: old,cell,raster
-#% description: Map or constant indicating how much each cell contributes to downstream flow (as a "percentage" from 0-100). If no map or value entered, routine will assume 100% downstream contribution
-#% required : no
-#% guisection: Hydrology
-#%end
-#%option
-#% key: convergence
-#% type: integer
-#% description: Value for the flow convergence variable in r.watershed. Small values make water spread out, high values make it converge in narrower channels.
-#% answer: 5
-#% options: 1,2,3,4,5,6,7,8,9,10
-#% required : no
-#% guisection: Hydrology
-#%end
-#%option
-#% key: cutoff1
-#% type: double
-#% description: Flow accumulation breakpoint value for shift from diffusion to overland flow
-#% answer: 0
-#% required : no
-#% guisection: Hydrology
-#%end
-#%option
-#% key: cutoff2
-#% type: double
-#% description: Flow accumulation breakpoint value for shift from overland flow to rill/gully flow (if value is the same as cutoff1, no sheetwash procesess will be modeled)
-#% answer: 100
-#% required : no
-#% guisection: Hydrology
-#%end
-#%option
-#% key: cutoff3
-#% type: double
-#% description: Flow accumulation breakpoint value for shift from rill/gully flow to stream flow (if value is the same as cutoff2, no rill procesess will be modeled)
-#% answer: 100
-#% required : no
-#% guisection: Hydrology
-#%end
-#%option
-#% key: smoothing
-#% type: string
-#% description: Amount of additional smoothing (answer "no" unless you notice large spikes in the erdep rate map)
-#% answer: no
-#% options: no,low,high
-#% required : yes
-#%end
-#%option
-#% key: prefx
-#% type: string
-#% description: Prefix for all output maps
-#% answer: levol_
-#% required : yes
-#%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
-#%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
-#%end
-#%option
-#% key: number
-#% type: integer
-#% description: Number of iterations (cycles) to run
-#% answer: 1
-#% required : yes
-#%end
-
-
-# #%option
-# #% key: alpha ###This may be added back in if I can find a good equation for mass movement
-# #% type: integer
-# #% description: Critical slope threshold for mass movement of sediment (in degrees above horizontal)
-# #% answer: 40
-# #% required : yes
-# #%end
-
-#%flag
-#% key: p
-#% description: -p Output a vector points map with sampled values of flow accumulation and curvatures suitable for determining cutoff values. NOTE: Overrides all other output options, and exits after completion. The output vector points map will be named "PREFIX_#_randomly_sampled_points".
-#% guisection: Optional
-#%end
-#%flag
-#% key: 1
-#% description: -1 Calculate streams as 1D difference instead of 2D divergence
-#% guisection: Landscape Evolution
-#%end
-#%flag
-#% key: c
-#% description: -c Calculate streams with a shear stress equation, rather than a stream-power equation
-#% guisection: Landscape Evolution
-#%end
-#%flag
-#% key: k
-#% description: -k Keep ALL temporary maps (overides flags -drst). This will make A LOT of maps!
-#% guisection: Optional
-#%end
-#%flag
-#% key: d
-#% description: -d Don't output yearly soil depth maps
-#% guisection: Optional
-#%end
-#%flag
-#% key: r
-#% description: -r Don't output yearly maps of the erosion/deposition rates ("ED_rate" map, in vertical meters)
-#% guisection: Optional
-#%end
-#%flag
-#% key: s
-#% description: -s Keep all slope maps
-#% guisection: Optional
-#%end
-#%flag
-#% key: t
-#% description: -t Keep yearly maps of the Transport Capacity at each cell ("Qs" maps)
-#% guisection: Optional
-#%end
-#%flag
-#% key: e
-#% description: -e Keep yearly maps of the Excess Transport Capacity (divergence) at each cell ("DeltaQs" maps)
-#% guisection: Optional
-#%end
-#%Option
-#% key: statsout
-#% type: string
-#% description: Name for the statsout text file (optional, if none provided, a default name will be used)
-#% required: no
-#% guisection: Optional
-#%end
-
-import sys
-import os
-import math
-import tempfile
-grass_install_tree = os.getenv('GISBASE')
-sys.path.append(grass_install_tree + os.sep + 'etc' + os.sep + 'python')
-import grass.script as grass
-
-# 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, s = master list of lists of climate data
-def main(m, o, p, q, r, s):
- #get the process id to tag any temporary maps we make for easy clean up in the loop
- pid = os.getpid()
- # Get variables from user input
- smoothing = options["smoothing"]
- years = options["number"]
- initbdrk = options["initbdrk"]
- outdem = options["outdem"]
- outsoil = options["outsoil"]
- K = options["k"]
- sdensity = options["sdensity"]
- C = options["c"]
- kappa = options["kappa"]
- cutoff1 = options["cutoff1"]
- cutoff2 = options["cutoff2"]
- cutoff3 = options["cutoff3"]
- flowcontrib = options["flowcontrib"]
- convergence = options["convergence"]
- manningn = options["manningn"]
- old_bdrk = options["initbdrk"]
- # these variables come in as a list of lists, so let's get this year's numbers out of them.
- rain = s[0][m]
- storms = s[2][m]
- #CHANGES
- R = s[1][m] / storms
- #R = s[1][m]
- stormlengthsecs = float(s[3][m])*3600.00 # number of seconds in the storm
- stormtimet = stormlengthsecs / (float(options["speed"]) * float(r)) # number of hydrologic instants in the storm
-# timet = stormlengthsecs/stormtimet # length of a single hydrologic instant in seconds, currently unused, but might be important in future versions
- Kt = options["kt"]
- loadexp = options["loadexp"]
- # Make some variables for temporary map names, labeled different depending on if we keep them or not
- if ( flags["k"] is True ):
- aspect = '%saspect%04d' % (p, o)
- flowacc = '%sflowacc%04d' % (p, o)
- flacclargenums = '%sflowacc_largenums%04d' % (p, o)
- flowdir = '%sflowdir%04d' % (p, o)
- pc = '%spc%04d' % (p, o)
- tc = '%stc%04d' % (p, o)
-# meancurv = '%smeancurv%04d' % (p, o) # This variable might be used if bedrock weathering is ever implemented
-# rate = '%srate%04d' % (p, o) # This variable might be used if bedrock weathering is ever implemented
- rainexcess = "%s_rainfall_excess_map_%04d"% (p, o)
- tempnetchange1 = '%sTEMPORARY_unsmoothed_ED_rate%04d' % (p, o)
- tempnetchange2 = '%sTEMPORARY_smoothed_ED_rate%04d' % (p, o)
- tmperosion = '%sTEMPORARY_erosion%04d' % (p, o)
- tmpdep = '%sTEMPORARY_deposition%04d' % (p, o)
- else:
- aspect = '%saspect%04d' % (pid, o)
- flowacc = '%sflowacc%04d' % (pid, o)
- flacclargenums = '%sflowacc_largenums%04d' % (pid, o)
- flowdir = '%sflowdir%04d' % (pid, o)
- pc = '%spc%04d' % (pid, o)
- tc = '%stc%04d' % (pid, o)
-# meancurv = '%smeancurv%04d' % (pid, o) # This variable might be used if bedrock weathering is ever implemented
-# rate = '%srate%04d' % (p, o) # This variable might be used if bedrock weathering is ever implemented
- rainexcess = "%s_rainfall_excess_map_%04d"% (pid, o)
- tempnetchange1 = '%sTEMPORARY_unsmoothed_ED_rate%04d' % (pid, o)
- tempnetchange2 = '%sTEMPORARY_smoothed_ED_rate%04d' % (pid, o)
- tmperosion = '%sTEMPORARY_erosion%04d' % (pid, o)
- tmpdep = '%sTEMPORARY_deposition%04d' % (pid, 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' % options["elev"]
- old_soil = "%s%s_init" % (prefx, options["outsoil"])
- grass.mapcalc('${old_soil}=${old_dem}-${old_bdrk}', overwrite = "True", quiet = "True", old_soil = old_soil, old_dem = old_dem, old_bdrk = old_bdrk)
- else :
- old_dem = '%s%s%04d' % (p, options["outdem"], m)
- old_soil = '%s%s%04d' % (p, options["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 = '%sED_rate' % p
- new_dem ='%s%s' % (p, outdem)
- new_soil = '%s%s' % (p, outsoil)
- else:
- slope = '%sslope%04d' % (p, o)
- netchange = '%sED_rate%04d' % (p, o)
- new_dem = '%s%s%04d' % (p, outdem, o)
- new_soil = '%s%s%04d' % (p, outsoil, o)
- #Check to see if we are going to only output diagnostics for determing cutoff values, and act accordingly
- if ( flags["p"] is True ):
- grass.message('GATHERING STATISTICS FOR DETERMINING CUTOFF VALUES\n-------------------------------------------------\n1) Calculating slope and curvatures')
- grass.run_command('r.slope.aspect', quiet = "True", elevation = old_dem, slope = slope, pcurv = pc, tcurv = tc)
- else:
- grass.message('\n##################################################\n\n*************************\n Iteration %s -- ' % o + 'step 1: calculating slope\n*************************\n')
- grass.run_command('r.slope.aspect', quiet = "True", elevation = old_dem, aspect = aspect, slope = slope)
- if ( flags["p"] is True ):
- grass.message('2) Calculating map of rainfall excess')
- else:
- grass.message('\n*************************\n Iteration %s -- ' % o + 'step 2: calculating accumulated flow depths\n*************************\n')
- grass.message('Calculating runoff excess rates (uplsope accumulated cells scaled to \"flowcontrib\" map')
-
- #make map of rainfall excess (proportion each cell contributes to downstrem flow) from flowcontrib. Note that if flowcontrib is a map, we are just making a copy of it. This map is a percentage, but has to be scale from 0-100, because r.watershed will only allow values greater than 1 as input in it's 'flow' variable. This creates a flow accumulation map with large numbers, but this map will be divided by 100 after it is made, which brings the values back down to what they should be.
- if flowcontrib == "":
- flowcontrib = 100
- grass.mapcalc('${rainexcess}=int(${flowcontrib})', quiet = "True", rainexcess = rainexcess, flowcontrib = flowcontrib)
- if ( os.getenv("GIS_FLAG_p") == "1" ):
- grass.message('3) Calculating accumulated flow (in numbers of upslope cells, scaled by runoff contribution')
- grass.run_command('r.watershed', quiet = "True", flags = 'a', elevation = old_dem, flow = rainexcess, accumulation = flacclargenums, drainage = flowdir, convergence = convergence)
- grass.mapcalc('${flowacc}=${flacclargenums}/100', quiet = "True", flowacc = flowacc, flacclargenums = flacclargenums)
- #again, do something different if we are only making an evaluation of cutoffs
- if ( flags["p"] is True ):
- grass.message('4) Determining number of sampling points using formula: "ln(#cells_in_input_map)*100"')
- flaccstats = grass.parse_command('r.univar', flags = 'g', map = flowacc)
- numpts = int(math.log(int(flaccstats['n']))*100)
- grass.message('5) Creating random points and sampling values of flow accumulation, curvatures, and slope.')
- vout = '%s%s_randomly_sampled_points' % (p, numpts)
- grass.run_command('r.random', quiet = "True", input = flowacc, cover = pc, n = numpts, vector_output = vout)
- grass.run_command('v.db.renamecol', quiet = "True", map = vout, column = 'value,Flow_acc')
- grass.run_command('v.db.renamecol', quiet = "True", map = vout, column = 'covervalue,Princ_curv')
- grass.run_command('v.db.addcol', quiet = "True", map = vout, columns = 'Tang_curv double precision, Slope double precision')
- grass.run_command('v.what.rast', quiet = "True", vector = vout, raster = tc, column = "Tang_curv")
- grass.run_command('v.what.rast', quiet = "True", vector = vout, raster = slope, column = "Slope")
- if ( flags["k"] is True ):
- grass.message('--Keeping the created maps (Flow Accumulation, Slope, Principle Curvature, Tangential Curvature)')
- else:
- grass.message('6) Cleaning up...')
- grass.run_command('g.remove', quiet = "True", flags = 'f', type = "rast", name = slope + "," + pc + "," + tc + "," + flowacc)
- grass.message('FINISHED. \nRandom sample points map "%s" created successfully.\n' % vout)
- sys.exit(0)
- grass.message('\n*************************\n Iteration %s -- ' % o + 'step 3: calculating sediment transport rates (units variable depending upon process) \n*************************\n')
- # This step calculates the force of the flowing water at every cell on the landscape using the proper transport process law for the specific point in the flow regime. For upper hillslopes (below cutoff point 1) this done by multiplying the diffusion coeficient by the accumulated flow/cell res width. For midslopes (between cutoff 1 and 2) this is done by multiplying slope by accumulated flow with the m and n exponents set to 1. For channel catchment heads (between cutoff 2 and 3), this is done by multiplying slope by accumulated flow with the m and n exponents set to 1.6 and 1.3 respectively. For Channelized flow in streams (above cutoff 3), this is done by calculating the reach average shear stress (hydraulic radius [here estimated for a cellular landscape simply as the depth of flow] times slope times accumulated flow [cells] times gravitatiopnal acceleration of water [9806.65 newtons], all raised to the appropriate exponant for the type of transport (bedload or suspe
nded load), and then divided by the resolution. Depth of flow is calculated as a mean "instantaneous depth" during any given rain event, here estimated by the maximum depth of an idealized unit hydrograph with base equal to the duration of the storm, and area equal to the total accumulated excess rainfall during the storm. Then finally calculates the stream power or sediment carrying capacity (qs) of the water flowing at each part of the map by multiplying the reach average shear stress (channelized flow in streams) or the estimated flow force (overland flow) by the transport coeficient (estimated by R*K*C for hillslopes or kt for streams). This is a "Transport Limited" equation, however, we add some constraints on detachment by checking to see if the sediment supply has been exhausted: if the current soil depth is 0 or negative (checking for a negative value is kind of an error trap) then we make the transport coefficient small (0.000001) to simulate erosion on bedrock. Bec
ause diffusion and USPED require 2D divergence later on, we calculate these as vectors in the X and Y directions. Stream flow only needs 1D difference, so it's calulated in the direction of flow.
- if ( flags["1"] is True ):
- #This is the version with 1D streams
- qs1 = '%sQs_1D_streams%04d' % (p, o)
- #CHANGES
- #choose shear stress or stream power
- #these are the stream-power versions * note that I'm converting the stream-power output (kg/m2) to same units as USPED (T/ha) by multiplying by ten. This ensures they are even going into the divergence calculation #Qs = Kt * n^-1 * 9810 * depth^1.6 * tan(slope)^1.5
- if ( flags['c'] is True ):
- grass.mapcalc('${qs1}=(${Kt} * exp(9806.65*(((${rain}/1000)*${flowacc})/(0.595*${stormtimet}))*tan(${slope}), ${loadexp}) )', quiet = "True", qs1 = qs1, flowacc = flowacc, stormtimet = stormtimet, rain = rain, slope = slope, loadexp = loadexp, Kt = Kt, sdensity = sdensity)
- else:
- grass.mapcalc('${qs1}=10 * ${Kt} * exp(${manningn}, -1) * 9810 * exp( ( ( (${rain}/1000)*${flowacc}) / (0.595*${stormtimet}) ), 1.6) * exp(tan(${slope}), 1.5)', quiet = "True", qs1 = qs1, flowacc = flowacc, stormtimet = stormtimet, rain = rain, slope = slope, loadexp = loadexp, Kt = Kt, sdensity = sdensity, manningn = manningn)
- qsx = "%sQsx_%04d" % (p,o)
- qsy = "%sQsy_%04d" % (p,o)
- grass.mapcalc("${qsx}=eval(a=(${kappa} * sin(${slope}) * cos(${aspect})), b=((${R}*${K}*${C}*${flowacc}*${res}*sin(${slope})) * cos(${aspect})), c=( (${R}*${K}*${C}*exp((${flowacc}*${res}),1.6000000)*exp(sin(${slope}),1.3000000)) * cos(${aspect})), if(${flowacc} <= ${cutoff1}, a, if(${flowacc} <= ${cutoff2} && ${flowacc} > ${cutoff1}, b, c)) )", quiet = "True", qsx = qsx, kappa = kappa, slope = slope, aspect = aspect, R = R, K = K, C =C, res = r, flowacc = flowacc, Kt = Kt, rain = rain, stormtimet = stormtimet, loadexp = loadexp, cutoff1 = cutoff1, cutoff2 = cutoff2, cutoff3 = cutoff3)
- grass.mapcalc("${qsy}=eval(a=(${kappa} * sin(${slope}) * sin(${aspect})), b=((${R}*${K}*${C}*${flowacc}*${res}*sin(${slope})) * sin(${aspect})), c=( (${R}*${K}*${C}*exp((${flowacc}*${res}),1.6000000)*exp(sin(${slope}),1.3000000)) * sin(${aspect})), if(${flowacc} <= ${cutoff1}, a, if(${flowacc} <= ${cutoff2} && ${flowacc} > ${cutoff1}, b, c)) )", quiet = "True", qsy = qsy, kappa = kappa, slope = slope, aspect = aspect, R = R, K = K, C =C, res = r, flowacc = flowacc, Kt = Kt, rain = rain, stormtimet = stormtimet, loadexp = loadexp, cutoff1 = cutoff1, cutoff2 = cutoff2, cutoff3 = cutoff3)
- else:
- #This is the normal version (with 2D streams)
- qsx = "%sQsx_%04d" % (p,o)
- qsy = "%sQsy_%04d" % (p,o)
- if ( flags['c'] is True ): #do the shear stress version. Note that I'm converting the stream-power output (kg/m2) to same units as USPED (T/ha) by multiplying by ten. This ensures they are even going into the divergence calculation
- grass.mapcalc("${qsx}=eval(a=(${kappa} * sin(${slope}) * cos(${aspect})), b=((${R}*${K}*${C}*${flowacc}*${res}*sin(${slope})) * cos(${aspect})), c=( (${R}*${K}*${C}*exp((${flowacc}*${res}),1.6000000)*exp(sin(${slope}),1.3000000)) * cos(${aspect})), d=10 * (${Kt} * exp(9806.65*(((${rain}/1000)*${flowacc})/(0.595*${stormtimet}))*tan(${slope}), ${loadexp}) ) * cos(${aspect}), if(${flowacc} >= ${cutoff3}, a, if(${flowacc} >= ${cutoff2} && ${flowacc} < ${cutoff3}, b, if(${flowacc} >= ${cutoff1} && ${flowacc} < ${cutoff2}, c, d))) )", quiet = "True", qsx = qsx, kappa = kappa, slope = slope, aspect = aspect, R = R, K = K, C =C, res = r, flowacc = flowacc, Kt = Kt, rain = rain, stormtimet = stormtimet, loadexp = loadexp, cutoff1 = cutoff1, cutoff2 = cutoff2, cutoff3 = cutoff3)
- grass.mapcalc("${qsy}=eval(a=(${kappa} * sin(${slope}) * sin(${aspect})), b=((${R}*${K}*${C}*${flowacc}*${res}*sin(${slope})) * sin(${aspect})), c=( (${R}*${K}*${C}*exp((${flowacc}*${res}),1.6000000)*exp(sin(${slope}),1.3000000)) * sin(${aspect})), d=10 * (${Kt} * exp(9806.65*(((${rain}/1000)*${flowacc})/(0.595*${stormtimet}))*tan(${slope}), ${loadexp}) ) * sin(${aspect}), if(${flowacc} >= ${cutoff3}, a, if(${flowacc} >= ${cutoff2} && ${flowacc} < ${cutoff3}, b, if(${flowacc} >= ${cutoff1} && ${flowacc} < ${cutoff2}, c, d))) )", quiet = "True", qsy = qsy, kappa = kappa, slope = slope, aspect = aspect, R = R, K = K, C =C, res = r, flowacc = flowacc, Kt = Kt, rain = rain, stormtimet = stormtimet, loadexp = loadexp, cutoff1 = cutoff1, cutoff2 = cutoff2, cutoff3 = cutoff3)
- else: #do the stream powered version. Note that I'm converting the stream-power output (kg/m2) to same units as USPED (T/ha) by multiplying by ten. This ensures they are even going into the divergence calculation #Qs = Kt * n^-1 * 9810 * depth^1.6 * tan(slope)^1.5
- grass.mapcalc("${qsx}=eval(a=(${kappa} * sin(${slope}) * cos(${aspect})), b=((${R}*${K}*${C}*${flowacc}*${res}*sin(${slope})) * cos(${aspect})), c=( (${R}*${K}*${C}*exp((${flowacc}*${res}),1.6000000)*exp(sin(${slope}),1.3000000)) * cos(${aspect})), d=10 *${Kt} * exp(${manningn}, -1) * 9810 * exp((((${rain}/1000)*${flowacc})/(0.595*${stormtimet})), 1.6) * exp(tan(${slope}), 1.5) * cos(${aspect}), if(${flowacc} <= ${cutoff1}, a, if(${flowacc} <= ${cutoff2} && ${flowacc} > ${cutoff1}, b, if(${flowacc} <= ${cutoff3} && ${flowacc} > ${cutoff2}, c, d))) )", quiet = "True", qsx = qsx, kappa = kappa, slope = slope, aspect = aspect, R = R, K = K, C =C, res = r, flowacc = flowacc, Kt = Kt, rain = rain, stormtimet = stormtimet, loadexp = loadexp, cutoff1 = cutoff1, cutoff2 = cutoff2, cutoff3 = cutoff3, manningn = manningn)
- grass.mapcalc("${qsy}=eval(a=(${kappa} * sin(${slope}) * sin(${aspect})), b=((${R}*${K}*${C}*${flowacc}*${res}*sin(${slope})) * sin(${aspect})), c=( (${R}*${K}*${C}*exp((${flowacc}*${res}),1.6000000)*exp(sin(${slope}),1.3000000)) * sin(${aspect})), d=10 * ${Kt} * exp(${manningn}, -1) * 9810 * exp((((${rain}/1000)*${flowacc})/(0.595*${stormtimet})), 1.6) * exp(tan(${slope}), 1.5) * sin(${aspect}), if(${flowacc} <= ${cutoff1}, a, if(${flowacc} <= ${cutoff2} && ${flowacc} > ${cutoff1}, b, if(${flowacc} <= ${cutoff3} && ${flowacc} > ${cutoff2}, c, d))) )", quiet = "True", qsy = qsy, kappa = kappa, slope = slope, aspect = aspect, R = R, K = K, C =C, res = r, flowacc = flowacc, Kt = Kt, rain = rain, stormtimet = stormtimet, loadexp = loadexp, cutoff1 = cutoff1, cutoff2 = cutoff2, cutoff3 = cutoff3, manningn = manningn)
- #make a map of the total TC for debugging purposes
- #TC = "%sTC_%04d" % (p,o)
- #grass.mapcalc("${TC}=eval(a=${kappa} * sin(${slope}), b=${R}*${K}*${C}*${flowacc}*${res}*sin(${slope}), c=${R}* ${K}* ${C}* exp( (${flowacc}*${res}),1.6000000) * exp(sin(${slope}),1.3000000), d=10 * ${Kt} * exp(${manningn}, -1) * 9810 * exp( ( ( (${rain}/1000)*${flowacc}) / (0.595*${stormtimet}) ), 1.6) * exp(tan(${slope}), 1.5), if(${flowacc} >= ${cutoff3}, a, if(${flowacc} >= ${cutoff2} && ${flowacc} < ${cutoff3}, b, if(${flowacc} >= ${cutoff1} && ${flowacc} < ${cutoff2}, c, d) ) ) )", quiet = "True", TC = TC, kappa = kappa, slope = slope, aspect = aspect, R = R, K = K, C =C, res = r, flowacc = flowacc, Kt = Kt, rain = rain, stormtimet = stormtimet, loadexp = loadexp, cutoff1 = cutoff1, cutoff2 = cutoff2, cutoff3 = cutoff3, manningn = manningn)
- #/CHANGES
-
- grass.message('\n*************************\n Iteration %s -- ' % o + 'step 4: calculating divergence/difference of sediment transport for each process and the actual amount of erosion or deposition in vertical meters/cell/year\n*************************\n\n')
- #Here is where we figure out the change in transport capacity, and thus the actual amount of erosion an deposition that would occur. There are two ways of doing this. On planar and convex surfaces (i.e., ridgetops, flats, hillslopes), it is better to take the 2D divergence of sediment flux (we use r.slope.aspect to calculate this), but on highly convex surfaces (i.e., in channels) it is better to take the 1D difference between one cell, and the cell that is immediately downstream from it. This all assumes that the system is always operating at Transport Capacity, or if it is not, then is still behaves as if it were (ie., that the actual differences in transported sediment between the cells would be proportional to the system operating at capacity). Thus, under this assumption, the divergence of capacity is equals to actual amount of sediment eroded/deposited.
- #This is the way we implemnt this: First calculate, we calculate the divergence/differnce for EACH of the different flow processes on the ENTIRE map (i.e., make one map per process, difference for streams, divergence for USPED and diffusion). Then, we cut out the pieces of each of these maps that correspond to the correct landforms from each specific process (based on the user-input cutoffs in flow accumulation), and patch them together into a single map (NOTE: see output unit conversions section below to see how we get all the units to line up during this process). This counters the "boundary effect" that happens when running the differential equations for divergence across the boundary of two different flow processes. Then we may still have to run a median smoother on the patched map to get rid of any latent spikes.
- if ( flags["1"] is True ):
- #This is the version with 1D streams
- qsd1 = '%sDelta_Qs_1D_streams%04d' % (p, o)
- grass.mapcalc('${qsd1}=if(${flowdir} == 7, (${qs1}[-1,-1]-${qs1}), if (${flowdir} == 6, (${qs1}[-1,0]-${qs1}), if (${flowdir} == 5, (${qs1}[-1,1]-${qs1}), if (${flowdir} == 4, (${qs1}[0,1]-${qs1}), if (${flowdir} == 3, (${qs1}[1,1]-${qs1}), if (${flowdir} == 2, (${qs1}[1,0]-${qs1}), if (${flowdir} == 1, (${qs1}[1,-1]-${qs1}), if (${flowdir} == 8, (${qs1}[0,-1]-${qs1}), ${qs1}))))))))', quiet = "True", qsd1 = qsd1, flowdir = flowdir, qs1 = qs1)
- qsxdx = '%sDelta_Qsx_%04d' % (p, o)
- qsydy = '%sDelta_Qsy_%04d' % (p, o)
- grass.run_command('r.slope.aspect', quiet = "True", elevation = qsx, dx = qsxdx)
- grass.run_command('r.slope.aspect', quiet = "True", elevation = qsy, dy = qsydy)
- else:
- #This is the normal version (with 2D streams)
- qsxdx = '%sDelta_Qsx_%04d' % (p, o)
- qsydy = '%sDelta_Qsy_%04d' % (p, o)
- grass.run_command('r.slope.aspect', quiet = "True", elevation = qsx, dx = qsxdx)
- grass.run_command('r.slope.aspect', quiet = "True", elevation = qsy, dy = qsydy)
-
- #This is the smoothing routine. First we calculate the rate of Erosion and Deposition by converting the Delta QS of the different processes to vertical meters by dividing by the soil denisity (with apropriate constants to get into the correct units, see UNIT CONVERSION note below), and for streams, also expand from the storm to the year level. All units of this initial (temporary) ED_rate map will be in m/cell/year.
- #CHANGES
- #OUTPUT UNIT CONVERSIONS: In the case of the diffusion equation, the output units are in vertical meters of sediment per cell per year, so these will be left alone. Everything else should be in units of T/cell per storm. So we just need to convert to kg/cell, divide by the soil density and multiply the number of storms
- if ( flags["1"] is True ):
- #This is the version with 1D streams
- grass.mapcalc('${tempnetchange1}=if(${flowacc} >= ${cutoff3}, ((${qsd1}*0.1)/${sdensity})*${storms}, if(${flowacc} >= ${cutoff1} && ${flowacc} < ${cutoff3}, (((${qsxdx}+${qsydy})*0.1)/${sdensity})*${storms}, ${qsxdx}+${qsydy}))', quiet = "True", tempnetchange1 = tempnetchange1, qsd1 = qsd1, qsxdx = qsxdx, qsydy = qsydy, flowacc = flowacc, cutoff1 = cutoff1, cutoff3 = cutoff3, sdensity = sdensity, storms = storms, stormtimet = stormtimet)
- else:
- #This is the normal version (with 2D streams)
- grass.mapcalc('${tempnetchange1}=if(${flowacc} >= ${cutoff1}, (((${qsxdx} + ${qsydy})*0.1)/${sdensity})*${storms}, ${qsxdx}+${qsydy})', quiet = "True", tempnetchange1 = tempnetchange1, qsxdx = qsxdx, qsydy = qsydy, flowacc = flowacc, cutoff1 = cutoff1, cutoff3 = cutoff3, sdensity = sdensity, storms = storms, stormtimet = stormtimet)
- #grass.mapcalc('${tempnetchange1}=if(${flowacc} >= ${cutoff3}, (((${qsxdx} + ${qsydy})*0.1)/${sdensity})*${storms}, if(${flowacc} >= ${cutoff1} && ${flowacc} < ${cutoff3}, ((${qsxdx}+${qsydy})*0.1)/${sdensity}, ${qsxdx}+${qsydy}))', quiet = "True", tempnetchange1 = tempnetchange1, qsxdx = qsxdx, qsydy = qsydy, flowacc = flowacc, cutoff1 = cutoff1, cutoff3 = cutoff3, sdensity = sdensity, storms = storms, stormtimet = stormtimet)
- #/CHANGES
-
- #Make some temp maps of just erosion rate and just deposition rate so we can grab some stats from them for the soft-knee limiting filter
- grass.message('Running soft-knee smoothing filter...')
- grass.mapcalc('${tmperosion}=if(${tempnetchange1} < -0, ${tempnetchange1}, null())', quiet = "True", tmperosion = tmperosion, tempnetchange1 = tempnetchange1)
- grass.mapcalc('${tmpdep}=if(${tempnetchange1} > 0, ${tempnetchange1}, null())', quiet = "True", tmpdep = tmpdep, tempnetchange1 = tempnetchange1)
- #Grab the stats from these temp files and save them to dictionaries
- erosstats = grass.parse_command('r.univar', flags = 'ge', percentile = '1', map = tmperosion)
- depostats = grass.parse_command('r.univar', flags = 'ge', percentile = '99', map = tmpdep)
- maximum = depostats['max']
- minimum = erosstats['min']
- erosbreak = float(erosstats['first_quartile'])
- deposbreak = float(depostats['third_quartile'])
- scalemin = float(erosstats['percentile_1'])
- scalemax = float(depostats['percentile_99'])
- #Use the stats we gathered to do some smoothing with a hi-cut and lo-cut filter (with soft-knee limiting) of the unsmoothed ED_rate map. Values from the 1st quartile of erosion to the minimum (i.e., the very large negative numbers) will be rescaled linearly from the 1st quartile to the 1st percentile value, and values from the 3rd quartile of deposition to the maximum (i.e., the very large positiive numbers) will be rescaled linearly from the 3rd quartile to the 99th percentile value. This brings any values that were really unreasonnable as originally calculated (spikes) into the range of what the maximum values should be on a normally distrubuted dataset, but does so with out a "brick wall" style of limiting, which would make all values above some cutoff equal to a theoretical maximum. By setting both maximum cutoff point AND a "soft" scaling point, this "soft-knee" style of limiting sill retains some of the original scaling at the high ends, which allows for the smooth
ed value of very high cells to still be relatively higher than values in other cells that were also above the scaling cutoff, but were not originally as high as those very high cells.
- grass.mapcalc('${tempnetchange2}=graph(${tempnetchange1}, ${minimum},${scalemin}, ${erosbreak},${erosbreak}, ${deposbreak},${deposbreak}, ${maximum},${scalemax})', quiet = "True", tempnetchange2 = tempnetchange2, tempnetchange1 =tempnetchange1, minimum = minimum, scalemin = scalemin, erosbreak = erosbreak, deposbreak = deposbreak, maximum = maximum, scalemax = scalemax)
- #Check if additional smoothing is requested.
- if smoothing == "no":
- grass.message('No additional modal smoothing was requested...')
- grass.run_command('g.rename', quiet = "True", rast = tempnetchange2 + ',' + netchange)
- elif smoothing == "low":
- grass.message('Enacting additional "low" smoothing: one pass of a 3x3 modal smoothing window.')
- grass.run_command('r.neighbors', quiet = "True", input = tempnetchange2, output = netchange, method = 'mode', size = '3')
- elif smoothing == "high":
- grass.message('Enacting additional "high" smoothing: one pass of a 5x5 modal smoothing window.')
- grass.run_command('r.neighbors', quiet = "True", input = tempnetchange2, output = netchange, method = 'mode', size = '5')
- else:
- grass.message('There was a problem reading the median-smoothing variable, so maps will not be median-smoothed.')
- grass.run_command('g.rename', quiet = "True", rast = tempnetchange2 + ',' + netchange)
- #Set the netchange map colors to the rules we've provided above
- grass.run_command('r.colors', quiet = "True", map = netchange, rules = nccolors.name)
- #Grab the stats from these new smoothed netchange maps and save them to dictionaries (Note that the temporary erosion and deposition maps made in this step are overwriting the two temporary maps made for gathering the stats for the soft-knee limiting filter)
- grass.mapcalc('${tmperosion}=if(${netchange} < -0, ${netchange}, null())', quiet = "True", overwrite = "True", tmperosion = tmperosion, netchange = netchange)
- grass.mapcalc('${tmpdep}=if(${netchange} > 0, ${netchange}, null())', quiet = "True", overwrite = "True", tmpdep = tmpdep, netchange = netchange)
- erosstats1 = grass.parse_command('r.univar', flags = 'ge', map = tmperosion)
- depostats1 = grass.parse_command('r.univar', flags = 'ge', map = tmpdep)
-
- grass.message('\n*************************\n Iteration %s -- ' % o + 'step 5: calculating terrain evolution and new soil depths\n *************************\n\n')
- #Set up a temp dem, and then do initial addition of ED change to old DEM. This mapcalc statement first checks the amount of erodable soil in a given cell against the amount of erosion calculated, and keeps the cell from eroding past this amount (if there is soil, then if the amount of erosion is more than the amount of soil, just remove all the soil and stop, else remove the amount of caclulated erosion. It also runs an error catch that checks to make sure that soil depth is not negative (could happen, I suppose), and if it is, corrects it). Finally, do patch-job to catch the shrinking edge problem (the edge cells have no upstream cell, so get turned null in the calculations in step 4)
- grass.mapcalc('${new_dem}=eval(x=if(${old_soil} > 0.0 && (-1*${netchange}) <= ${old_soil}, ${netchange}, if((-1*${netchange}) > ${old_soil}, (-1*${old_soil}), 0)), y=(${old_dem} + x), if(isnull(y), ${old_dem}, y))', quiet = "True", new_dem = new_dem, old_soil = old_soil, old_dem = old_dem, netchange = netchange)
- #Set colors for elevation map to match other dems
- grass.run_command('r.colors', quiet = "True", map = new_dem, rast = options["elev"])
- grass.mapcalc('${new_soil}=if ((${new_dem} - ${initbdrk}) < 0, 0, (${new_dem} - ${initbdrk}))', quiet = "True", new_soil = new_soil, new_dem = new_dem, initbdrk = initbdrk)
- grass.run_command('r.colors', quiet = "True", map = new_soil, rules = sdcolors.name)
- grass.message('\n*************************\n Iteration %s -- ' % o + 'step 6: writing stats to output file\n *************************\n\n')
- #Finish gathering stats (just need the soil depth stats now)
- soilstats = grass.parse_command('r.univar', flags = 'ge', map = new_soil, percentile = '99')
- #Write stats to a new line in the stats file
- #HEADER of the file should be: ',,Mean Values,,,,Standard Deviations,,,,Totals,,,Additional Stats\nIteration,,Mean Erosion,Mean Deposition,Mean Soil Depth,,Standard Deviation Erosion,Standard Deviation Deposition,Standard Deviation Soil Depth,,Total Sediment Eroded,Total Sediment Deposited,,Minimum Erosion,First Quartile Erosion,Median Erosion,Third Quartile Erosion,Maximum Erosion,Original Un-smoothed Maximum Erosion,,Minimum Deposition,First Quartile Deposition,Median Deposition,Third Quartile Deposition,Maximum Deposition,Original Un-smoothed Maximum Deposition,,Minimum Soil Depth,First Quartile Soil Depth,Median Soil Depth,Third Quartile Soil Depth,Maximum Soil Depth'
- grass.message('Outputing stats to textfile: ' + q)
- f.write('\n%s' % o + ',,' + erosstats1['mean'] + ',' + depostats1['mean'] + ',' + soilstats['mean'] + ',,' + erosstats1['stddev'] + ',' + depostats1['stddev'] + ',' + soilstats['stddev'] + ',,' + erosstats1['sum'] + ',' + depostats1['sum'] + ',,' + erosstats1['max'] + ',' + erosstats1['third_quartile'] + ',' + erosstats1['median'] + ',' + erosstats1['first_quartile'] + ',' + erosstats1['min'] + ',' + minimum + ',,' + depostats1['min'] + ',' + depostats1['first_quartile'] + ',' + depostats1['median'] + ',' + depostats1['third_quartile'] + ',' + depostats1['max'] + ',' + maximum + ',,' + soilstats['min'] + ',' + soilstats['first_quartile'] + ',' + soilstats['median'] + ',' + soilstats['third_quartile'] + ',' + soilstats['max'])
-
- #Clean up temporary maps
- if flags["k"] is True:
- grass.message('\nTemporary maps will NOT be deleted!!!!\n')
- else:
- grass.message('\nCleaning up temporary maps...\n\n')
- #first remove all the easy temporary maps labeled with "pid"
- grass.run_command("g.remove", quiet = "True", flags = 'f', type = 'rast', pattern = '%s*' % pid)
- #now check all the flag options, and build a list of maps to delete
- mapstoremove = []
- if flags["s"] is True:
- grass.message('Keeping Slope map.')
- else:
- mapstoremove.append(slope)
- if flags["d"] is True :
- grass.message('Not keeping Soil Depth map.')
- mapstoremove.append(old_soil)
- #check if this is the last year and remove the "new-soil" map too
- if ( o == int(options["number"])):
- mapstoremove.append(new_soil)
- else:
- #check if this is the first year, and if so, remove the temporary "soildepths_init" map
- if ( o <= 1 ):
- mapstoremove.append("%s%s_init" % (prefx, options["outsoil"]))
- if flags["e"] is True :
- grass.message('Keeping Excess Transport Capacity (divergence) maps for all processes.')
- else:
- mapstoremove.extend([qsxdx, qsydy])
- if flags["1"] is True :
- mapstoremove.append(qsd1)
- if flags["t"] is True :
- grass.message('Keeping Transport Capacity maps for all processes.')
- else:
- mapstoremove.extend([qsx, qsy])
- if flags["1"] is True :
- mapstoremove.append(qs1)
- if flags["r"] is True :
- grass.message('Not keeping an Erosion and Deposition rate map.')
- mapstoremove.append(netchange)
- if len(mapstoremove) == 0:
- pass
- else:
- grass.run_command('g.remove', quiet = "True", flags = 'f', type = "rast", name = ','.join(mapstoremove))
- sdcolors.close()
- nccolors.close()
- grass.message('\n*************************\nDone with Iteration %s ' % o + '\n*************************\n')
- return(0)
-
-#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__":
- options, flags = grass.parser()
- # Set up some basic variables
- years = options["number"]
- prefx = options["prefx"]
- #these values could be read in from a climate file, so check that, and act accordingly. Either way, the result will be some lists with the same number of entries as there are iterations.
- rain2 = []
- try:
- rain1 = float(options["rain"])
- for year in range(int(years)):
- rain2.append(rain1)
- except:
- with open(options["rain"], 'rU') as f:
- for line in f:
- rain2.append(line.split(",")[0])
- #check for text header and remove if present
- try:
- float(rain2[0])
- except:
- del rain2[0]
- #throw a warning if there aren't enough values in the column
- if len(rain2) != int(years):
- grass.fatal("Number of rows of rainfall data in your climate file\n do not match the number of iterations you wish to run.\n Please ensure that these numbers match and try again")
- sys.exit(1)
- R2 = []
- try:
- R1 = float(options["r"])
- for year in range(int(years)):
- R2.append(R1)
- except:
- with open(options["r"], 'rU') as f:
- for line in f:
- R2.append(line.split(",")[1])
- #check for text header and remove if present
- try:
- float(R2[0])
- except:
- del R2[0]
- #throw a warning if there aren't enough values in the column
- if len(R2) != int(years):
- grass.fatal("Number of rows of R-Factor data in your climate file\n do not match the number of iterations you wish to run.\n Please ensure that these numbers match and try again")
- sys.exit(1)
- storms2 = []
- try:
- storms1 = float(options["storms"])
- for year in range(int(years)):
- storms2.append(storms1)
- except:
- with open(options["storms"], 'rU') as f:
- for line in f:
- storms2.append(line.split(",")[2])
- #check for text header and remove if present
- try:
- float(storms2[0])
- except:
- del storms2[0]
- #throw a warning if there aren't enough values in the column
- if len(storms2) != int(years):
- grass.fatal("Number of rows of storm frequency data in your climate file\n do not match the number of iterations you wish to run.\n Please ensure that these numbers match and try again")
- sys.exit(1)
- stormlength2 = []
- try:
- stormlength1 = float(options["stormlength"])
- for year in range(int(years)):
- stormlength2.append(stormlength1)
- except:
- with open(options["stormlength"], 'rU') as f:
- for line in f:
- stormlength2.append(line.split(",")[3])
- #check for text header and remove if present
- try:
- float(stormlength2[0])
- except:
- del stormlength2[0]
- #throw a warning if there aren't enough values in the column
- if len(stormlength2) != int(years):
- grass.fatal("Number of rows of storm length data in your climate file\n do not match the number of iterations you wish to run.\n Please ensure that these numbers match and try again")
- sys.exit(1)
- #Now gather these four lists into one master list, to make it easier to pass on to main()
- masterlist = [rain2,R2,storms2,stormlength2]
- #Make the statsout file with correct column headers
- if options["statsout"] == "":
- env = grass.gisenv()
- mapset = env['MAPSET']
- statsout = '%s_%slsevol_stats.csv' % (mapset, prefx)
- else:
- statsout = options["statsout"]
- if os.path.isfile(statsout):
- f = file(statsout, 'a')
- else:
- f = file(statsout, 'wt')
- f.write('These statistics are in units of vertical meters (depth) per cell\n,,Mean Values,,,,Standard Deviations,,,,Totals,,,Additional Stats\nIteration,,Mean Erosion,Mean Deposition,Mean Soil Depth,,Standard Deviation Erosion,Standard Deviation Deposition,Standard Deviation Soil Depth,,Total Sediment Eroded,Total Sediment Deposited,,Minimum Erosion,First Quartile Erosion,Median Erosion,Third Quartile Erosion,Maximum Erosion,Original Un-smoothed Maximum Erosion,,Minimum Deposition,First Quartile Deposition,Median Deposition,Third Quartile Deposition,Maximum Deposition,Original Un-smoothed Maximum Deposition,,Minimum Soil Depth,First Quartile Soil Depth,Median Soil Depth,Third Quartile Soil Depth,Maximum Soil Depth')
- if flags["p"] is True :
- grass.message('Making sample points map for determining cutoffs.')
- else:
- grass.message('\n##################################################\n##################################################\n\n STARTING SIMULATION\n\nBeginning iteration sequence. This may take some time.\nProcess is not finished until you see the message: \'Done with everything\'\n _____________________________________________________________\n_____________________________________________________________\n')
- grass.message("Total number of iterations to be run is %s" % years)
- #Get the region settings
- region1 = grass.region()
- # This is the loop!
- for x in range(int(years)):
- grass.message("Iteration = %s" % (x + 1))
- main(x, (x + 1), prefx, statsout, region1['nsres'], masterlist);
- #Since we are now done with the loop, close the stats file.
- f.close()
- grass.message('\nIterations complete!\n\nDone with everything')
- sys.exit(0)
-
-
-
-
Deleted: grass-addons/grass7/raster/r.landscape.evol/r.landscape.evol description.odf.odt
===================================================================
(Binary files differ)
Deleted: grass-addons/grass7/raster/r.landscape.evol/r.landscape.evol description.odt
===================================================================
(Binary files differ)
Copied: grass-addons/grass7/raster/r.landscape.evol/r.landscape.evol.html (from rev 64600, grass-addons/grass7/raster/r.landscape.evol/description.html)
===================================================================
--- grass-addons/grass7/raster/r.landscape.evol/r.landscape.evol.html (rev 0)
+++ grass-addons/grass7/raster/r.landscape.evol/r.landscape.evol.html 2015-02-13 09:49:02 UTC (rev 64601)
@@ -0,0 +1,366 @@
+<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 either 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><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>and <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. Note that some
+experimentation is required in order to find the best possible values
+for these cutoffs, and the -p flag will provide some output data that
+may be useful for this. <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>
+is the total annual precipitation measured in meters (or the average
+annual rainfall in meters per year). <b>raindays</b>
+is the total number of days on which it rained in one year (or an
+average value of days per year). <b>infilt</b>
+is the proportion of rainfall that infiltrates into the soil and thus
+does not contribute to runoff (values are between 0 and 1). <b>Kt</b>
+is the stream transport efficiency variable that describes the
+cohesiveness of the stream channel beds (0.001 for normal
+gravel/sandy/silt channel bed to 0.000001 for a bedrock channel bed).
+<b>loadexp</b> 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). <b>alpha</b> is
+the critical slope threshold above which the model will simulate the
+cumulative effects of mass wasting (landsliding). These
+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>By default, <em>r.watershed</em> is used to calculate flow
+accumulation modeling using the MFD algorithm included in GRASS 6.4
+and higher. This can be made backwards compatible by checking the -f
+flag, which will use <i>r.terraflow </i>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.
+<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).
+
+<h3>Calculating Erosion and Deposition</h3>
+<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>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)
+<center>
+<img src="r_landscape_evol_equation1.gif"><br>
+</center>
+<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>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>),
+and R-factor (<i>R</i>) components
+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>For areas of the DEM that have flow accumulation values greater
+than <b>cutoff3 </b>(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 (<FONT FACE="Times New Roman, serif">τ</FONT>),
+here estimated for a cellular landscape simply as:
+<center>
+<p><img src="r_landscape_evol_equation2.gif"><br>
+</center>
+<p> Where: <i>9806.65</i>
+is a constant related to the gravitational acceleration of water, <i>B</i>
+is the slope of the cell in degrees, and <i>D</i>
+is the instantaneous depth of flowing water in the cell. <i>D
+</i>is
+here assumed to be roughly equivalent to the depth of flow during the
+average minute of rainfall, calculated by:
+<center>
+<img src="r_landscape_evol_equation3.gif"><br>
+</center>
+<p>Where: <i>R</i><sub><i>m</i></sub>
+is the total annual precipitation in meters, <i>i</i>
+is the proportion of rainfall that infiltrates rather than runs
+off, <i>A</i>
+is the uplsope accumulated area per unit contour width at the cell,
+<i>R</i><sub><i>d</i></sub>
+is the number of days on which it rained in a one year period, and
+<i>1440</i>
+is a constant relating to the number of minutes in a day.
+<p>Then the transport capacity is calculated by:
+<center>
+<img src="r_landscape_evol_equation4.gif"><br>
+</center>
+<p>Where: <i>K</i><sub><i>t</i></sub>
+is the transport efficiency factor related to the character of the
+stream bed (0.001 for normal sediment to 0.000001 for bedrock), and <i>n</i>
+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).
+<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:
+<center>
+<img src="r_landscape_evol_equation5.gif"><br>
+</center>
+<p>where ED is net erosion or
+deposition rate for sediment and <em><FONT FACE="Times New Roman, serif">α</FONT></em>
+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 <em><U>change</U></em>
+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.
+<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 neighborhood 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>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 bedrock
+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>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:
+<center>
+<img src="r_landscape_evol_equation6.gif"><br>
+</center>
+<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 raster cell resolution x <b>sdensity</b>)" to
+create a map of soil erosion/deposition rates across the landscape.
+
+<h3>Determining Cutoff Points</h3>
+<p>
+To get started with <em>r.landscape.evol</em>, you need to determine the appropriate values for <b>cutoff1</b>, <b>cutoff2</b>, and <b>cutoff3</b>, which are transition points between different types of erosive processes. These are in units of flow accumulation scaled to actual surface flow as determined in r.watershed from the values of rainfall and flow hindrance from vegetation. To do this, you should parameterize the module as best as possible, EXCEPT for the three "cutoffs". Then, run the module with the <b>-p</b> flag, which will make a random points vector file with the values of scaled flow accumulation (scaled to actual rainfall and vegetation), profile curvature, and tangential curvature in the associated table. Plotting the log of the scaled flow accumulation against each of these two curvatures will help you to determine reasonable values for the cutoffs, as each transition should show a unique relationship between curvature and flow accumulations. See
the figures below for examples:
+
+<center>
+<img src="r_landscape_evol_Flow_acc_vs_curvature.png" width="1000" height="500" alt="Log Scaled Flow Accumulation versus Topographic Curvatures"><br>
+
+Log Scaled Flow Accumulation versus Topographic Curvatures.<br><br>
+
+<img src="r_landscape_evol_Map1.png" width="500" height="284" alt="Map showing the spatial patterns of the cutoffs determined from the previous figure"><br>
+
+Map showing the spatial patterns of the cutoffs determined from the previous figure.<br><br>
+</center>
+
+<p>
+<h3>Note About Climate Parameters</h3>
+<p>
+r.landscape.evol accepts an external "climate file", which should be a comma separated plain text file with four columns in the order of, "<b>rain</b>,<b>R</b>,<b>storms</b>,<b>stormlength</b>" (without headers). Each of these columns must exist, although there need not be values in every column (i.e., you can enter a single value for any of these parameters in the command line, and combine that with populated columns for the other values). Note that the climate file must have the same number of rows as there are iterations of the simulation (<b>years</b>).
+
+<h2>SEE ALSO</h2>
+<ul>
+ <li><p>The <a href="http://medland.asu.edu/">MEDLAND</a>
+ project at Arizona State University
+
+ <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>
+
+ <li><p>Mitasova, H., C. M. Barton, I. I. Ullah, J. Hofierka, and R. S. Harmon 2013 GIS-based soil erosion modeling. In Remote Sensing and GIScience in Geomorphology, edited by J. Shroder and M. P. Bishop. 3:228-258. San Diego: Academic Press.
+
+
+</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>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>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>Howard, A. D. 1980. Thresholds in river regimes. Thresholds
+in geomorphology, 227-258.
+
+<p>Mitas, L. and H. Mitasova 1998 Distributed soil erosion simulation
+for effective erosion prevention. Water Resources Research
+34(3):505-516.
+
+<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>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>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>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>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>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>Rawls, W. J. 1983 Estimating soil bulk denisty from particle size
+analysis and organic matter content. Soil Science 135(2):123.
+
+<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-51. vol. 703. US Department of
+Agriculture, Washington, DC.
+
+<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>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>Tucker, G. E. and G. R Hancock 2010 Modelling landscape
+evolution. Earth Surface Processes
+and Landforms 35(1): 28-50.
+
+<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>Wischmeier, W. H. 1976 Use and Misuse of the Universal Soil Loss
+Equation. Journal of Soil and Water Conservation 31:5-9.
+
+<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>Wischmeier, W. H. and D. D. Smith 1978 Predicting Rainfall-Erosion
+Losses - A Guide to Conservation Planning. USDA Agriculture Handbook
+282.
+
+
+<h2>AUTHORS</h2>
+Isaac I. Ullah, C. Michael Barton, and Helena Mitasova
+
+<p><i>Last changed: $Date: 2014-010-02 05:03:58 -0400 (Thur, 02 Oct 2014) $</i>
Copied: grass-addons/grass7/raster/r.landscape.evol/r.landscape.evol.py (from rev 64600, grass-addons/grass7/raster/r.landscape.evol/r.landscape.evol)
===================================================================
--- grass-addons/grass7/raster/r.landscape.evol/r.landscape.evol.py (rev 0)
+++ grass-addons/grass7/raster/r.landscape.evol/r.landscape.evol.py 2015-02-13 09:49:02 UTC (rev 64601)
@@ -0,0 +1,687 @@
+#!/usr/bin/python
+
+############################################################################
+#
+# MODULE: r.landscape.evol.py
+# AUTHOR(S): Isaac Ullah and Michael Barton
+# COPYRIGHT: (C) 2012 GRASS Development Team/Isaac Ullah
+#
+# description: Simulates the cumulative effect of erosion and deposition on a landscape over time. This module uses appropriate flow on different landforms by default; however, singular flow regimes can be chosen by manipulating the cutoff points. This module requires GRASS 6.4 or greater. THIS SCRIPT WILL PRODUCE MANY TEMPORARY MAPS AND REQUIRES A LOT OF FREE FILE SPACE!
+
+# This program is free software; you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# 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: Simulates the cumulative effect of erosion and deposition on a landscape over time. This module uses appropriate flow on different landforms by default; however, singular flow regimes can be chosen by manipulating the cutoff points. This module requires GRASS 6.4 or greater. THIS SCRIPT WILL PRODUCE MANY TEMPORARY MAPS AND REQUIRES A LOT OF FREE FILE SPACE!
+#%End
+#%option
+#% key: elev
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: Input elevation map (DEM of surface)
+#% required : yes
+#%end
+#%option
+#% key: initbdrk
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: Bedrock elevations map (DEM of bedrock)
+#% answer:
+#% required : yes
+#%end
+#%option
+#% key: k
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: Soil erodability index (K factor) map or constant
+#% answer: 0.42
+#% required : no
+#% guisection: Landscape Evolution
+#%end
+#%option
+#% key: sdensity
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: Soil density map or constant [T/m3] for conversion from mass to volume
+#% answer: 1.2184
+#% required : no
+#% guisection: Landscape Evolution
+#%end
+#%option
+#% key: kt
+#% type: double
+#% description: Stream transport efficiency variable (0.001 for a soft substrate, 0.0001 for a normal substrate, 0.00001 for a hard substrate, 0.000001 for a very hard substrate)
+#% answer: 0.0001
+#% required : no
+#% options : 0.001,0.0001,0.00001,0.000001
+#% guisection: Landscape Evolution
+#%end
+#%option
+#% key: loadexp
+#% type: double
+#% description: Stream transport type variable (1.5 for mainly bedload transport, 2.5 for mainly suspended load transport)
+#% answer: 1.5
+#% options: 1.5,2.5
+#% required : no
+#% guisection: Landscape Evolution
+#%end
+#%option
+#% key: kappa
+#% type: double
+#% description: Hillslope diffusion (Kappa) rate map or constant [m/kyr]
+#% answer: 1
+#% required : no
+#% guisection: Landscape Evolution
+#%end
+#%option
+#% key: c
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: Landcover index (C factor) map or constant
+#% answer: 0.005
+#% required : no
+#% guisection: Landscape Evolution
+#%end
+#%option
+#% key: rain
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: Precip totals for the average storm [mm] (or path to climate file of comma separated values of "rain,R,storms,stormlength", with a new line for each year of the simulation)
+#% answer: 20.61
+#% guisection: Climate
+#%end
+#%option
+#% key: r
+#% type: string
+#% description: Rainfall (R factor) constant (AVERAGE FOR WHOLE MAP AREA) (or path to climate file of comma separated values of "rain,R,storms,stormlength", with a new line for each year of the simulation)
+#% answer: 4.54
+#% guisection: Climate
+#%end
+#%option
+#% key: storms
+#% type: string
+#% description: Number of storms per year (integer) (or path to climate file of comma separated values of "rain,R,storms,stormlength", with a new line for each year of the simulation)
+#% answer: 25
+#% guisection: Climate
+#%end
+#%option
+#% key: stormlength
+#% type: string
+#% description: Average length of the storm [h] (or path to climate file of comma separated values of "rain,R,storms,stormlength", with a new line for each year of the simulation)
+#% answer: 24.0
+#% guisection: Climate
+#%end
+#%option
+#% key: speed
+#% type: double
+#% description: Average velocity of flowing water in the drainage [m/s]
+#% answer: 1.4
+#% required : no
+#% guisection: Hydrology
+#%end
+#%option
+#% key: manningn
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: Map or constant of the value of Manning's "N" value for channelized flow.
+#% answer: 0.05
+#% required : no
+#% guisection: Hydrology
+#%end
+#%option
+#% key: flowcontrib
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: Map or constant indicating how much each cell contributes to downstream flow (as a "percentage" from 0-100). If no map or value entered, routine will assume 100% downstream contribution
+#% required : no
+#% guisection: Hydrology
+#%end
+#%option
+#% key: convergence
+#% type: integer
+#% description: Value for the flow convergence variable in r.watershed. Small values make water spread out, high values make it converge in narrower channels.
+#% answer: 5
+#% options: 1,2,3,4,5,6,7,8,9,10
+#% required : no
+#% guisection: Hydrology
+#%end
+#%option
+#% key: cutoff1
+#% type: double
+#% description: Flow accumulation breakpoint value for shift from diffusion to overland flow
+#% answer: 0
+#% required : no
+#% guisection: Hydrology
+#%end
+#%option
+#% key: cutoff2
+#% type: double
+#% description: Flow accumulation breakpoint value for shift from overland flow to rill/gully flow (if value is the same as cutoff1, no sheetwash procesess will be modeled)
+#% answer: 100
+#% required : no
+#% guisection: Hydrology
+#%end
+#%option
+#% key: cutoff3
+#% type: double
+#% description: Flow accumulation breakpoint value for shift from rill/gully flow to stream flow (if value is the same as cutoff2, no rill procesess will be modeled)
+#% answer: 100
+#% required : no
+#% guisection: Hydrology
+#%end
+#%option
+#% key: smoothing
+#% type: string
+#% description: Amount of additional smoothing (answer "no" unless you notice large spikes in the erdep rate map)
+#% answer: no
+#% options: no,low,high
+#% required : yes
+#%end
+#%option
+#% key: prefx
+#% type: string
+#% description: Prefix for all output maps
+#% answer: levol_
+#% required : yes
+#%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
+#%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
+#%end
+#%option
+#% key: number
+#% type: integer
+#% description: Number of iterations (cycles) to run
+#% answer: 1
+#% required : yes
+#%end
+
+
+# #%option
+# #% key: alpha ###This may be added back in if I can find a good equation for mass movement
+# #% type: integer
+# #% description: Critical slope threshold for mass movement of sediment (in degrees above horizontal)
+# #% answer: 40
+# #% required : yes
+# #%end
+
+#%flag
+#% key: p
+#% description: -p Output a vector points map with sampled values of flow accumulation and curvatures suitable for determining cutoff values. NOTE: Overrides all other output options, and exits after completion. The output vector points map will be named "PREFIX_#_randomly_sampled_points".
+#% guisection: Optional
+#%end
+#%flag
+#% key: 1
+#% description: -1 Calculate streams as 1D difference instead of 2D divergence
+#% guisection: Landscape Evolution
+#%end
+#%flag
+#% key: c
+#% description: -c Calculate streams with a shear stress equation, rather than a stream-power equation
+#% guisection: Landscape Evolution
+#%end
+#%flag
+#% key: k
+#% description: -k Keep ALL temporary maps (overides flags -drst). This will make A LOT of maps!
+#% guisection: Optional
+#%end
+#%flag
+#% key: d
+#% description: -d Don't output yearly soil depth maps
+#% guisection: Optional
+#%end
+#%flag
+#% key: r
+#% description: -r Don't output yearly maps of the erosion/deposition rates ("ED_rate" map, in vertical meters)
+#% guisection: Optional
+#%end
+#%flag
+#% key: s
+#% description: -s Keep all slope maps
+#% guisection: Optional
+#%end
+#%flag
+#% key: t
+#% description: -t Keep yearly maps of the Transport Capacity at each cell ("Qs" maps)
+#% guisection: Optional
+#%end
+#%flag
+#% key: e
+#% description: -e Keep yearly maps of the Excess Transport Capacity (divergence) at each cell ("DeltaQs" maps)
+#% guisection: Optional
+#%end
+#%Option
+#% key: statsout
+#% type: string
+#% description: Name for the statsout text file (optional, if none provided, a default name will be used)
+#% required: no
+#% guisection: Optional
+#%end
+
+import sys
+import os
+import math
+import tempfile
+grass_install_tree = os.getenv('GISBASE')
+sys.path.append(grass_install_tree + os.sep + 'etc' + os.sep + 'python')
+import grass.script as grass
+
+# 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, s = master list of lists of climate data
+def main(m, o, p, q, r, s):
+ #get the process id to tag any temporary maps we make for easy clean up in the loop
+ pid = os.getpid()
+ # Get variables from user input
+ smoothing = options["smoothing"]
+ years = options["number"]
+ initbdrk = options["initbdrk"]
+ outdem = options["outdem"]
+ outsoil = options["outsoil"]
+ K = options["k"]
+ sdensity = options["sdensity"]
+ C = options["c"]
+ kappa = options["kappa"]
+ cutoff1 = options["cutoff1"]
+ cutoff2 = options["cutoff2"]
+ cutoff3 = options["cutoff3"]
+ flowcontrib = options["flowcontrib"]
+ convergence = options["convergence"]
+ manningn = options["manningn"]
+ old_bdrk = options["initbdrk"]
+ # these variables come in as a list of lists, so let's get this year's numbers out of them.
+ rain = s[0][m]
+ storms = s[2][m]
+ #CHANGES
+ R = s[1][m] / storms
+ #R = s[1][m]
+ stormlengthsecs = float(s[3][m])*3600.00 # number of seconds in the storm
+ stormtimet = stormlengthsecs / (float(options["speed"]) * float(r)) # number of hydrologic instants in the storm
+# timet = stormlengthsecs/stormtimet # length of a single hydrologic instant in seconds, currently unused, but might be important in future versions
+ Kt = options["kt"]
+ loadexp = options["loadexp"]
+ # Make some variables for temporary map names, labeled different depending on if we keep them or not
+ if ( flags["k"] is True ):
+ aspect = '%saspect%04d' % (p, o)
+ flowacc = '%sflowacc%04d' % (p, o)
+ flacclargenums = '%sflowacc_largenums%04d' % (p, o)
+ flowdir = '%sflowdir%04d' % (p, o)
+ pc = '%spc%04d' % (p, o)
+ tc = '%stc%04d' % (p, o)
+# meancurv = '%smeancurv%04d' % (p, o) # This variable might be used if bedrock weathering is ever implemented
+# rate = '%srate%04d' % (p, o) # This variable might be used if bedrock weathering is ever implemented
+ rainexcess = "%s_rainfall_excess_map_%04d"% (p, o)
+ tempnetchange1 = '%sTEMPORARY_unsmoothed_ED_rate%04d' % (p, o)
+ tempnetchange2 = '%sTEMPORARY_smoothed_ED_rate%04d' % (p, o)
+ tmperosion = '%sTEMPORARY_erosion%04d' % (p, o)
+ tmpdep = '%sTEMPORARY_deposition%04d' % (p, o)
+ else:
+ aspect = '%saspect%04d' % (pid, o)
+ flowacc = '%sflowacc%04d' % (pid, o)
+ flacclargenums = '%sflowacc_largenums%04d' % (pid, o)
+ flowdir = '%sflowdir%04d' % (pid, o)
+ pc = '%spc%04d' % (pid, o)
+ tc = '%stc%04d' % (pid, o)
+# meancurv = '%smeancurv%04d' % (pid, o) # This variable might be used if bedrock weathering is ever implemented
+# rate = '%srate%04d' % (p, o) # This variable might be used if bedrock weathering is ever implemented
+ rainexcess = "%s_rainfall_excess_map_%04d"% (pid, o)
+ tempnetchange1 = '%sTEMPORARY_unsmoothed_ED_rate%04d' % (pid, o)
+ tempnetchange2 = '%sTEMPORARY_smoothed_ED_rate%04d' % (pid, o)
+ tmperosion = '%sTEMPORARY_erosion%04d' % (pid, o)
+ tmpdep = '%sTEMPORARY_deposition%04d' % (pid, 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' % options["elev"]
+ old_soil = "%s%s_init" % (prefx, options["outsoil"])
+ grass.mapcalc('${old_soil}=${old_dem}-${old_bdrk}', overwrite = "True", quiet = "True", old_soil = old_soil, old_dem = old_dem, old_bdrk = old_bdrk)
+ else :
+ old_dem = '%s%s%04d' % (p, options["outdem"], m)
+ old_soil = '%s%s%04d' % (p, options["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 = '%sED_rate' % p
+ new_dem ='%s%s' % (p, outdem)
+ new_soil = '%s%s' % (p, outsoil)
+ else:
+ slope = '%sslope%04d' % (p, o)
+ netchange = '%sED_rate%04d' % (p, o)
+ new_dem = '%s%s%04d' % (p, outdem, o)
+ new_soil = '%s%s%04d' % (p, outsoil, o)
+ #Check to see if we are going to only output diagnostics for determing cutoff values, and act accordingly
+ if ( flags["p"] is True ):
+ grass.message('GATHERING STATISTICS FOR DETERMINING CUTOFF VALUES\n-------------------------------------------------\n1) Calculating slope and curvatures')
+ grass.run_command('r.slope.aspect', quiet = "True", elevation = old_dem, slope = slope, pcurv = pc, tcurv = tc)
+ else:
+ grass.message('\n##################################################\n\n*************************\n Iteration %s -- ' % o + 'step 1: calculating slope\n*************************\n')
+ grass.run_command('r.slope.aspect', quiet = "True", elevation = old_dem, aspect = aspect, slope = slope)
+ if ( flags["p"] is True ):
+ grass.message('2) Calculating map of rainfall excess')
+ else:
+ grass.message('\n*************************\n Iteration %s -- ' % o + 'step 2: calculating accumulated flow depths\n*************************\n')
+ grass.message('Calculating runoff excess rates (uplsope accumulated cells scaled to \"flowcontrib\" map')
+
+ #make map of rainfall excess (proportion each cell contributes to downstrem flow) from flowcontrib. Note that if flowcontrib is a map, we are just making a copy of it. This map is a percentage, but has to be scale from 0-100, because r.watershed will only allow values greater than 1 as input in it's 'flow' variable. This creates a flow accumulation map with large numbers, but this map will be divided by 100 after it is made, which brings the values back down to what they should be.
+ if flowcontrib == "":
+ flowcontrib = 100
+ grass.mapcalc('${rainexcess}=int(${flowcontrib})', quiet = "True", rainexcess = rainexcess, flowcontrib = flowcontrib)
+ if ( os.getenv("GIS_FLAG_p") == "1" ):
+ grass.message('3) Calculating accumulated flow (in numbers of upslope cells, scaled by runoff contribution')
+ grass.run_command('r.watershed', quiet = "True", flags = 'a', elevation = old_dem, flow = rainexcess, accumulation = flacclargenums, drainage = flowdir, convergence = convergence)
+ grass.mapcalc('${flowacc}=${flacclargenums}/100', quiet = "True", flowacc = flowacc, flacclargenums = flacclargenums)
+ #again, do something different if we are only making an evaluation of cutoffs
+ if ( flags["p"] is True ):
+ grass.message('4) Determining number of sampling points using formula: "ln(#cells_in_input_map)*100"')
+ flaccstats = grass.parse_command('r.univar', flags = 'g', map = flowacc)
+ numpts = int(math.log(int(flaccstats['n']))*100)
+ grass.message('5) Creating random points and sampling values of flow accumulation, curvatures, and slope.')
+ vout = '%s%s_randomly_sampled_points' % (p, numpts)
+ grass.run_command('r.random', quiet = "True", input = flowacc, cover = pc, n = numpts, vector_output = vout)
+ grass.run_command('v.db.renamecol', quiet = "True", map = vout, column = 'value,Flow_acc')
+ grass.run_command('v.db.renamecol', quiet = "True", map = vout, column = 'covervalue,Princ_curv')
+ grass.run_command('v.db.addcol', quiet = "True", map = vout, columns = 'Tang_curv double precision, Slope double precision')
+ grass.run_command('v.what.rast', quiet = "True", vector = vout, raster = tc, column = "Tang_curv")
+ grass.run_command('v.what.rast', quiet = "True", vector = vout, raster = slope, column = "Slope")
+ if ( flags["k"] is True ):
+ grass.message('--Keeping the created maps (Flow Accumulation, Slope, Principle Curvature, Tangential Curvature)')
+ else:
+ grass.message('6) Cleaning up...')
+ grass.run_command('g.remove', quiet = "True", flags = 'f', type = "rast", name = slope + "," + pc + "," + tc + "," + flowacc)
+ grass.message('FINISHED. \nRandom sample points map "%s" created successfully.\n' % vout)
+ sys.exit(0)
+ grass.message('\n*************************\n Iteration %s -- ' % o + 'step 3: calculating sediment transport rates (units variable depending upon process) \n*************************\n')
+ # This step calculates the force of the flowing water at every cell on the landscape using the proper transport process law for the specific point in the flow regime. For upper hillslopes (below cutoff point 1) this done by multiplying the diffusion coeficient by the accumulated flow/cell res width. For midslopes (between cutoff 1 and 2) this is done by multiplying slope by accumulated flow with the m and n exponents set to 1. For channel catchment heads (between cutoff 2 and 3), this is done by multiplying slope by accumulated flow with the m and n exponents set to 1.6 and 1.3 respectively. For Channelized flow in streams (above cutoff 3), this is done by calculating the reach average shear stress (hydraulic radius [here estimated for a cellular landscape simply as the depth of flow] times slope times accumulated flow [cells] times gravitatiopnal acceleration of water [9806.65 newtons], all raised to the appropriate exponant for the type of transport (bedload or suspe
nded load), and then divided by the resolution. Depth of flow is calculated as a mean "instantaneous depth" during any given rain event, here estimated by the maximum depth of an idealized unit hydrograph with base equal to the duration of the storm, and area equal to the total accumulated excess rainfall during the storm. Then finally calculates the stream power or sediment carrying capacity (qs) of the water flowing at each part of the map by multiplying the reach average shear stress (channelized flow in streams) or the estimated flow force (overland flow) by the transport coeficient (estimated by R*K*C for hillslopes or kt for streams). This is a "Transport Limited" equation, however, we add some constraints on detachment by checking to see if the sediment supply has been exhausted: if the current soil depth is 0 or negative (checking for a negative value is kind of an error trap) then we make the transport coefficient small (0.000001) to simulate erosion on bedrock. Bec
ause diffusion and USPED require 2D divergence later on, we calculate these as vectors in the X and Y directions. Stream flow only needs 1D difference, so it's calulated in the direction of flow.
+ if ( flags["1"] is True ):
+ #This is the version with 1D streams
+ qs1 = '%sQs_1D_streams%04d' % (p, o)
+ #CHANGES
+ #choose shear stress or stream power
+ #these are the stream-power versions * note that I'm converting the stream-power output (kg/m2) to same units as USPED (T/ha) by multiplying by ten. This ensures they are even going into the divergence calculation #Qs = Kt * n^-1 * 9810 * depth^1.6 * tan(slope)^1.5
+ if ( flags['c'] is True ):
+ grass.mapcalc('${qs1}=(${Kt} * exp(9806.65*(((${rain}/1000)*${flowacc})/(0.595*${stormtimet}))*tan(${slope}), ${loadexp}) )', quiet = "True", qs1 = qs1, flowacc = flowacc, stormtimet = stormtimet, rain = rain, slope = slope, loadexp = loadexp, Kt = Kt, sdensity = sdensity)
+ else:
+ grass.mapcalc('${qs1}=10 * ${Kt} * exp(${manningn}, -1) * 9810 * exp( ( ( (${rain}/1000)*${flowacc}) / (0.595*${stormtimet}) ), 1.6) * exp(tan(${slope}), 1.5)', quiet = "True", qs1 = qs1, flowacc = flowacc, stormtimet = stormtimet, rain = rain, slope = slope, loadexp = loadexp, Kt = Kt, sdensity = sdensity, manningn = manningn)
+ qsx = "%sQsx_%04d" % (p,o)
+ qsy = "%sQsy_%04d" % (p,o)
+ grass.mapcalc("${qsx}=eval(a=(${kappa} * sin(${slope}) * cos(${aspect})), b=((${R}*${K}*${C}*${flowacc}*${res}*sin(${slope})) * cos(${aspect})), c=( (${R}*${K}*${C}*exp((${flowacc}*${res}),1.6000000)*exp(sin(${slope}),1.3000000)) * cos(${aspect})), if(${flowacc} <= ${cutoff1}, a, if(${flowacc} <= ${cutoff2} && ${flowacc} > ${cutoff1}, b, c)) )", quiet = "True", qsx = qsx, kappa = kappa, slope = slope, aspect = aspect, R = R, K = K, C =C, res = r, flowacc = flowacc, Kt = Kt, rain = rain, stormtimet = stormtimet, loadexp = loadexp, cutoff1 = cutoff1, cutoff2 = cutoff2, cutoff3 = cutoff3)
+ grass.mapcalc("${qsy}=eval(a=(${kappa} * sin(${slope}) * sin(${aspect})), b=((${R}*${K}*${C}*${flowacc}*${res}*sin(${slope})) * sin(${aspect})), c=( (${R}*${K}*${C}*exp((${flowacc}*${res}),1.6000000)*exp(sin(${slope}),1.3000000)) * sin(${aspect})), if(${flowacc} <= ${cutoff1}, a, if(${flowacc} <= ${cutoff2} && ${flowacc} > ${cutoff1}, b, c)) )", quiet = "True", qsy = qsy, kappa = kappa, slope = slope, aspect = aspect, R = R, K = K, C =C, res = r, flowacc = flowacc, Kt = Kt, rain = rain, stormtimet = stormtimet, loadexp = loadexp, cutoff1 = cutoff1, cutoff2 = cutoff2, cutoff3 = cutoff3)
+ else:
+ #This is the normal version (with 2D streams)
+ qsx = "%sQsx_%04d" % (p,o)
+ qsy = "%sQsy_%04d" % (p,o)
+ if ( flags['c'] is True ): #do the shear stress version. Note that I'm converting the stream-power output (kg/m2) to same units as USPED (T/ha) by multiplying by ten. This ensures they are even going into the divergence calculation
+ grass.mapcalc("${qsx}=eval(a=(${kappa} * sin(${slope}) * cos(${aspect})), b=((${R}*${K}*${C}*${flowacc}*${res}*sin(${slope})) * cos(${aspect})), c=( (${R}*${K}*${C}*exp((${flowacc}*${res}),1.6000000)*exp(sin(${slope}),1.3000000)) * cos(${aspect})), d=10 * (${Kt} * exp(9806.65*(((${rain}/1000)*${flowacc})/(0.595*${stormtimet}))*tan(${slope}), ${loadexp}) ) * cos(${aspect}), if(${flowacc} >= ${cutoff3}, a, if(${flowacc} >= ${cutoff2} && ${flowacc} < ${cutoff3}, b, if(${flowacc} >= ${cutoff1} && ${flowacc} < ${cutoff2}, c, d))) )", quiet = "True", qsx = qsx, kappa = kappa, slope = slope, aspect = aspect, R = R, K = K, C =C, res = r, flowacc = flowacc, Kt = Kt, rain = rain, stormtimet = stormtimet, loadexp = loadexp, cutoff1 = cutoff1, cutoff2 = cutoff2, cutoff3 = cutoff3)
+ grass.mapcalc("${qsy}=eval(a=(${kappa} * sin(${slope}) * sin(${aspect})), b=((${R}*${K}*${C}*${flowacc}*${res}*sin(${slope})) * sin(${aspect})), c=( (${R}*${K}*${C}*exp((${flowacc}*${res}),1.6000000)*exp(sin(${slope}),1.3000000)) * sin(${aspect})), d=10 * (${Kt} * exp(9806.65*(((${rain}/1000)*${flowacc})/(0.595*${stormtimet}))*tan(${slope}), ${loadexp}) ) * sin(${aspect}), if(${flowacc} >= ${cutoff3}, a, if(${flowacc} >= ${cutoff2} && ${flowacc} < ${cutoff3}, b, if(${flowacc} >= ${cutoff1} && ${flowacc} < ${cutoff2}, c, d))) )", quiet = "True", qsy = qsy, kappa = kappa, slope = slope, aspect = aspect, R = R, K = K, C =C, res = r, flowacc = flowacc, Kt = Kt, rain = rain, stormtimet = stormtimet, loadexp = loadexp, cutoff1 = cutoff1, cutoff2 = cutoff2, cutoff3 = cutoff3)
+ else: #do the stream powered version. Note that I'm converting the stream-power output (kg/m2) to same units as USPED (T/ha) by multiplying by ten. This ensures they are even going into the divergence calculation #Qs = Kt * n^-1 * 9810 * depth^1.6 * tan(slope)^1.5
+ grass.mapcalc("${qsx}=eval(a=(${kappa} * sin(${slope}) * cos(${aspect})), b=((${R}*${K}*${C}*${flowacc}*${res}*sin(${slope})) * cos(${aspect})), c=( (${R}*${K}*${C}*exp((${flowacc}*${res}),1.6000000)*exp(sin(${slope}),1.3000000)) * cos(${aspect})), d=10 *${Kt} * exp(${manningn}, -1) * 9810 * exp((((${rain}/1000)*${flowacc})/(0.595*${stormtimet})), 1.6) * exp(tan(${slope}), 1.5) * cos(${aspect}), if(${flowacc} <= ${cutoff1}, a, if(${flowacc} <= ${cutoff2} && ${flowacc} > ${cutoff1}, b, if(${flowacc} <= ${cutoff3} && ${flowacc} > ${cutoff2}, c, d))) )", quiet = "True", qsx = qsx, kappa = kappa, slope = slope, aspect = aspect, R = R, K = K, C =C, res = r, flowacc = flowacc, Kt = Kt, rain = rain, stormtimet = stormtimet, loadexp = loadexp, cutoff1 = cutoff1, cutoff2 = cutoff2, cutoff3 = cutoff3, manningn = manningn)
+ grass.mapcalc("${qsy}=eval(a=(${kappa} * sin(${slope}) * sin(${aspect})), b=((${R}*${K}*${C}*${flowacc}*${res}*sin(${slope})) * sin(${aspect})), c=( (${R}*${K}*${C}*exp((${flowacc}*${res}),1.6000000)*exp(sin(${slope}),1.3000000)) * sin(${aspect})), d=10 * ${Kt} * exp(${manningn}, -1) * 9810 * exp((((${rain}/1000)*${flowacc})/(0.595*${stormtimet})), 1.6) * exp(tan(${slope}), 1.5) * sin(${aspect}), if(${flowacc} <= ${cutoff1}, a, if(${flowacc} <= ${cutoff2} && ${flowacc} > ${cutoff1}, b, if(${flowacc} <= ${cutoff3} && ${flowacc} > ${cutoff2}, c, d))) )", quiet = "True", qsy = qsy, kappa = kappa, slope = slope, aspect = aspect, R = R, K = K, C =C, res = r, flowacc = flowacc, Kt = Kt, rain = rain, stormtimet = stormtimet, loadexp = loadexp, cutoff1 = cutoff1, cutoff2 = cutoff2, cutoff3 = cutoff3, manningn = manningn)
+ #make a map of the total TC for debugging purposes
+ #TC = "%sTC_%04d" % (p,o)
+ #grass.mapcalc("${TC}=eval(a=${kappa} * sin(${slope}), b=${R}*${K}*${C}*${flowacc}*${res}*sin(${slope}), c=${R}* ${K}* ${C}* exp( (${flowacc}*${res}),1.6000000) * exp(sin(${slope}),1.3000000), d=10 * ${Kt} * exp(${manningn}, -1) * 9810 * exp( ( ( (${rain}/1000)*${flowacc}) / (0.595*${stormtimet}) ), 1.6) * exp(tan(${slope}), 1.5), if(${flowacc} >= ${cutoff3}, a, if(${flowacc} >= ${cutoff2} && ${flowacc} < ${cutoff3}, b, if(${flowacc} >= ${cutoff1} && ${flowacc} < ${cutoff2}, c, d) ) ) )", quiet = "True", TC = TC, kappa = kappa, slope = slope, aspect = aspect, R = R, K = K, C =C, res = r, flowacc = flowacc, Kt = Kt, rain = rain, stormtimet = stormtimet, loadexp = loadexp, cutoff1 = cutoff1, cutoff2 = cutoff2, cutoff3 = cutoff3, manningn = manningn)
+ #/CHANGES
+
+ grass.message('\n*************************\n Iteration %s -- ' % o + 'step 4: calculating divergence/difference of sediment transport for each process and the actual amount of erosion or deposition in vertical meters/cell/year\n*************************\n\n')
+ #Here is where we figure out the change in transport capacity, and thus the actual amount of erosion an deposition that would occur. There are two ways of doing this. On planar and convex surfaces (i.e., ridgetops, flats, hillslopes), it is better to take the 2D divergence of sediment flux (we use r.slope.aspect to calculate this), but on highly convex surfaces (i.e., in channels) it is better to take the 1D difference between one cell, and the cell that is immediately downstream from it. This all assumes that the system is always operating at Transport Capacity, or if it is not, then is still behaves as if it were (ie., that the actual differences in transported sediment between the cells would be proportional to the system operating at capacity). Thus, under this assumption, the divergence of capacity is equals to actual amount of sediment eroded/deposited.
+ #This is the way we implemnt this: First calculate, we calculate the divergence/differnce for EACH of the different flow processes on the ENTIRE map (i.e., make one map per process, difference for streams, divergence for USPED and diffusion). Then, we cut out the pieces of each of these maps that correspond to the correct landforms from each specific process (based on the user-input cutoffs in flow accumulation), and patch them together into a single map (NOTE: see output unit conversions section below to see how we get all the units to line up during this process). This counters the "boundary effect" that happens when running the differential equations for divergence across the boundary of two different flow processes. Then we may still have to run a median smoother on the patched map to get rid of any latent spikes.
+ if ( flags["1"] is True ):
+ #This is the version with 1D streams
+ qsd1 = '%sDelta_Qs_1D_streams%04d' % (p, o)
+ grass.mapcalc('${qsd1}=if(${flowdir} == 7, (${qs1}[-1,-1]-${qs1}), if (${flowdir} == 6, (${qs1}[-1,0]-${qs1}), if (${flowdir} == 5, (${qs1}[-1,1]-${qs1}), if (${flowdir} == 4, (${qs1}[0,1]-${qs1}), if (${flowdir} == 3, (${qs1}[1,1]-${qs1}), if (${flowdir} == 2, (${qs1}[1,0]-${qs1}), if (${flowdir} == 1, (${qs1}[1,-1]-${qs1}), if (${flowdir} == 8, (${qs1}[0,-1]-${qs1}), ${qs1}))))))))', quiet = "True", qsd1 = qsd1, flowdir = flowdir, qs1 = qs1)
+ qsxdx = '%sDelta_Qsx_%04d' % (p, o)
+ qsydy = '%sDelta_Qsy_%04d' % (p, o)
+ grass.run_command('r.slope.aspect', quiet = "True", elevation = qsx, dx = qsxdx)
+ grass.run_command('r.slope.aspect', quiet = "True", elevation = qsy, dy = qsydy)
+ else:
+ #This is the normal version (with 2D streams)
+ qsxdx = '%sDelta_Qsx_%04d' % (p, o)
+ qsydy = '%sDelta_Qsy_%04d' % (p, o)
+ grass.run_command('r.slope.aspect', quiet = "True", elevation = qsx, dx = qsxdx)
+ grass.run_command('r.slope.aspect', quiet = "True", elevation = qsy, dy = qsydy)
+
+ #This is the smoothing routine. First we calculate the rate of Erosion and Deposition by converting the Delta QS of the different processes to vertical meters by dividing by the soil denisity (with apropriate constants to get into the correct units, see UNIT CONVERSION note below), and for streams, also expand from the storm to the year level. All units of this initial (temporary) ED_rate map will be in m/cell/year.
+ #CHANGES
+ #OUTPUT UNIT CONVERSIONS: In the case of the diffusion equation, the output units are in vertical meters of sediment per cell per year, so these will be left alone. Everything else should be in units of T/cell per storm. So we just need to convert to kg/cell, divide by the soil density and multiply the number of storms
+ if ( flags["1"] is True ):
+ #This is the version with 1D streams
+ grass.mapcalc('${tempnetchange1}=if(${flowacc} >= ${cutoff3}, ((${qsd1}*0.1)/${sdensity})*${storms}, if(${flowacc} >= ${cutoff1} && ${flowacc} < ${cutoff3}, (((${qsxdx}+${qsydy})*0.1)/${sdensity})*${storms}, ${qsxdx}+${qsydy}))', quiet = "True", tempnetchange1 = tempnetchange1, qsd1 = qsd1, qsxdx = qsxdx, qsydy = qsydy, flowacc = flowacc, cutoff1 = cutoff1, cutoff3 = cutoff3, sdensity = sdensity, storms = storms, stormtimet = stormtimet)
+ else:
+ #This is the normal version (with 2D streams)
+ grass.mapcalc('${tempnetchange1}=if(${flowacc} >= ${cutoff1}, (((${qsxdx} + ${qsydy})*0.1)/${sdensity})*${storms}, ${qsxdx}+${qsydy})', quiet = "True", tempnetchange1 = tempnetchange1, qsxdx = qsxdx, qsydy = qsydy, flowacc = flowacc, cutoff1 = cutoff1, cutoff3 = cutoff3, sdensity = sdensity, storms = storms, stormtimet = stormtimet)
+ #grass.mapcalc('${tempnetchange1}=if(${flowacc} >= ${cutoff3}, (((${qsxdx} + ${qsydy})*0.1)/${sdensity})*${storms}, if(${flowacc} >= ${cutoff1} && ${flowacc} < ${cutoff3}, ((${qsxdx}+${qsydy})*0.1)/${sdensity}, ${qsxdx}+${qsydy}))', quiet = "True", tempnetchange1 = tempnetchange1, qsxdx = qsxdx, qsydy = qsydy, flowacc = flowacc, cutoff1 = cutoff1, cutoff3 = cutoff3, sdensity = sdensity, storms = storms, stormtimet = stormtimet)
+ #/CHANGES
+
+ #Make some temp maps of just erosion rate and just deposition rate so we can grab some stats from them for the soft-knee limiting filter
+ grass.message('Running soft-knee smoothing filter...')
+ grass.mapcalc('${tmperosion}=if(${tempnetchange1} < -0, ${tempnetchange1}, null())', quiet = "True", tmperosion = tmperosion, tempnetchange1 = tempnetchange1)
+ grass.mapcalc('${tmpdep}=if(${tempnetchange1} > 0, ${tempnetchange1}, null())', quiet = "True", tmpdep = tmpdep, tempnetchange1 = tempnetchange1)
+ #Grab the stats from these temp files and save them to dictionaries
+ erosstats = grass.parse_command('r.univar', flags = 'ge', percentile = '1', map = tmperosion)
+ depostats = grass.parse_command('r.univar', flags = 'ge', percentile = '99', map = tmpdep)
+ maximum = depostats['max']
+ minimum = erosstats['min']
+ erosbreak = float(erosstats['first_quartile'])
+ deposbreak = float(depostats['third_quartile'])
+ scalemin = float(erosstats['percentile_1'])
+ scalemax = float(depostats['percentile_99'])
+ #Use the stats we gathered to do some smoothing with a hi-cut and lo-cut filter (with soft-knee limiting) of the unsmoothed ED_rate map. Values from the 1st quartile of erosion to the minimum (i.e., the very large negative numbers) will be rescaled linearly from the 1st quartile to the 1st percentile value, and values from the 3rd quartile of deposition to the maximum (i.e., the very large positiive numbers) will be rescaled linearly from the 3rd quartile to the 99th percentile value. This brings any values that were really unreasonnable as originally calculated (spikes) into the range of what the maximum values should be on a normally distrubuted dataset, but does so with out a "brick wall" style of limiting, which would make all values above some cutoff equal to a theoretical maximum. By setting both maximum cutoff point AND a "soft" scaling point, this "soft-knee" style of limiting sill retains some of the original scaling at the high ends, which allows for the smooth
ed value of very high cells to still be relatively higher than values in other cells that were also above the scaling cutoff, but were not originally as high as those very high cells.
+ grass.mapcalc('${tempnetchange2}=graph(${tempnetchange1}, ${minimum},${scalemin}, ${erosbreak},${erosbreak}, ${deposbreak},${deposbreak}, ${maximum},${scalemax})', quiet = "True", tempnetchange2 = tempnetchange2, tempnetchange1 =tempnetchange1, minimum = minimum, scalemin = scalemin, erosbreak = erosbreak, deposbreak = deposbreak, maximum = maximum, scalemax = scalemax)
+ #Check if additional smoothing is requested.
+ if smoothing == "no":
+ grass.message('No additional modal smoothing was requested...')
+ grass.run_command('g.rename', quiet = "True", rast = tempnetchange2 + ',' + netchange)
+ elif smoothing == "low":
+ grass.message('Enacting additional "low" smoothing: one pass of a 3x3 modal smoothing window.')
+ grass.run_command('r.neighbors', quiet = "True", input = tempnetchange2, output = netchange, method = 'mode', size = '3')
+ elif smoothing == "high":
+ grass.message('Enacting additional "high" smoothing: one pass of a 5x5 modal smoothing window.')
+ grass.run_command('r.neighbors', quiet = "True", input = tempnetchange2, output = netchange, method = 'mode', size = '5')
+ else:
+ grass.message('There was a problem reading the median-smoothing variable, so maps will not be median-smoothed.')
+ grass.run_command('g.rename', quiet = "True", rast = tempnetchange2 + ',' + netchange)
+ #Set the netchange map colors to the rules we've provided above
+ grass.run_command('r.colors', quiet = "True", map = netchange, rules = nccolors.name)
+ #Grab the stats from these new smoothed netchange maps and save them to dictionaries (Note that the temporary erosion and deposition maps made in this step are overwriting the two temporary maps made for gathering the stats for the soft-knee limiting filter)
+ grass.mapcalc('${tmperosion}=if(${netchange} < -0, ${netchange}, null())', quiet = "True", overwrite = "True", tmperosion = tmperosion, netchange = netchange)
+ grass.mapcalc('${tmpdep}=if(${netchange} > 0, ${netchange}, null())', quiet = "True", overwrite = "True", tmpdep = tmpdep, netchange = netchange)
+ erosstats1 = grass.parse_command('r.univar', flags = 'ge', map = tmperosion)
+ depostats1 = grass.parse_command('r.univar', flags = 'ge', map = tmpdep)
+
+ grass.message('\n*************************\n Iteration %s -- ' % o + 'step 5: calculating terrain evolution and new soil depths\n *************************\n\n')
+ #Set up a temp dem, and then do initial addition of ED change to old DEM. This mapcalc statement first checks the amount of erodable soil in a given cell against the amount of erosion calculated, and keeps the cell from eroding past this amount (if there is soil, then if the amount of erosion is more than the amount of soil, just remove all the soil and stop, else remove the amount of caclulated erosion. It also runs an error catch that checks to make sure that soil depth is not negative (could happen, I suppose), and if it is, corrects it). Finally, do patch-job to catch the shrinking edge problem (the edge cells have no upstream cell, so get turned null in the calculations in step 4)
+ grass.mapcalc('${new_dem}=eval(x=if(${old_soil} > 0.0 && (-1*${netchange}) <= ${old_soil}, ${netchange}, if((-1*${netchange}) > ${old_soil}, (-1*${old_soil}), 0)), y=(${old_dem} + x), if(isnull(y), ${old_dem}, y))', quiet = "True", new_dem = new_dem, old_soil = old_soil, old_dem = old_dem, netchange = netchange)
+ #Set colors for elevation map to match other dems
+ grass.run_command('r.colors', quiet = "True", map = new_dem, rast = options["elev"])
+ grass.mapcalc('${new_soil}=if ((${new_dem} - ${initbdrk}) < 0, 0, (${new_dem} - ${initbdrk}))', quiet = "True", new_soil = new_soil, new_dem = new_dem, initbdrk = initbdrk)
+ grass.run_command('r.colors', quiet = "True", map = new_soil, rules = sdcolors.name)
+ grass.message('\n*************************\n Iteration %s -- ' % o + 'step 6: writing stats to output file\n *************************\n\n')
+ #Finish gathering stats (just need the soil depth stats now)
+ soilstats = grass.parse_command('r.univar', flags = 'ge', map = new_soil, percentile = '99')
+ #Write stats to a new line in the stats file
+ #HEADER of the file should be: ',,Mean Values,,,,Standard Deviations,,,,Totals,,,Additional Stats\nIteration,,Mean Erosion,Mean Deposition,Mean Soil Depth,,Standard Deviation Erosion,Standard Deviation Deposition,Standard Deviation Soil Depth,,Total Sediment Eroded,Total Sediment Deposited,,Minimum Erosion,First Quartile Erosion,Median Erosion,Third Quartile Erosion,Maximum Erosion,Original Un-smoothed Maximum Erosion,,Minimum Deposition,First Quartile Deposition,Median Deposition,Third Quartile Deposition,Maximum Deposition,Original Un-smoothed Maximum Deposition,,Minimum Soil Depth,First Quartile Soil Depth,Median Soil Depth,Third Quartile Soil Depth,Maximum Soil Depth'
+ grass.message('Outputing stats to textfile: ' + q)
+ f.write('\n%s' % o + ',,' + erosstats1['mean'] + ',' + depostats1['mean'] + ',' + soilstats['mean'] + ',,' + erosstats1['stddev'] + ',' + depostats1['stddev'] + ',' + soilstats['stddev'] + ',,' + erosstats1['sum'] + ',' + depostats1['sum'] + ',,' + erosstats1['max'] + ',' + erosstats1['third_quartile'] + ',' + erosstats1['median'] + ',' + erosstats1['first_quartile'] + ',' + erosstats1['min'] + ',' + minimum + ',,' + depostats1['min'] + ',' + depostats1['first_quartile'] + ',' + depostats1['median'] + ',' + depostats1['third_quartile'] + ',' + depostats1['max'] + ',' + maximum + ',,' + soilstats['min'] + ',' + soilstats['first_quartile'] + ',' + soilstats['median'] + ',' + soilstats['third_quartile'] + ',' + soilstats['max'])
+
+ #Clean up temporary maps
+ if flags["k"] is True:
+ grass.message('\nTemporary maps will NOT be deleted!!!!\n')
+ else:
+ grass.message('\nCleaning up temporary maps...\n\n')
+ #first remove all the easy temporary maps labeled with "pid"
+ grass.run_command("g.remove", quiet = "True", flags = 'f', type = 'rast', pattern = '%s*' % pid)
+ #now check all the flag options, and build a list of maps to delete
+ mapstoremove = []
+ if flags["s"] is True:
+ grass.message('Keeping Slope map.')
+ else:
+ mapstoremove.append(slope)
+ if flags["d"] is True :
+ grass.message('Not keeping Soil Depth map.')
+ mapstoremove.append(old_soil)
+ #check if this is the last year and remove the "new-soil" map too
+ if ( o == int(options["number"])):
+ mapstoremove.append(new_soil)
+ else:
+ #check if this is the first year, and if so, remove the temporary "soildepths_init" map
+ if ( o <= 1 ):
+ mapstoremove.append("%s%s_init" % (prefx, options["outsoil"]))
+ if flags["e"] is True :
+ grass.message('Keeping Excess Transport Capacity (divergence) maps for all processes.')
+ else:
+ mapstoremove.extend([qsxdx, qsydy])
+ if flags["1"] is True :
+ mapstoremove.append(qsd1)
+ if flags["t"] is True :
+ grass.message('Keeping Transport Capacity maps for all processes.')
+ else:
+ mapstoremove.extend([qsx, qsy])
+ if flags["1"] is True :
+ mapstoremove.append(qs1)
+ if flags["r"] is True :
+ grass.message('Not keeping an Erosion and Deposition rate map.')
+ mapstoremove.append(netchange)
+ if len(mapstoremove) == 0:
+ pass
+ else:
+ grass.run_command('g.remove', quiet = "True", flags = 'f', type = "rast", name = ','.join(mapstoremove))
+ sdcolors.close()
+ nccolors.close()
+ grass.message('\n*************************\nDone with Iteration %s ' % o + '\n*************************\n')
+ return(0)
+
+#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__":
+ options, flags = grass.parser()
+ # Set up some basic variables
+ years = options["number"]
+ prefx = options["prefx"]
+ #these values could be read in from a climate file, so check that, and act accordingly. Either way, the result will be some lists with the same number of entries as there are iterations.
+ rain2 = []
+ try:
+ rain1 = float(options["rain"])
+ for year in range(int(years)):
+ rain2.append(rain1)
+ except:
+ with open(options["rain"], 'rU') as f:
+ for line in f:
+ rain2.append(line.split(",")[0])
+ #check for text header and remove if present
+ try:
+ float(rain2[0])
+ except:
+ del rain2[0]
+ #throw a warning if there aren't enough values in the column
+ if len(rain2) != int(years):
+ grass.fatal("Number of rows of rainfall data in your climate file\n do not match the number of iterations you wish to run.\n Please ensure that these numbers match and try again")
+ sys.exit(1)
+ R2 = []
+ try:
+ R1 = float(options["r"])
+ for year in range(int(years)):
+ R2.append(R1)
+ except:
+ with open(options["r"], 'rU') as f:
+ for line in f:
+ R2.append(line.split(",")[1])
+ #check for text header and remove if present
+ try:
+ float(R2[0])
+ except:
+ del R2[0]
+ #throw a warning if there aren't enough values in the column
+ if len(R2) != int(years):
+ grass.fatal("Number of rows of R-Factor data in your climate file\n do not match the number of iterations you wish to run.\n Please ensure that these numbers match and try again")
+ sys.exit(1)
+ storms2 = []
+ try:
+ storms1 = float(options["storms"])
+ for year in range(int(years)):
+ storms2.append(storms1)
+ except:
+ with open(options["storms"], 'rU') as f:
+ for line in f:
+ storms2.append(line.split(",")[2])
+ #check for text header and remove if present
+ try:
+ float(storms2[0])
+ except:
+ del storms2[0]
+ #throw a warning if there aren't enough values in the column
+ if len(storms2) != int(years):
+ grass.fatal("Number of rows of storm frequency data in your climate file\n do not match the number of iterations you wish to run.\n Please ensure that these numbers match and try again")
+ sys.exit(1)
+ stormlength2 = []
+ try:
+ stormlength1 = float(options["stormlength"])
+ for year in range(int(years)):
+ stormlength2.append(stormlength1)
+ except:
+ with open(options["stormlength"], 'rU') as f:
+ for line in f:
+ stormlength2.append(line.split(",")[3])
+ #check for text header and remove if present
+ try:
+ float(stormlength2[0])
+ except:
+ del stormlength2[0]
+ #throw a warning if there aren't enough values in the column
+ if len(stormlength2) != int(years):
+ grass.fatal("Number of rows of storm length data in your climate file\n do not match the number of iterations you wish to run.\n Please ensure that these numbers match and try again")
+ sys.exit(1)
+ #Now gather these four lists into one master list, to make it easier to pass on to main()
+ masterlist = [rain2,R2,storms2,stormlength2]
+ #Make the statsout file with correct column headers
+ if options["statsout"] == "":
+ env = grass.gisenv()
+ mapset = env['MAPSET']
+ statsout = '%s_%slsevol_stats.csv' % (mapset, prefx)
+ else:
+ statsout = options["statsout"]
+ if os.path.isfile(statsout):
+ f = file(statsout, 'a')
+ else:
+ f = file(statsout, 'wt')
+ f.write('These statistics are in units of vertical meters (depth) per cell\n,,Mean Values,,,,Standard Deviations,,,,Totals,,,Additional Stats\nIteration,,Mean Erosion,Mean Deposition,Mean Soil Depth,,Standard Deviation Erosion,Standard Deviation Deposition,Standard Deviation Soil Depth,,Total Sediment Eroded,Total Sediment Deposited,,Minimum Erosion,First Quartile Erosion,Median Erosion,Third Quartile Erosion,Maximum Erosion,Original Un-smoothed Maximum Erosion,,Minimum Deposition,First Quartile Deposition,Median Deposition,Third Quartile Deposition,Maximum Deposition,Original Un-smoothed Maximum Deposition,,Minimum Soil Depth,First Quartile Soil Depth,Median Soil Depth,Third Quartile Soil Depth,Maximum Soil Depth')
+ if flags["p"] is True :
+ grass.message('Making sample points map for determining cutoffs.')
+ else:
+ grass.message('\n##################################################\n##################################################\n\n STARTING SIMULATION\n\nBeginning iteration sequence. This may take some time.\nProcess is not finished until you see the message: \'Done with everything\'\n _____________________________________________________________\n_____________________________________________________________\n')
+ grass.message("Total number of iterations to be run is %s" % years)
+ #Get the region settings
+ region1 = grass.region()
+ # This is the loop!
+ for x in range(int(years)):
+ grass.message("Iteration = %s" % (x + 1))
+ main(x, (x + 1), prefx, statsout, region1['nsres'], masterlist);
+ #Since we are now done with the loop, close the stats file.
+ f.close()
+ grass.message('\nIterations complete!\n\nDone with everything')
+ sys.exit(0)
+
+
+
+
Copied: grass-addons/grass7/raster/r.landscape.evol/r.landscape.evol_description.odf.odt (from rev 64600, grass-addons/grass7/raster/r.landscape.evol/r.landscape.evol description.odf.odt)
===================================================================
(Binary files differ)
Copied: grass-addons/grass7/raster/r.landscape.evol/r.landscape.evol_description.odt (from rev 64600, grass-addons/grass7/raster/r.landscape.evol/r.landscape.evol description.odt)
===================================================================
(Binary files differ)
Deleted: grass-addons/grass7/raster/r.viewshed.cva/r.viewshed.cva
===================================================================
--- grass-addons/grass7/raster/r.viewshed.cva/r.viewshed.cva 2015-02-13 03:10:26 UTC (rev 64600)
+++ grass-addons/grass7/raster/r.viewshed.cva/r.viewshed.cva 2015-02-13 09:49:02 UTC (rev 64601)
@@ -1,183 +0,0 @@
-#!/usr/bin/python
-#
-############################################################################
-#
-# MODULE: r.viewshed.cva.py
-# AUTHOR(S): Isaac Ullah
-# PURPOSE: Undertakes a "cumulative viewshed analysis" using a vector points map as input "viewing" locations, using r.viewshed to calculate the individual viewsheds.
-# COPYRIGHT: (C) 2015 by Isaac Ullah
-# REFERENCES: r.viewshed
-# This program is free software under the GNU General Public
-# License (>=v2). Read the file COPYING that comes with GRASS
-# for details.
-#
-#############################################################################
-
-
-#%Module
-#% description: Undertakes a "cumulative viewshed analysis" using a vector points map as input "viewing" locations, using r.viewshed to calculate the individual viewsheds.
-#%End
-
-#%option
-#% key: elev
-#% type: string
-#% gisprompt: new,cell,raster
-#% description: Input elevation map (DEM)
-#% required : yes
-#%END
-
-#%option
-#% key: output
-#% type: string
-#% gisprompt: new,cell,raster
-#% description: Output CVA raster
-#% required : yes
-#%END
-
-#%option
-#% key: vect
-#% type: string
-#% gisprompt: new,vector,vector
-#% description: Name of input vector points map containg the set of sites for this analysis.
-#% required : yes
-#%END
-
-#%option
-#% key: observer_elevation
-#% type: string
-#% description: Height of observation points off the ground
-#%answer: 0.0
-#% required : yes
-#%END
-
-#%option
-#% key: target_elevation
-#% type: string
-#% description: Height of target areas off the ground
-#%answer: 1.75
-#% required : yes
-#%END
-
-#%option
-#% key: max_distance
-#% type: string
-#% description: Maximum viewing distance (-1 = infinity)
-#%answer: -1
-#% required : yes
-#%END
-
-#%option
-#% key: memory
-#% type: string
-#% description: Amount of memory to use (in MB)
-#%answer: 1500
-#% required : yes
-#%END
-
-#%option
-#% key: refraction_coeff
-#% type: string
-#% description: Refraction coefficient (with flag -r)
-#%answer: 0.14286
-#% required : no
-#%END
-
-
-#%flag
-#% key: k
-#% description: -k Keep all interim viewshed maps produced by the routine (maps will be named "vshed_'name'", where 'name' is the value in "name_column" for each input point)
-#%END
-
-#%flag
-#% key: c
-#% description: -c Consider the curvature of the earth (current ellipsoid)
-#%END
-
-#%flag
-#% key: r
-#% description: -r Consider the effect of atmospheric refraction
-#%END
-
-#%flag
-#% key: b
-#% description: -b Output format is {0 (invisible) 1 (visible)}
-#%END
-
-#%flag
-#% key: e
-#% description: -e Output format is invisible = NULL, else current elev - viewpoint_elev
-#%END
-
-
-
-import sys
-import os
-grass_install_tree = os.getenv('GISBASE')
-sys.path.append(grass_install_tree + os.sep + 'etc' + os.sep + 'python')
-import grass.script as grass
-
-
-#main block of code starts here
-def main():
- #bring in input variables
- elev = options["elev"]
- vect = options["vect"]
- observer_elevation =options["observer_elevation"]
- target_elevation = options['target_elevation']
- max_distance = options["max_distance"]
- memory = options["memory"]
- refraction_coeff = options["refraction_coeff"]
- out = options["output"]
- #assemble flag string
- if flags['r'] is True:
- f1 = "r"
- else:
- f1 = ""
- if flags['c'] is True:
- f2 = "c"
- else:
- f2 = ""
- if flags['b'] is True:
- f3 = "b"
- else:
- f3 = ""
- if flags['e'] is True:
- f4 = "e"
- else:
- f4 = ""
- flagstring = f1 + f2 + f3 +f4
- #make a tempfile, and write out the coords from the vector map.
- tmp1 = grass.tempfile()
- grass.run_command("v.out.ascii", flags = 'r', input = vect, type = "point", output = tmp1, format = "point", separator = ",")
- # note that the "r" flag will constrain to points in the current geographic region.
- grass.message("Note that the routine is constrained to points in the current geographic region.")
- #read the temp file back in, and parse it up.
- f = open(tmp1, 'r')
- masterlist = []
- for line in f.readlines():
- masterlist.append(line.strip('\n').split(','))
- f.close() #close the file
- #now, loop through the master list and run r.viewshed for each of the sites, and append the viewsheds to a list (so we can work with them later)
- vshed_list = []
- for site in masterlist:
- grass.message('Calculating viewshed for location %s,%s (point name = %s)\n' % (site[0], site[1], site[2]))
- tempry = "vshed_%s" % site[2]
- vshed_list.append(tempry)
- grass.run_command("r.viewshed", quiet = "True", overwrite = grass.overwrite(), flags = flagstring, input = elev, output = tempry, coordinates = site[0] + "," + site[1], observer_elevation = observer_elevation, target_elevation = target_elevation, max_distance = max_distance, memory = memory, refraction_coeff = refraction_coeff)
- #now make a mapcalc statement to add all the viewsheds together to make the outout cumulative viewsheds map
- grass.message("Calculating \"Cumulative Viewshed\" map")
- #grass.mapcalc("${output}=${command_string}", quiet = "True", output = out, command_string = ("+").join(vshed_list))
- grass.run_command("r.series", quiet = "True", overwrite = grass.overwrite(), input = (",").join(vshed_list), output = out, method = "count")
- #Clean up temporary maps, if requested
- if os.getenv('GIS_FLAG_k') == '1':
- grass.message("Temporary viewshed maps will not removed")
- else:
- grass.message("Removing temporary viewshed maps")
- grass.run_command("g.remove", quiet = "True", flags = 'f', type = 'raster', name = (",").join(vshed_list))
- 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__":
- options, flags = grass.parser()
- main()
- exit(0)
Copied: grass-addons/grass7/raster/r.viewshed.cva/r.viewshed.cva.py (from rev 64600, grass-addons/grass7/raster/r.viewshed.cva/r.viewshed.cva)
===================================================================
--- grass-addons/grass7/raster/r.viewshed.cva/r.viewshed.cva.py (rev 0)
+++ grass-addons/grass7/raster/r.viewshed.cva/r.viewshed.cva.py 2015-02-13 09:49:02 UTC (rev 64601)
@@ -0,0 +1,183 @@
+#!/usr/bin/python
+#
+############################################################################
+#
+# MODULE: r.viewshed.cva.py
+# AUTHOR(S): Isaac Ullah
+# PURPOSE: Undertakes a "cumulative viewshed analysis" using a vector points map as input "viewing" locations, using r.viewshed to calculate the individual viewsheds.
+# COPYRIGHT: (C) 2015 by Isaac Ullah
+# REFERENCES: r.viewshed
+# This program is free software under the GNU General Public
+# License (>=v2). Read the file COPYING that comes with GRASS
+# for details.
+#
+#############################################################################
+
+
+#%Module
+#% description: Undertakes a "cumulative viewshed analysis" using a vector points map as input "viewing" locations, using r.viewshed to calculate the individual viewsheds.
+#%End
+
+#%option
+#% key: elev
+#% type: string
+#% gisprompt: new,cell,raster
+#% description: Input elevation map (DEM)
+#% required : yes
+#%END
+
+#%option
+#% key: output
+#% type: string
+#% gisprompt: new,cell,raster
+#% description: Output CVA raster
+#% required : yes
+#%END
+
+#%option
+#% key: vect
+#% type: string
+#% gisprompt: new,vector,vector
+#% description: Name of input vector points map containg the set of sites for this analysis.
+#% required : yes
+#%END
+
+#%option
+#% key: observer_elevation
+#% type: string
+#% description: Height of observation points off the ground
+#%answer: 0.0
+#% required : yes
+#%END
+
+#%option
+#% key: target_elevation
+#% type: string
+#% description: Height of target areas off the ground
+#%answer: 1.75
+#% required : yes
+#%END
+
+#%option
+#% key: max_distance
+#% type: string
+#% description: Maximum viewing distance (-1 = infinity)
+#%answer: -1
+#% required : yes
+#%END
+
+#%option
+#% key: memory
+#% type: string
+#% description: Amount of memory to use (in MB)
+#%answer: 1500
+#% required : yes
+#%END
+
+#%option
+#% key: refraction_coeff
+#% type: string
+#% description: Refraction coefficient (with flag -r)
+#%answer: 0.14286
+#% required : no
+#%END
+
+
+#%flag
+#% key: k
+#% description: -k Keep all interim viewshed maps produced by the routine (maps will be named "vshed_'name'", where 'name' is the value in "name_column" for each input point)
+#%END
+
+#%flag
+#% key: c
+#% description: -c Consider the curvature of the earth (current ellipsoid)
+#%END
+
+#%flag
+#% key: r
+#% description: -r Consider the effect of atmospheric refraction
+#%END
+
+#%flag
+#% key: b
+#% description: -b Output format is {0 (invisible) 1 (visible)}
+#%END
+
+#%flag
+#% key: e
+#% description: -e Output format is invisible = NULL, else current elev - viewpoint_elev
+#%END
+
+
+
+import sys
+import os
+grass_install_tree = os.getenv('GISBASE')
+sys.path.append(grass_install_tree + os.sep + 'etc' + os.sep + 'python')
+import grass.script as grass
+
+
+#main block of code starts here
+def main():
+ #bring in input variables
+ elev = options["elev"]
+ vect = options["vect"]
+ observer_elevation =options["observer_elevation"]
+ target_elevation = options['target_elevation']
+ max_distance = options["max_distance"]
+ memory = options["memory"]
+ refraction_coeff = options["refraction_coeff"]
+ out = options["output"]
+ #assemble flag string
+ if flags['r'] is True:
+ f1 = "r"
+ else:
+ f1 = ""
+ if flags['c'] is True:
+ f2 = "c"
+ else:
+ f2 = ""
+ if flags['b'] is True:
+ f3 = "b"
+ else:
+ f3 = ""
+ if flags['e'] is True:
+ f4 = "e"
+ else:
+ f4 = ""
+ flagstring = f1 + f2 + f3 +f4
+ #make a tempfile, and write out the coords from the vector map.
+ tmp1 = grass.tempfile()
+ grass.run_command("v.out.ascii", flags = 'r', input = vect, type = "point", output = tmp1, format = "point", separator = ",")
+ # note that the "r" flag will constrain to points in the current geographic region.
+ grass.message("Note that the routine is constrained to points in the current geographic region.")
+ #read the temp file back in, and parse it up.
+ f = open(tmp1, 'r')
+ masterlist = []
+ for line in f.readlines():
+ masterlist.append(line.strip('\n').split(','))
+ f.close() #close the file
+ #now, loop through the master list and run r.viewshed for each of the sites, and append the viewsheds to a list (so we can work with them later)
+ vshed_list = []
+ for site in masterlist:
+ grass.message('Calculating viewshed for location %s,%s (point name = %s)\n' % (site[0], site[1], site[2]))
+ tempry = "vshed_%s" % site[2]
+ vshed_list.append(tempry)
+ grass.run_command("r.viewshed", quiet = "True", overwrite = grass.overwrite(), flags = flagstring, input = elev, output = tempry, coordinates = site[0] + "," + site[1], observer_elevation = observer_elevation, target_elevation = target_elevation, max_distance = max_distance, memory = memory, refraction_coeff = refraction_coeff)
+ #now make a mapcalc statement to add all the viewsheds together to make the outout cumulative viewsheds map
+ grass.message("Calculating \"Cumulative Viewshed\" map")
+ #grass.mapcalc("${output}=${command_string}", quiet = "True", output = out, command_string = ("+").join(vshed_list))
+ grass.run_command("r.series", quiet = "True", overwrite = grass.overwrite(), input = (",").join(vshed_list), output = out, method = "count")
+ #Clean up temporary maps, if requested
+ if os.getenv('GIS_FLAG_k') == '1':
+ grass.message("Temporary viewshed maps will not removed")
+ else:
+ grass.message("Removing temporary viewshed maps")
+ grass.run_command("g.remove", quiet = "True", flags = 'f', type = 'raster', name = (",").join(vshed_list))
+ 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__":
+ options, flags = grass.parser()
+ main()
+ exit(0)
More information about the grass-commit
mailing list