[GRASS-SVN] r64448 - in grass-addons/grass7/raster: . r.flexure
svn_grass at osgeo.org
svn_grass at osgeo.org
Tue Feb 3 15:48:12 PST 2015
Author: awickert
Date: 2015-02-03 15:48:11 -0800 (Tue, 03 Feb 2015)
New Revision: 64448
Added:
grass-addons/grass7/raster/r.flexure/
grass-addons/grass7/raster/r.flexure/Makefile
grass-addons/grass7/raster/r.flexure/r.flexure.html
grass-addons/grass7/raster/r.flexure/r.flexure.py
Log:
Adding flexural isostasy raster interface
Added: grass-addons/grass7/raster/r.flexure/Makefile
===================================================================
--- grass-addons/grass7/raster/r.flexure/Makefile (rev 0)
+++ grass-addons/grass7/raster/r.flexure/Makefile 2015-02-03 23:48:11 UTC (rev 64448)
@@ -0,0 +1,7 @@
+MODULE_TOPDIR = ../..
+
+PGM = r.flexure
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script
Added: grass-addons/grass7/raster/r.flexure/r.flexure.html
===================================================================
--- grass-addons/grass7/raster/r.flexure/r.flexure.html (rev 0)
+++ grass-addons/grass7/raster/r.flexure/r.flexure.html 2015-02-03 23:48:11 UTC (rev 64448)
@@ -0,0 +1,60 @@
+<h2>DESCRIPTION</h2>
+
+<em>r.flexure</em> computes how the rigid outer shell of a planet deforms elastically in response to surface-normal loads by solving equations for plate bending. This phenomenon is known as "flexural isostasy" and can be useful in cases of glacier/ice-cap/ice-sheet loading, sedimentary basin filling, mountain belt growth, volcano emplacement, sea-level change, and other geologic processes. <em>r.flexure</em> and <a href="v.flexure"><em>v.flexure</em></a> are the GRASS GIS interfaces to the the model <a href=http://csdms.colorado.edu/wiki/Model:gFlex><b>gFlex</b></a>.
+
+As both <em>r.flexure</em> and <a href="v.flexure"><em>v.flexure</em></a> are interfaces to gFlex, this must be downloaded and installed. The most recent versions of <b>gFlex</b> are available from <a href=https://github.com/awickert/gFlex>https://github.com/awickert/gFlex</a>, and installation instructions are avaliable on that page via the <i>README.md</i> file.
+
+<h2>NOTES</h2>
+
+The parameter <b>method</b> sets whether the solution is Finite Difference ("FD") or Superposition of Analytical Solutions ("SAS"). The Finite difference method is typically faster for large arrays, and allows lithospheric elastic thickness to be varied laterally, following the solution of van Wees and Cloetingh (1994). However, it is quite memory-intensive, so unless the user has a computer with a very large amount of memory and quite a lot of time to wait, they should ensure that they use a grid spacing that is appropriate to solve the problem at hand. Flexural isostatic solutions act to smooth inputs over a given flexural wavelength (see , so if an appropriate solution resolution is chosen, the calculated flexural response can be interpolated to a higher resolution without fear of aliasing.
+<p>
+The flexural solution is generated for the current computational region, so be sure to check <em>g.region</em> before running the model!
+<p>
+<!--Add equation or a link to paper, etc.-->
+<b>q0</b> is a 2-D array of loads in a GRASS raster. These are in units of stress, and equal the density of the material times the acceleration due to gravity times the thickness of the column. This is not affected by what you choose for <b>g</b>, later: it is pre-calculated by the user.
+<p>
+<b>te</b>, written in standard text as T<sub>e</sub> is the lithospheric elastic thickness.
+<p>
+Several boundary conditions are available, and these depend on if the solution method is finite differece (FD) or superposition of analytical solutions (SAS). In the latter, it is assumed that there are no loads outside of those that are explicitly listed, so the boundary conditions are "NoOutsideLoads". As this is the implicit case, the boundary conditions all default to this.
+<p>
+The finite difference boundary conditions are a bit more complicated, but are largely self-explanitory:
+<p>
+ <dl>
+ <dt><b>Dirichlet0</b></dt>
+ <dd>0-displacement boundary condition</dd>
+ <dt><b>0Moment0Shear</b></dt>
+ <dd>"Broken plate" boundary condition: second and third derivatives of vertical displacement are 0. This is like the end of a diving board.</dd>
+ <dt><b>0Slope0Shear</b></dt>
+ <dd>First and third derivatives of vertical displacement are zero. While this does not lend itsellf so easily to physical meaning, it is helpful to aid in efforts to make boundary condition effects disappear (i.e. to emulate the NoOutsideLoads cases)</dd>
+ <dt><b>Mirror</b></dt>
+ <dd>Load and elastic thickness structures reflected at boundary.</dd>
+ <dt><b>Periodic</b></dt>
+ <dd>"Wrap-around" boundary condition: must be applied to both North and South and/or both East and West. This causes, for example, the edge of the eastern and western limits of the domain to act like they are next to each other in an infinite loop.</dd>
+</dl>
+<p>
+All of these boundary conditions may be combined in an way, with the exception of the note for periodic boundary conditions. If one does not want the boundary conditions to affect the solutions, it is recommended that one places the boundaries at lesat one flexural wavelength away from the load.
+<p>
+<em>r.flexure</em> may be run in latitude/longitude coordinates (with the "-l" flag), but its grid constraint is that it can have only one <i>dx</i> and one <i>dy</i> for the entire domain. Thus, it chooses the average <i>dx</i> at the midpoint between the northernmost and southernmost latitudes for which the calculations are made. This assumption can break down at the poles, where the East–West dimension rapidly diminishes.
+<p>
+<!--Add figures and example from paper once that is out-->
+<!--And examples-->
+<!--Add runtime fig.-->
+The <a href="http://csdms.colorado.edu">Community Surface Dynamics Modeling System</a>, into which <b>gFlex</b> is integrated, is a community-driven effort to build an open-source modeling infrastructure for Earth-surface processes.
+
+<h2>SEE ALSO</h2>
+
+<em>
+<a href="v.flexure">v.flexure</a>
+</em>
+
+<h2>REFERENCES</h2>
+
+Wickert, A. D., G. E. Tucker, E. W. H. Hutton, B. Yan, and S. D. Peckham (2011), <a href="http://csdms.colorado.edu/wiki/CSDMS_2011_annual_meeting_Andrew_Wickert">Feedbacks between surface processes and flexural isostasy: a motivation for coupling models</a>, in <i>CSDMS 2011 Meeting: Impact of time and process scales</i>, Student Keynote, Boulder, CO.
+<p>
+Wees, J. D., Cloetingh, S. and van Wees, J. D.: A Finite-Difference Technique to Incorporate Spatial Variations In Rigidity and Planar Faults Into 3-D Models For Lithospheric Flexure, Geophys. J. Int., 117(1), 179–195, doi:10.1111/j.1365-246X.1994.tb03311.x, 1994.
+
+<h2>AUTHOR</h2>
+
+Andrew D. Wickert<br>
+
+<p><i>Last changed: $Date 2015-02-04 00:29 +0100 (Wed, 4 Feb 2015)$</i>
Added: grass-addons/grass7/raster/r.flexure/r.flexure.py
===================================================================
--- grass-addons/grass7/raster/r.flexure/r.flexure.py (rev 0)
+++ grass-addons/grass7/raster/r.flexure/r.flexure.py 2015-02-03 23:48:11 UTC (rev 64448)
@@ -0,0 +1,293 @@
+#! /usr/bin/python
+############################################################################
+#
+# MODULE: r.flexure
+#
+# AUTHOR(S): Andrew Wickert
+#
+# PURPOSE: Calculate flexure of the lithosphere under a specified
+# set of loads and with a given elastic thickness (scalar
+# or array)
+#
+# COPYRIGHT: (c) 2012, 2014, 2015 Andrew Wickert
+#
+# This program is free software under the GNU General Public
+# License (>=v2). Read the file COPYING that comes with GRASS
+# for details.
+#
+#############################################################################
+#
+# REQUIREMENTS:
+# - gFlex: http://csdms.colorado.edu/wiki/gFlex
+# (should be downloaded automatically along with the module)
+# github repository: https://github.com/awickert/gFlex
+
+# More information
+# Started 11 March 2012 as a GRASS interface for Flexure (now gFlex)
+# Revised 15--?? November 2014 after significantly improving the model
+# by Andy Wickert
+
+#%module
+#% description: Computes lithospheric flexural isostasy
+#% keyword: raster
+#% keyword: geophysics
+#%end
+#%flag
+#% key: l
+#% description: Allows running in lat/lon: dx is f(lat) at grid N-S midpoint
+#%end
+#%option
+#% key: method
+#% type: string
+#% description: Solution method: Finite Diff. or Superpos. of analytical sol'ns
+#% options: FD, SAS
+#% required : yes
+#%end
+#%option
+#% key: q0
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: Raster map of loads (thickness * density * g) [Pa]
+#% required : yes
+#%end
+#%option
+#% key: te
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: Elastic thicnkess: scalar or raster; unis chosen in "te_units"
+#% required : yes
+#%end
+#%option
+#% key: te_units
+#% type: string
+#% description: Units for elastic thickness
+#% options: m, km
+#% required : yes
+#%end
+#%option
+#% key: output
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: Output raster map of vertical deflections [m]
+#% required : yes
+#%end
+#%option
+#% key: solver
+#% type: string
+#% description: Solver type
+#% options: direct, iterative
+#% answer: direct
+#% required : no
+#%end
+#%option
+#% key: tolerance
+#% type: double
+#% description: Convergence tolerance (between iterations) for iterative solver
+#% answer: 1E-3
+#% required : no
+#%end
+#%option
+#% key: northbc
+#% type: string
+#% description: Northern boundary condition
+#% options: Dirichlet0, 0Moment0Shear, 0Slope0Shear, Mirror, Periodic, NoOutsideLoads
+#% answer: NoOutsideLoads
+#% required : no
+#%end
+#%option
+#% key: southbc
+#% type: string
+#% description: Southern boundary condition
+#% options: Dirichlet0, 0Moment0Shear, 0Slope0Shear, Mirror, Periodic, NoOutsideLoads
+#% answer: NoOutsideLoads
+#% required : no
+#%end
+#%option
+#% key: westbc
+#% type: string
+#% description: Western boundary condition
+#% options: Dirichlet0, 0Moment0Shear, 0Slope0Shear, Mirror, Periodic, NoOutsideLoads
+#% answer: NoOutsideLoads
+#% required : no
+#%end
+#%option
+#% key: eastbc
+#% type: string
+#% description: Eastern boundary condition
+#% options: Dirichlet0, 0Moment0Shear, 0Slope0Shear, Mirror, Periodic, NoOutsideLoads
+#% answer: NoOutsideLoads
+#% required : no
+#%end
+#%option
+#% key: g
+#% type: double
+#% description: gravitational acceleration at surface [m/s^2]
+#% answer: 9.8
+#% required : no
+#%end
+#%option
+#% key: ym
+#% type: double
+#% description: Young's Modulus [Pa]
+#% answer: 65E9
+#% required : no
+#%end
+#%option
+#% key: nu
+#% type: double
+#% description: Poisson's ratio
+#% answer: 0.25
+#% required : no
+#%end
+#%option
+#% key: rho_fill
+#% type: double
+#% description: Density of material that fills flexural depressions [kg/m^3]
+#% answer: 0
+#% required : no
+#%end
+#%option
+#% key: rho_m
+#% type: double
+#% description: Mantle density [kg/m^3]
+#% answer: 3300
+#% required : no
+#%end
+
+##################
+# IMPORT MODULES #
+##################
+
+# PYTHON
+import numpy as np
+# GRASS
+import grass.script as grass
+import grass.script.array as garray
+# GFLEX
+try:
+ import gflex
+except:
+ print ""
+ print "MODULE IMPORT ERROR."
+ print "In order to run r.flexure or g.flexure, you must download and install"
+ print "gFlex. The most recent development version is available from"
+ print "https://github.com/awickert/gFlex."
+ print "Installation instructions are available on the page."
+ grass.fatal("Software dependency must be installed.")
+
+
+############################
+# PASS VARIABLES AND SOLVE #
+############################
+
+def main():
+ """
+ Gridded flexural isostatic solutions
+ """
+
+ # This code is for 2D flexural isostasy
+ flex = gflex.F2D()
+ # And show that it is coming from GRASS GIS
+ flex.grass = True
+
+ # Flags
+ latlon_override = flags['l']
+
+ # Inputs
+ # Solution selection
+ flex.Method = options['method']
+ if flex.Method == 'FD':
+ flex.Solver = options['solver']
+ if flex.Solver:
+ flex.ConvergenceTolerance = options['tolerance']
+ # Always use the van Wees and Cloetingh (1994) solution type.
+ # It is the best.
+ flex.PlateSolutionType = 'vWC1994'
+ # Parameters that are often changed for the solution
+ qs = options['q0']
+ flex.qs = garray.array()
+ flex.qs.read(qs)
+ # Elastic thickness
+ try:
+ flex.Te = float(options['te'])
+ except:
+ flex.Te = garray.array() # FlexureTe is the one that is used by Flexure
+ flex.Te.read(options['te'])
+ flex.Te = np.array(flex.Te)
+ if options['te_units'] == 'km':
+ flex.Te *= 1000
+ elif options['te_units'] == 'm':
+ pass
+ # No "else"; shouldn't happen
+ flex.rho_fill = float(options['rho_fill'])
+ # Parameters that often stay at their default values
+ flex.g = float(options['g'])
+ flex.E = float(options['ym']) # Can't just use "E" because reserved for "east", I think
+ flex.nu = float(options['nu'])
+ flex.rho_m = float(options['rho_m'])
+ # Solver type and iteration tolerance
+ flex.Solver = options['solver']
+ flex.ConvergenceTolerance = float(options['tolerance'])
+ # Boundary conditions
+ flex.BC_N = options['northbc']
+ flex.BC_S = options['southbc']
+ flex.BC_W = options['westbc']
+ flex.BC_E = options['eastbc']
+
+ # Set verbosity
+ if grass.verbosity() >= 2:
+ flex.Verbose = True
+ if grass.verbosity() >= 3:
+ flex.Debug = True
+ elif grass.verbosity() == 0:
+ flex.Quiet = True
+
+ # First check if output exists
+ if len(grass.parse_command('g.list', type='rast', pattern=options['output'])):
+ if not grass.overwrite():
+ grass.fatal("Raster map '" + options['output'] + "' already exists. Use '--o' to overwrite.")
+
+ # Get grid spacing from GRASS
+ # Check if lat/lon and proceed as directed
+ if grass.region_env()[6] == '3':
+ if latlon_override:
+ if flex.Verbose:
+ print "Latitude/longitude grid."
+ print "Based on r_Earth = 6371 km"
+ print "Setting y-resolution [m] to 111,195 * [degrees]"
+ flex.dy = grass.region()['nsres']*111195.
+ NSmid = (grass.region()['n'] + grass.region()['s'])/2.
+ dx_at_mid_latitude = (3.14159/180.) * 6371000. * np.cos(np.deg2rad(NSmid))
+ if flex.Verbose:
+ print "Setting x-resolution [m] to "+"%.2f" %dx_at_mid_latitude+" * [degrees]"
+ flex.dx = grass.region()['ewres']*dx_at_mid_latitude
+ else:
+ grass.fatal("Need the '-l' flag to enable lat/lon solution approximation.")
+ # Otherwise straightforward
+ else:
+ flex.dx = grass.region()['ewres']
+ flex.dy = grass.region()['nsres']
+
+ # CALCULATE!
+ flex.initialize()
+ flex.run()
+ flex.finalize()
+
+ # Write to GRASS
+ # Create a new garray buffer and write to it
+ outbuffer = garray.array() # Instantiate output buffer
+ outbuffer[...] = flex.w
+ outbuffer.write(options['output'], overwrite=grass.overwrite) # Write it with the desired name
+ # And create a nice colormap!
+ grass.run_command('r.colors', map=options['output'], color='differences', quiet=True)
+
+ # Reinstate this with a flag or output filename
+ # But I think better to let interpolation happen a posteriori
+ # So the user knows what the solution is and what it isn't
+ #grass.run_command('r.resamp.interp', input=output, output=output + '_interp', method='lanczos', overwrite=True, quiet=True)
+ #grass.run_command('r.colors', map=output + '_interp', color='rainbow', quiet=True)#, flags='e')
+
+if __name__ == "__main__":
+ options, flags = grass.parser()
+ main()
+
More information about the grass-commit
mailing list