[GRASS-SVN] r64449 - in grass-addons/grass7/vector: . v.flexure
svn_grass at osgeo.org
svn_grass at osgeo.org
Tue Feb 3 15:58:16 PST 2015
Author: awickert
Date: 2015-02-03 15:58:16 -0800 (Tue, 03 Feb 2015)
New Revision: 64449
Added:
grass-addons/grass7/vector/v.flexure/
grass-addons/grass7/vector/v.flexure/Makefile
grass-addons/grass7/vector/v.flexure/v.flexure.html
grass-addons/grass7/vector/v.flexure/v.flexure.py
Log:
Vector interface to the open-source lithospheric flexure model, gFlex
Added: grass-addons/grass7/vector/v.flexure/Makefile
===================================================================
--- grass-addons/grass7/vector/v.flexure/Makefile (rev 0)
+++ grass-addons/grass7/vector/v.flexure/Makefile 2015-02-03 23:58:16 UTC (rev 64449)
@@ -0,0 +1,7 @@
+MODULE_TOPDIR = ../..
+
+PGM = v.flexure
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script
Added: grass-addons/grass7/vector/v.flexure/v.flexure.html
===================================================================
--- grass-addons/grass7/vector/v.flexure/v.flexure.html (rev 0)
+++ grass-addons/grass7/vector/v.flexure/v.flexure.html 2015-02-03 23:58:16 UTC (rev 64449)
@@ -0,0 +1,33 @@
+<h2>DESCRIPTION</h2>
+
+<em>v.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>v.flexure</em> and <a href="r.flexure"><em>r.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>v.flexure</em> and <a href="r.flexure"><em>r.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 <em>README.md</em> file.
+
+<h2>NOTES</h2>
+
+<b>q0</b> is a vector points file containing the loads in units of force. Typically, this will be a representation of a distributed field of loads as a set of points, so the user will implicitly include the area over which a stress (vertical load) acts into the quantities in the database table of <b>q0</b>.
+<p>
+<b>te</b>, written in standard text as T<sub>e</sub> is the lithospheric elastic thickness.
+<p>
+<b>output</b> is provided as a grid of vector points corresponding to the GRASS region when this command is invoked. Be sure to use <em>g.region</em> to properly set the input region! <b>raster_output</b> is the same output, except converted to a raster grid at the same resolution as the current computational region. If you have a grid spacing that is much smaller than a flexural wavelength, it is possible to interpolate the vector output to a much finer resolution than this raster output provides.
+<p>
+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>
+<a href="v.surf.bspline">v.surf.bspline</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.
+
+<h2>AUTHORS</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/vector/v.flexure/v.flexure.py
===================================================================
--- grass-addons/grass7/vector/v.flexure/v.flexure.py (rev 0)
+++ grass-addons/grass7/vector/v.flexure/v.flexure.py 2015-02-03 23:58:16 UTC (rev 64449)
@@ -0,0 +1,268 @@
+#! /usr/bin/python
+############################################################################
+#
+# MODULE: v.flexure
+#
+# AUTHOR(S): Andrew Wickert
+#
+# PURPOSE: Calculate flexure of the lithosphere under a specified
+# set of loads and with a given elastic thickness (scalar)
+#
+# COPYRIGHT: (c) 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 20 Jan 2015 to add GRASS GIS support for distributed point loads
+# and their effects on lithospheric flexure
+
+#%module
+#% description: Lithospheric flexure: gridded deflections from scattered point loads
+#% keyword: vector
+#% keyword: geophysics
+#%end
+#%option
+#% key: q0
+#% type: string
+#% gisprompt: old,vector,vector
+#% description: Vector map of loads (thickness * area * density * g) [N]
+#% required : yes
+#%end
+#%option
+#% key: column
+#% type: string
+#% description: Column containing load values [N]
+#% required : yes
+#%end
+#%option
+#% key: te
+#% type: double
+#% description: Elastic thicnkess: scalar; 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,vector,vector
+#% description: Output vector points map of vertical deflections [m]
+#% required : yes
+#%end
+#%option
+#% key: raster_output
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: Output raster map of vertical deflections [m]
+#% 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
+from grass.pygrass import vector
+# 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.")
+
+
+####################
+# UTILITY FUNCTION #
+####################
+
+def get_points_xy(vect_name):
+ """
+ to find x and y using pygrass, see my (A. Wickert's) StackOverflow answer:
+ http://gis.stackexchange.com/questions/28061/how-to-access-vector-coordinates-in-grass-gis-from-python
+ """
+ points = vector.VectorTopo(vect_name)
+ points.open('r')
+ coords = []
+ for i in range(len(points)):
+ coords.append(points.read(i+1).coords())
+ coords = np.array(coords)
+ return coords[:,0], coords[:,1] # x, y
+
+
+############################
+# PASS VARIABLES AND SOLVE #
+############################
+
+def main():
+ """
+ Superposition of analytical solutions in gFlex for flexural isostasy in
+ GRASS GIS
+ """
+
+ ##########
+ # SET-UP #
+ ##########
+
+ # This code is for 2D flexural isostasy
+ flex = gflex.F2D()
+ # And show that it is coming from GRASS GIS
+ flex.grass = True
+
+ # Method
+ flex.Method = 'SAS_NG'
+
+ # Parameters that are often changed for the solution
+ ######################################################
+
+ # x, y, q
+ flex.x, flex.y = get_points_xy(options['q0'])
+ # xw, yw: gridded output
+ if len(grass.parse_command('g.list', type='vect', pattern=options['output'])):
+ if not grass.overwrite():
+ grass.fatal("Vector map '" + options['output'] + "' already exists. Use '--o' to overwrite.")
+ # Just check raster at the same time if it exists
+ if len(grass.parse_command('g.list', type='rast', pattern=options['raster_output'])):
+ if not grass.overwrite():
+ grass.fatal("Raster map '" + options['raster_output'] + "' already exists. Use '--o' to overwrite.")
+ grass.run_command('v.mkgrid', map=options['output'], type='point', overwrite=grass.overwrite(), quiet=True)
+ grass.run_command('v.db.addcolumn', map=options['output'], columns='w double precision', quiet=True)
+ flex.xw, flex.yw = get_points_xy(options['output']) # gridded output coordinates
+ vect_db = grass.vector_db_select(options['q0'])
+ col_names = np.array(vect_db['columns'])
+ q_col = (col_names == options['column'])
+ if np.sum(q_col):
+ col_values = np.array(vect_db['values'].values()).astype(float)
+ flex.q = col_values[:, q_col].squeeze() # Make it 1D for consistency w/ x, y
+ else:
+ grass.fatal("provided column name, "+options['column']+" does not match\nany column in "+options['q0']+".")
+ # Elastic thickness
+ flex.Te = float(options['te'])
+ if options['te_units'] == 'km':
+ flex.Te *= 1000
+ elif options['te_units'] == 'm':
+ pass
+ else:
+ grass.fatal("Inappropriate te_units. How? Options should be limited by GRASS.")
+ 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'])
+
+ # Set verbosity
+ if grass.verbosity() >= 2:
+ flex.Verbose = True
+ if grass.verbosity() >= 3:
+ flex.Debug = True
+ elif grass.verbosity() == 0:
+ flex.Quiet = True
+
+ # Check if lat/lon and let user know if verbosity is True
+ if grass.region_env()[6] == '3':
+ flex.latlon = True
+ flex.PlanetaryRadius = float(grass.parse_command('g.proj', flags='j')['+a'])
+ if flex.Verbose:
+ print "Latitude/longitude grid."
+ print "Based on r_Earth = 6371 km"
+ print "Computing distances between load points using great circle paths"
+
+ ##########
+ # SOLVE! #
+ ##########
+
+ flex.initialize()
+ flex.run()
+ flex.finalize()
+
+ # Now to use lower-level GRASS vector commands to work with the database
+ # table and update its entries
+ # See for help:
+ # http://nbviewer.ipython.org/github/zarch/workshop-pygrass/blob/master/02_Vector.ipynb
+ w = vector.VectorTopo(options['output'])
+ w.open('rw') # Get ready to read and write
+ wdb = w.dblinks[0]
+ wtable = wdb.table()
+ col = int((np.array(wtable.columns.names()) == 'w').nonzero()[0]) # update this column
+ for i in range(1, len(w)+1):
+ # ignoring 1st column: assuming it will be category (always true here)
+ wnewvalues = w[i].attrs.values()[1:col] + tuple([flex.w[i-1]]) + w[i].attrs.values()[col+1:]
+ wtable.update(key=i, values=wnewvalues)
+ wtable.conn.commit() # Save this
+ w.close(build=False) # don't build here b/c it is always verbose
+ grass.run_command('v.build', map=options['output'], quiet=True)
+
+ # And raster export
+ # "w" vector defined by raster resolution, so can do direct v.to.rast
+ # though if this option isn't selected, the user can do a finer-grained
+ # interpolation, which shouldn't introduce much error so long as these
+ # outputs are spaced at << 1 flexural wavelength.
+ if options['raster_output']:
+ grass.run_command('v.to.rast', input=options['output'], output=options['raster_output'], use='attr', attribute_column='w', type='point', overwrite=grass.overwrite(), quiet=True)
+ # And create a nice colormap!
+ grass.run_command('r.colors', map=options['raster_output'], color='differences', quiet=True)
+
+if __name__ == "__main__":
+ options, flags = grass.parser()
+ main()
+
More information about the grass-commit
mailing list