[GRASS-SVN] r71511 - in grass-addons/grass7/vector: . v.gsflow.export v.gsflow.gravres v.gsflow.grid v.gsflow.hruparams v.gsflow.reaches v.gsflow.segments v.stream.inbasin v.stream.network
svn_grass at osgeo.org
svn_grass at osgeo.org
Wed Sep 27 23:08:49 PDT 2017
Author: awickert
Date: 2017-09-27 23:08:49 -0700 (Wed, 27 Sep 2017)
New Revision: 71511
Added:
grass-addons/grass7/vector/v.gsflow.export/
grass-addons/grass7/vector/v.gsflow.export/Makefile
grass-addons/grass7/vector/v.gsflow.export/v.gsflow.export.html
grass-addons/grass7/vector/v.gsflow.export/v.gsflow.export.py
grass-addons/grass7/vector/v.gsflow.gravres/
grass-addons/grass7/vector/v.gsflow.gravres/Makefile
grass-addons/grass7/vector/v.gsflow.gravres/v.gsflow.gravres.html
grass-addons/grass7/vector/v.gsflow.gravres/v.gsflow.gravres.py
grass-addons/grass7/vector/v.gsflow.grid/
grass-addons/grass7/vector/v.gsflow.grid/Makefile
grass-addons/grass7/vector/v.gsflow.grid/v.gsflow.grid.html
grass-addons/grass7/vector/v.gsflow.grid/v.gsflow.grid.py
grass-addons/grass7/vector/v.gsflow.hruparams/
grass-addons/grass7/vector/v.gsflow.hruparams/Makefile
grass-addons/grass7/vector/v.gsflow.hruparams/v.gsflow.hruparams.html
grass-addons/grass7/vector/v.gsflow.hruparams/v.gsflow.hruparams.py
grass-addons/grass7/vector/v.gsflow.reaches/
grass-addons/grass7/vector/v.gsflow.reaches/Makefile
grass-addons/grass7/vector/v.gsflow.reaches/v.gsflow.reaches.html
grass-addons/grass7/vector/v.gsflow.reaches/v.gsflow.reaches.py
grass-addons/grass7/vector/v.gsflow.segments/
grass-addons/grass7/vector/v.gsflow.segments/Makefile
grass-addons/grass7/vector/v.gsflow.segments/v.gsflow.segments.html
grass-addons/grass7/vector/v.gsflow.segments/v.gsflow.segments.py
grass-addons/grass7/vector/v.stream.inbasin/
grass-addons/grass7/vector/v.stream.inbasin/Makefile
grass-addons/grass7/vector/v.stream.inbasin/v.stream.inbasin.html
grass-addons/grass7/vector/v.stream.inbasin/v.stream.inbasin.py
grass-addons/grass7/vector/v.stream.network/
grass-addons/grass7/vector/v.stream.network/Makefile
grass-addons/grass7/vector/v.stream.network/v.stream.network.html
grass-addons/grass7/vector/v.stream.network/v.stream.network.py
Log:
Stream networks and GSFLOW
Added: grass-addons/grass7/vector/v.gsflow.export/Makefile
===================================================================
--- grass-addons/grass7/vector/v.gsflow.export/Makefile (rev 0)
+++ grass-addons/grass7/vector/v.gsflow.export/Makefile 2017-09-28 06:08:49 UTC (rev 71511)
@@ -0,0 +1,8 @@
+MODULE_TOPDIR = ../..
+
+PGM = v.gsflow.export
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script
+
Added: grass-addons/grass7/vector/v.gsflow.export/v.gsflow.export.html
===================================================================
--- grass-addons/grass7/vector/v.gsflow.export/v.gsflow.export.html (rev 0)
+++ grass-addons/grass7/vector/v.gsflow.export/v.gsflow.export.html 2017-09-28 06:08:49 UTC (rev 71511)
@@ -0,0 +1,25 @@
+<h2>DESCRIPTION</h2>
+
+<em>v.gsflow.export</em> produces output tables and maps that are used for the USGS combined groundwater (MODFLOW) and surface-water (PRMS) model GSFLOW. These are read in through a sequence of <a href=https://github.com/UMN-Hydro/GSFLOW_pre-processor>"GSFLOW pre-processor" scripts available at the University of Minnesota's Hydro GitHub page</a>.
+
+<h2>REFERENCES</h2>
+
+Ng, G.-H. C., A. D. Wickert, R. L. McLaughlin, J. La Frenierre, S. Liess, and L. Sabeeri (2016),
+Modeling the role of groundwater and vegetation in the hydrological response of tropical
+glaciated watersheds to climate change, in AGU Fall Meeting Abstracts, H13L–1590, San
+Francisco, CA.
+
+<h2>SEE ALSO</h2>
+<a href="v.gsflow.gravres">v.gsflow.gravres</a>
+<a href="v.gsflow.grid">v.gsflow.grid</a>
+<a href="v.gsflow.hruparams">v.gsflow.hruparams</a>
+<a href="v.gsflow.reaches">v.gsflow.reaches</a>
+<a href="v.gsflow.segments">v.gsflow.segments</a>
+<a href="v.stream.inbasin">v.stream.inbasin</a>
+<a href="v.stream.network">v.stream.network</a>
+
+<h2>AUTHORS</h2>
+
+Andrew D. Wickert<br>
+
+<p><i>Last changed: $Date 2016-09-27$</i>
Added: grass-addons/grass7/vector/v.gsflow.export/v.gsflow.export.py
===================================================================
--- grass-addons/grass7/vector/v.gsflow.export/v.gsflow.export.py (rev 0)
+++ grass-addons/grass7/vector/v.gsflow.export/v.gsflow.export.py 2017-09-28 06:08:49 UTC (rev 71511)
@@ -0,0 +1,256 @@
+#!/usr/bin/env python
+############################################################################
+#
+# MODULE: v.gsflow.export
+#
+# AUTHOR(S): Andrew Wickert
+#
+# PURPOSE: Export database tables and pour point for GSFLOW input and control files
+#
+# COPYRIGHT: (c) 2016-2017 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:
+# - uses inputs from the v.gsflow series
+
+# More information
+# Started December 2016
+
+#%module
+#% description: Export databse tables and pour point for GSFLOW input and control files
+#% keyword: vector
+#% keyword: stream network
+#% keyword: hydrology
+#% keyword: GSFLOW
+#%end
+
+#%option G_OPT_V_INPUT
+#% key: reaches_input
+#% label: Reaches where stream segments overlap MODFLOW grid cells
+#% required: no
+#% guidependency: layer,column
+#%end
+
+#%option G_OPT_V_INPUT
+#% key: segments_input
+#% label: Stream segments, coincident with HRUs
+#% required: no
+#% guidependency: layer,column
+#%end
+
+#%option G_OPT_V_INPUT
+#% key: gravres_input
+#% label: Union of MODFLOW grid and HRUs
+#% required: no
+#% guidependency: layer,column
+#%end
+
+#%option G_OPT_V_INPUT
+#% key: hru_input
+#% label: PRMS-style hydrologic response units
+#% required: no
+#% guidependency: layer,column
+#%end
+
+#%option G_OPT_V_INPUT
+#% key: pour_point_input
+#% label: Pour point at the outlet of the basin
+#% required: no
+#% guidependency: layer,column
+#%end
+
+#%option G_OPT_V_OUTPUT
+#% key: reaches_output
+#% label: Reaches table, no file ext
+#% required: no
+#% guidependency: layer,column
+#%end
+
+#%option G_OPT_V_OUTPUT
+#% key: segments_output
+#% label: Segments table, no file ext
+#% required: no
+#% guidependency: layer,column
+#%end
+
+#%option G_OPT_V_OUTPUT
+#% key: gravres_output
+#% label: Gravity Reservoir output table for GSFLOW input, no file ext
+#% required: no
+#% guidependency: layer,column
+#%end
+
+#%option G_OPT_V_OUTPUT
+#% key: hru_output
+#% label: HRU table, no file ext
+#% required: no
+#% guidependency: layer,column
+#%end
+
+#%option G_OPT_V_OUTPUT
+#% key: pour_point_output
+#% label: Pour point coordinates for GSFLOW input, no file ext
+#% required: no
+#% guidependency: layer,column
+#%end
+
+##################
+# IMPORT MODULES #
+##################
+
+# PYTHON
+import numpy as np
+from matplotlib import pyplot as plt
+import sys
+# GRASS
+from grass.pygrass.modules.shortcuts import general as g
+from grass.pygrass.modules.shortcuts import raster as r
+from grass.pygrass.modules.shortcuts import vector as v
+from grass.pygrass.modules.shortcuts import miscellaneous as m
+from grass.pygrass.gis import region
+from grass.pygrass import vector # Change to "v"?
+from grass.script import vector_db_select
+from grass.pygrass.vector import Vector, VectorTopo
+from grass.pygrass.raster import RasterRow
+from grass.pygrass import utils
+from grass import script as gscript
+
+###############
+# MAIN MODULE #
+###############
+
+def get_columns_in_order(vect, cols, nodata_value=-999):
+ colNames = np.array(gscript.vector_db_select(vect, layer=1)['columns'])
+ colValues = np.array(gscript.vector_db_select(vect, layer=1)\
+ ['values'].values())
+ outlist = []
+ for col in cols:
+ newcol = colValues[:,colNames == col].squeeze()
+ # If column does not exist, populate with nodata value
+ # Strangely, has shape (nrows, 0)
+ if np.prod(newcol.shape) == 0:
+ newcol = (-999*np.ones(colValues.shape[0])).astype(int).astype(str)
+ outlist.append(newcol)
+ return outlist
+
+def main():
+ """
+ Build gravity reservoirs in GSFLOW: combines MODFLOW grid and HRU sub-basins
+ These define the PRMS soil zone that connects to MODFLOW cells
+ """
+
+ ##################
+ # OPTION PARSING #
+ ##################
+
+ # Parser
+ options, flags = gscript.parser()
+
+ # Input
+ reaches = options['reaches_input']
+ segments = options['segments_input']
+ gravity_reservoirs = options['gravres_input']
+ HRUs = options['hru_input']
+ pour_point = options['pour_point_input']
+
+ # Output
+ out_reaches = options['reaches_output']
+ out_segments = options['segments_output']
+ out_gravity_reservoirs = options['gravres_output']
+ out_HRUs = options['hru_output']
+ out_pour_point = options['pour_point_output']
+
+ ##############
+ # PROCESSING #
+ ##############
+
+ # Reaches
+ ##########
+ if (len(reaches) > 0) and (len(out_reaches) > 0):
+ columns_in_order = ['KRCH', 'IRCH', 'JRCH', 'ISEG', 'IREACH', 'RCHLEN',
+ 'STRTOP', 'SLOPE', 'STRTHICK', 'STRHC1', 'THTS',
+ 'THTI', 'EPS', 'UHC']
+ outcols = get_columns_in_order(reaches, columns_in_order)
+ outarray = np.array(outcols).transpose()
+ outtable = np.vstack((columns_in_order, outarray))
+ np.savetxt(out_reaches+'.txt', outtable, fmt='%s', delimiter=',')
+ elif (len(reaches) > 0) or (len(out_reaches) > 0):
+ grass.fatal(_("You must inlcude both input and output reaches"))
+
+ # Segments
+ ###########
+ if (len(segments) > 0) and (len(out_segments) > 0):
+ columns_in_order = ['NSEG', 'ICALC', 'OUTSEG', 'IUPSEG', 'IPRIOR',
+ 'NSTRPTS', 'FLOW', 'RUNOFF', 'ETSW', 'PPTSW',
+ 'ROUGHCH', 'ROUGHBK', 'CDPTH', 'FDPTH', 'AWDTH',
+ 'BWDTH']
+ outcols = get_columns_in_order(segments, columns_in_order)
+ outarray = np.array(outcols).transpose()
+ outtable = np.vstack((columns_in_order, outarray))
+ np.savetxt(out_segments+'_4A_INFORMATION.txt', outtable, fmt='%s',
+ delimiter=',')
+
+ columns_in_order = ['HCOND1', 'THICKM1', 'ELEVUP', 'WIDTH1', 'DEPTH1',
+ 'THTS1', 'THTI1', 'EPS1', 'UHC1']
+ outcols = get_columns_in_order(segments, columns_in_order)
+ outarray = np.array(outcols).transpose()
+ outtable = np.vstack((columns_in_order, outarray))
+ np.savetxt(out_segments+'_4B_UPSTREAM.txt', outtable, fmt='%s',
+ delimiter=',')
+
+ columns_in_order = ['HCOND2', 'THICKM2', 'ELEVDN', 'WIDTH2', 'DEPTH2',
+ 'THTS2', 'THTI2', 'EPS2', 'UHC2']
+ outcols = get_columns_in_order(segments, columns_in_order)
+ outarray = np.array(outcols).transpose()
+ outtable = np.vstack((columns_in_order, outarray))
+ np.savetxt(out_segments+'_4C_DOWNSTREAM.txt', outtable, fmt='%s',
+ delimiter=',')
+ elif (len(segments) > 0) or (len(out_segments) > 0):
+ grass.fatal(_("You must inlcude both input and output segments"))
+
+ # Gravity reservoirs
+ #####################
+ if (len(gravity_reservoirs) > 0) and (len(out_gravity_reservoirs) > 0):
+ columns_in_order = ['gvr_hru_id', 'gvr_hru_pct', 'gvr_cell_id',
+ 'gvr_cell_pct']
+ outcols = get_columns_in_order(gravity_reservoirs, columns_in_order)
+ outarray = np.array(outcols).transpose()
+ outtable = np.vstack((columns_in_order, outarray))
+ np.savetxt(out_gravity_reservoirs+'.txt', outtable, fmt='%s',
+ delimiter=',')
+ elif (len(gravity_reservoirs) > 0) or (len(out_gravity_reservoirs) > 0):
+ grass.fatal(_("You must inlcude both input and output \
+ gravity reservoirs"))
+
+ # HRUs
+ #######
+ if (len(HRUs) > 0) and (len(out_HRUs) > 0):
+ columns_in_order = ['hru_area', 'hru_aspect', 'hru_elev', 'hru_lat',
+ 'hru_slope', 'hru_segment', 'hru_strmseg_down_id']
+ outcols = get_columns_in_order(HRUs, columns_in_order)
+ outarray = np.array(outcols).transpose()
+ outtable = np.vstack((columns_in_order, outarray))
+ np.savetxt(HRUs+'.txt', outtable, fmt='%s', delimiter=',')
+ elif (len(HRUs) > 0) or (len(out_HRUs) > 0):
+ grass.fatal(_("You must inlcude both input and output HRUs"))
+
+ # Pour Point
+ #############
+ if (len(pour_point) > 0) and (len(out_pour_point) > 0):
+ _y, _x = np.squeeze(gscript.db_select(sql='SELECT row,col FROM '+
+ pour_point))
+ outstr = 'discharge_pt: row_i '+_y+' col_i '+_x
+ outfile = file(out_pour_point+'.txt', 'w')
+ outfile.write(outstr)
+ outfile.close()
+ elif (len(pour_point) > 0) or (len(out_pour_point) > 0):
+ grass.fatal(_("You must inlcude both input and output pour points"))
+
+
+if __name__ == "__main__":
+ main()
Property changes on: grass-addons/grass7/vector/v.gsflow.export/v.gsflow.export.py
___________________________________________________________________
Added: svn:mime-type
+ text/x-python
Added: svn:eol-style
+ native
Added: grass-addons/grass7/vector/v.gsflow.gravres/Makefile
===================================================================
--- grass-addons/grass7/vector/v.gsflow.gravres/Makefile (rev 0)
+++ grass-addons/grass7/vector/v.gsflow.gravres/Makefile 2017-09-28 06:08:49 UTC (rev 71511)
@@ -0,0 +1,8 @@
+MODULE_TOPDIR = ../..
+
+PGM = v.gsflow.gravres
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script
+
Added: grass-addons/grass7/vector/v.gsflow.gravres/v.gsflow.gravres.html
===================================================================
--- grass-addons/grass7/vector/v.gsflow.gravres/v.gsflow.gravres.html (rev 0)
+++ grass-addons/grass7/vector/v.gsflow.gravres/v.gsflow.gravres.html 2017-09-28 06:08:49 UTC (rev 71511)
@@ -0,0 +1,25 @@
+<h2>DESCRIPTION</h2>
+
+<em>v.gsflow.gravres</em> produces "gravity reservoirs" that are used for the USGS combined groundwater (MODFLOW) and surface-water (PRMS) model GSFLOW. Each "gravity reservoir" is the union between a hydrologic response unit (HRU) and part or all of a co-located MODFLOW grid cell. These link surface-water and groundwater flow.
+
+<h2>REFERENCES</h2>
+
+Ng, G.-H. C., A. D. Wickert, R. L. McLaughlin, J. La Frenierre, S. Liess, and L. Sabeeri (2016),
+Modeling the role of groundwater and vegetation in the hydrological response of tropical
+glaciated watersheds to climate change, in AGU Fall Meeting Abstracts, H13L–1590, San
+Francisco, CA.
+
+<h2>SEE ALSO</h2>
+<a href="v.gsflow.export">v.gsflow.export</a>
+<a href="v.gsflow.grid">v.gsflow.grid</a>
+<a href="v.gsflow.hruparams">v.gsflow.hruparams</a>
+<a href="v.gsflow.reaches">v.gsflow.reaches</a>
+<a href="v.gsflow.segments">v.gsflow.segments</a>
+<a href="v.stream.inbasin">v.stream.inbasin</a>
+<a href="v.stream.network">v.stream.network</a>
+
+<h2>AUTHORS</h2>
+
+Andrew D. Wickert<br>
+
+<p><i>Last changed: $Date 2016-09-27$</i>
Added: grass-addons/grass7/vector/v.gsflow.gravres/v.gsflow.gravres.py
===================================================================
--- grass-addons/grass7/vector/v.gsflow.gravres/v.gsflow.gravres.py (rev 0)
+++ grass-addons/grass7/vector/v.gsflow.gravres/v.gsflow.gravres.py 2017-09-28 06:08:49 UTC (rev 71511)
@@ -0,0 +1,126 @@
+#!/usr/bin/env python
+############################################################################
+#
+# MODULE: v.gsflow.gravres
+#
+# AUTHOR(S): Andrew Wickert
+#
+# PURPOSE: Build gravity reservoirs -- intersections of HRU sub-basins and
+# MODFLOW grid cells -- for GSFLOW
+#
+# COPYRIGHT: (c) 2016-2017 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:
+# - uses inputs from r.stream.basins & r.to.vect
+# - uses inputs from v.gsflow.grid
+
+# More information
+# Started December 2016
+
+#%module
+#% description: Set parameters for GSFLOW Hydrologic Response Units (HRUs)
+#% keyword: vector
+#% keyword: stream network
+#% keyword: hydrology
+#% keyword: GSFLOW
+#%end
+
+#%option G_OPT_V_INPUT
+#% key: hru_input
+#% label: Sub-basins
+#% required: yes
+#% guidependency: layer,column
+#%end
+
+#%option G_OPT_V_INPUT
+#% key: grid_input
+#% label: MODFLOW grid
+#% required: yes
+#% guidependency: layer,column
+#%end
+
+#%option G_OPT_V_OUTPUT
+#% key: output
+#% label: Gravity reservoirs: union (AND) of sub-basins and MODFLOW grid
+#% required: yes
+#% guidependency: layer,column
+#%end
+
+##################
+# IMPORT MODULES #
+##################
+# PYTHON
+import numpy as np
+from matplotlib import pyplot as plt
+import sys
+# GRASS
+from grass.pygrass.modules.shortcuts import general as g
+from grass.pygrass.modules.shortcuts import raster as r
+from grass.pygrass.modules.shortcuts import vector as v
+from grass.pygrass.modules.shortcuts import miscellaneous as m
+from grass.pygrass.gis import region
+from grass.pygrass import vector # Change to "v"?
+from grass.script import vector_db_select
+from grass.pygrass.vector import Vector, VectorTopo
+from grass.pygrass.raster import RasterRow
+from grass.pygrass import utils
+from grass import script as gscript
+
+###############
+# MAIN MODULE #
+###############
+
+def main():
+ """
+ Build gravity reservoirs in GSFLOW: combines MODFLOW grid and HRU sub-basins
+ These define the PRMS soil zone that connects to MODFLOW cells
+ """
+
+ ##################
+ # OPTION PARSING #
+ ##################
+
+ # I/O
+ options, flags = gscript.parser()
+
+ # I/O
+ HRUs = options['hru_input']
+ grid = options['grid_input']
+ segments = options['output']
+ #col = options['col']
+ gravity_reservoirs = options['output']
+
+ ############
+ # ANALYSIS #
+ ############
+
+ """
+ # Basin areas
+ v.db_addcolumn(map=basins, columns=col)
+ v.to_db(map=basins, option='area', units='meters', columns=col)
+ """
+
+ # Create gravity reservoirs -- overlay cells=grid and HRUs
+ v.overlay(ainput=HRUs, binput=grid, operator='and', output=gravity_reservoirs, overwrite=gscript.overwrite())
+ v.db_dropcolumn(map=gravity_reservoirs, columns='a_cat,a_label,b_cat', quiet=True)
+ # Cell and HRU ID's
+ v.db_renamecolumn(map=gravity_reservoirs, column=('a_id', 'gvr_hru_id'), quiet=True)
+ v.db_renamecolumn(map=gravity_reservoirs, column=('b_id', 'gvr_cell_id'), quiet=True)
+ # Percent areas
+ v.db_renamecolumn(map=gravity_reservoirs, column=('a_hru_area_m2', 'hru_area_m2'), quiet=True)
+ v.db_renamecolumn(map=gravity_reservoirs, column=('b_area_m2', 'cell_area_m2'), quiet=True)
+ v.db_addcolumn(map=gravity_reservoirs, columns='area_m2', quiet=True)
+ v.to_db(map=gravity_reservoirs, option='area', units='meters', columns='area_m2', quiet=True)
+ v.db_addcolumn(map=gravity_reservoirs, columns='gvr_cell_pct double precision, gvr_hru_pct double precision', quiet=True)
+ v.db_update(map=gravity_reservoirs, column='gvr_cell_pct', query_column='100*area_m2/cell_area_m2', quiet=True)
+ v.db_update(map=gravity_reservoirs, column='gvr_hru_pct', query_column='100*area_m2/hru_area_m2', quiet=True)
+
+
+if __name__ == "__main__":
+ main()
Added: grass-addons/grass7/vector/v.gsflow.grid/Makefile
===================================================================
--- grass-addons/grass7/vector/v.gsflow.grid/Makefile (rev 0)
+++ grass-addons/grass7/vector/v.gsflow.grid/Makefile 2017-09-28 06:08:49 UTC (rev 71511)
@@ -0,0 +1,8 @@
+MODULE_TOPDIR = ../..
+
+PGM = v.gsflow.grid
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script
+
Added: grass-addons/grass7/vector/v.gsflow.grid/v.gsflow.grid.html
===================================================================
--- grass-addons/grass7/vector/v.gsflow.grid/v.gsflow.grid.html (rev 0)
+++ grass-addons/grass7/vector/v.gsflow.grid/v.gsflow.grid.html 2017-09-28 06:08:49 UTC (rev 71511)
@@ -0,0 +1,23 @@
+<h2>DESCRIPTION</h2>
+
+<em>v.gsflow.grid</em> generates the MODFLOW grid and associated supporting basin mask and topography files that are used for the USGS combined groundwater (MODFLOW) and surface-water (PRMS) model GSFLOW.
+
+<h2>REFERENCES</h2>
+
+Ng, G.-H. C., A. D. Wickert, R. L. McLaughlin, J. La Frenierre, S. Liess, and L. Sabeeri (2016),
+Modeling the role of groundwater and vegetation in the hydrological response of tropical
+glaciated watersheds to climate change, in AGU Fall Meeting Abstracts, H13L–1590, San
+Francisco, CA.
+
+<h2>AUTHORS</h2>
+
+Andrew D. Wickert<br>LSO</h2>
+<a href="v.gsflow.export">v.gsflow.export</a>
+<a href="v.gsflow.gravres">v.gsflow.gravres</a>
+<a href="v.gsflow.hruparams">v.gsflow.hruparams</a>
+<a href="v.gsflow.reaches">v.gsflow.reaches</a>
+<a href="v.gsflow.segments">v.gsflow.segments</a>
+<a href="v.stream.inbasin">v.stream.inbasin</a>
+<a href="v.stream.network">v.stream.network</a>
+
+<p><i>Last changed: $Date 2016-09-27$</i>
Added: grass-addons/grass7/vector/v.gsflow.grid/v.gsflow.grid.py
===================================================================
--- grass-addons/grass7/vector/v.gsflow.grid/v.gsflow.grid.py (rev 0)
+++ grass-addons/grass7/vector/v.gsflow.grid/v.gsflow.grid.py 2017-09-28 06:08:49 UTC (rev 71511)
@@ -0,0 +1,150 @@
+#!/usr/bin/env python
+############################################################################
+#
+# MODULE: v.gsflow.reaches
+#
+# AUTHOR(S): Andrew Wickert
+#
+# PURPOSE: Builds grid for the MODFLOW component of GSFLOW
+#
+# COPYRIGHT: (c) 2016-2017 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:
+# - uses inputs from r.stream.extract
+
+# More information
+# Started December 2016
+
+#%module
+#% description: Builds grid for the MODFLOW component of GSFLOW
+#% keyword: vector
+#% keyword: stream network
+#% keyword: hydrology
+#% keyword: GSFLOW
+#%end
+
+#%option G_OPT_V_INPUT
+#% key: basin
+#% label: Study basin, over which to build a MODFLOW grid
+#% required: yes
+#%end
+
+#%option G_OPT_V_INPUT
+#% key: pour_point
+#% label: Pour point, to which row and col (MODFLOW) will be added
+#% required: yes
+#%end
+
+#%option
+#% key: dx
+#% label: Cell size (x / E / zonal), in map units
+#% required: yes
+#%end
+
+#%option
+#% key: dy
+#% label: Cell size (y / N / meridional), in map units
+#% required: yes
+#%end
+
+#%option G_OPT_V_OUTPUT
+#% key: output
+#% label: MODFLOW grid
+#% required: yes
+#%end
+
+#%option G_OPT_R_OUTPUT
+#% key: mask_output
+#% label: Basin mask: inside (1) or outside (0) the watershed?
+#% required: no
+#%end
+
+##################
+# IMPORT MODULES #
+##################
+# PYTHON
+import numpy as np
+from matplotlib import pyplot as plt
+import sys
+import warnings
+# GRASS
+from grass.pygrass.modules.shortcuts import general as g
+from grass.pygrass.modules.shortcuts import raster as r
+from grass.pygrass.modules.shortcuts import vector as v
+from grass.pygrass.modules.shortcuts import miscellaneous as m
+from grass.pygrass.gis import region
+from grass.pygrass import vector # Change to "v"?
+from grass.script import vector_db_select
+from grass.pygrass.vector import Vector, VectorTopo
+from grass.pygrass.raster import RasterRow
+from grass.pygrass import utils
+from grass import script as gscript
+
+###############
+# MAIN MODULE #
+###############
+
+def main():
+ """
+ Builds a grid for the MODFLOW component of the USGS hydrologic model,
+ GSFLOW.
+ """
+
+ options, flags = gscript.parser()
+ basin = options['basin']
+ dx = options['dx']
+ dy = options['dy']
+ grid = options['output']
+ mask = options['mask_output']
+ pp = options['pour_point']
+
+ # Create grid
+ gscript.use_temp_region()
+ g.region(vector=basin, ewres=dx, nsres=dy)
+ v.mkgrid(map=grid, overwrite=gscript.overwrite())
+
+ # Cell numbers (row, column, continuous ID)
+ v.db_addcolumn(map=grid, columns='id int', quiet=True)
+ colNames = np.array(gscript.vector_db_select(grid, layer=1)['columns'])
+ colValues = np.array(gscript.vector_db_select(grid, layer=1)['values'].values())
+ cats = colValues[:,colNames == 'cat'].astype(int).squeeze()
+ rows = colValues[:,colNames == 'row'].astype(int).squeeze()
+ cols = colValues[:,colNames == 'col'].astype(int).squeeze()
+ nrows = np.max(rows)
+ ncols = np.max(cols)
+ cats = np.ravel([cats])
+ _id = np.ravel([ncols * (rows - 1) + cols])
+ _id_cat = []
+ for i in range(len(_id)):
+ _id_cat.append( (_id[i], cats[i]) )
+ gridTopo = VectorTopo(grid)
+ gridTopo.open('rw')
+ cur = gridTopo.table.conn.cursor()
+ cur.executemany("update "+grid+" set id=? where cat=?", _id_cat)
+ gridTopo.table.conn.commit()
+ gridTopo.close()
+
+ # Cell area
+ v.db_addcolumn(map=grid, columns='area_m2', quiet=True)
+ v.to_db(map=grid, option='area', units='meters', columns='area_m2', quiet=True)
+
+ # Basin mask
+ if len(mask) > 0:
+ v.to_rast(input=basin, output=mask, use='val', value=1, quiet=True)
+
+ # Pour point
+ if len(pp) > 0:
+ v.db_addcolumn(map=pp, columns=('row integer','col integer'), quiet=True)
+ v.build(map=pp, quiet=True)
+ v.what_vect(map=pp, query_map=grid, column='row', query_column='row', quiet=True)
+ v.what_vect(map=pp, query_map=grid, column='col', query_column='col', quiet=True)
+
+
+if __name__ == "__main__":
+ main()
Added: grass-addons/grass7/vector/v.gsflow.hruparams/Makefile
===================================================================
--- grass-addons/grass7/vector/v.gsflow.hruparams/Makefile (rev 0)
+++ grass-addons/grass7/vector/v.gsflow.hruparams/Makefile 2017-09-28 06:08:49 UTC (rev 71511)
@@ -0,0 +1,8 @@
+MODULE_TOPDIR = ../..
+
+PGM = v.gsflow.hruparams
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script
+
Added: grass-addons/grass7/vector/v.gsflow.hruparams/v.gsflow.hruparams.html
===================================================================
--- grass-addons/grass7/vector/v.gsflow.hruparams/v.gsflow.hruparams.html (rev 0)
+++ grass-addons/grass7/vector/v.gsflow.hruparams/v.gsflow.hruparams.html 2017-09-28 06:08:49 UTC (rev 71511)
@@ -0,0 +1,25 @@
+<h2>DESCRIPTION</h2>
+
+<em>v.gsflow.hruparameters</em> uploads required parameters to run the USGS combined groundwater (MODFLOW) and surface-water (PRMS) model GSFLOW to each of the hydrologic response units (HRU) in the study watershed.
+
+<h2>REFERENCES</h2>
+
+Ng, G.-H. C., A. D. Wickert, R. L. McLaughlin, J. La Frenierre, S. Liess, and L. Sabeeri (2016),
+Modeling the role of groundwater and vegetation in the hydrological response of tropical
+glaciated watersheds to climate change, in AGU Fall Meeting Abstracts, H13L–1590, San
+Francisco, CA.
+
+<h2>SEE ALSO</h2>
+<a href="v.gsflow.export">v.gsflow.export</a>
+<a href="v.gsflow.gravres">v.gsflow.gravres</a>
+<a href="v.gsflow.grid">v.gsflow.grid</a>
+<a href="v.gsflow.reaches">v.gsflow.reaches</a>
+<a href="v.gsflow.segments">v.gsflow.segments</a>
+<a href="v.stream.inbasin">v.stream.inbasin</a>
+<a href="v.stream.network">v.stream.network</a>
+
+<h2>AUTHORS</h2>
+
+Andrew D. Wickert<br>
+
+<p><i>Last changed: $Date 2016-09-27$</i>
Added: grass-addons/grass7/vector/v.gsflow.hruparams/v.gsflow.hruparams.py
===================================================================
--- grass-addons/grass7/vector/v.gsflow.hruparams/v.gsflow.hruparams.py (rev 0)
+++ grass-addons/grass7/vector/v.gsflow.hruparams/v.gsflow.hruparams.py 2017-09-28 06:08:49 UTC (rev 71511)
@@ -0,0 +1,313 @@
+#!/usr/bin/env python
+############################################################################
+#
+# MODULE: v.gsflow.hruparams
+#
+# AUTHOR(S): Andrew Wickert
+#
+# PURPOSE: Set parameters for GSFLOW Hydrologic Response Units (HRUs)
+#
+# COPYRIGHT: (c) 2016-2017 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:
+# - uses inputs from v.stream.network
+# - uses inputs from r.slope.aspect
+
+# More information
+# Started December 2016
+
+#%module
+#% description: Set parameters for GSFLOW Hydrologic Response Units (HRUs)
+#% keyword: vector
+#% keyword: stream network
+#% keyword: hydrology
+#% keyword: GSFLOW
+#%end
+
+#%option G_OPT_R_INPUT
+#% key: elevation
+#% label: Elevation raster
+#% required: yes
+#% guidependency: layer,column
+#%end
+
+#%option G_OPT_V_INPUT
+#% key: input
+#% label: Sub-basins to become HRUs
+#% required: yes
+#% guidependency: layer,column
+#%end
+
+#%option G_OPT_V_OUTPUT
+#% key: output
+#% label: HRUs
+#% required: yes
+#% guidependency: layer,column
+#%end
+
+#%option
+#% key: slope
+#% type: string
+#% description: Slope [unitless]: r.slope.aspect format=percent zscale=0.01
+#% required: yes
+#%end
+
+#%option
+#% key: aspect
+#% type: string
+#% description: Aspect from r.slope.aspect
+#% required: yes
+#%end
+
+##################
+# IMPORT MODULES #
+##################
+# PYTHON
+import numpy as np
+from matplotlib import pyplot as plt
+import sys
+# GRASS
+from grass.pygrass.modules.shortcuts import general as g
+from grass.pygrass.modules.shortcuts import raster as r
+from grass.pygrass.modules.shortcuts import vector as v
+from grass.pygrass.modules.shortcuts import miscellaneous as m
+from grass.pygrass.gis import region
+from grass.pygrass import vector # Change to "v"?
+from grass.script import vector_db_select
+from grass.pygrass.vector import Vector, VectorTopo
+from grass.pygrass.raster import RasterRow
+from grass.pygrass import utils
+from grass import script as gscript
+
+###############
+# MAIN MODULE #
+###############
+
+def main():
+ """
+ Adds GSFLOW parameters to a set of HRU sub-basins
+ """
+
+ ##################
+ # OPTION PARSING #
+ ##################
+
+ options, flags = gscript.parser()
+ basins = options['input']
+ HRU = options['output']
+ slope = options['slope']
+ aspect = options['aspect']
+ elevation = options['elevation']
+
+ ################################
+ # CREATE HRUs FROM SUB-BASINS #
+ ################################
+
+ g.copy(vector=(basins,HRU), overwrite=gscript.overwrite())
+
+ ############################################
+ # ATTRIBUTE COLUMNS (IN ORDER FROM MANUAL) #
+ ############################################
+
+ # HRU
+ hru_columns = []
+ # Self ID
+ hru_columns.append('id integer') # nhru
+ # Basic Physical Attributes (Geometry)
+ hru_columns.append('hru_area double precision') # acres (!!!!)
+ hru_columns.append('hru_area_m2 double precision') # [not for GSFLOW: for me!]
+ hru_columns.append('hru_aspect double precision') # Mean aspect [degrees]
+ hru_columns.append('hru_elev double precision') # Mean elevation
+ hru_columns.append('hru_lat double precision') # Latitude of centroid
+ hru_columns.append('hru_slope double precision') # Mean slope [percent]
+ # Basic Physical Attributes (Other)
+ #hru_columns.append('hru_type integer') # 0=inactive; 1=land; 2=lake; 3=swale; almost all will be 1
+ #hru_columns.append('elev_units integer') # 0=feet; 1=meters. 0=default. I think I will set this to 1 by default.
+ # Measured input
+ hru_columns.append('outlet_sta integer') # Index of streamflow station at basin outlet:
+ # station number if it has one, 0 if not
+ # Note that the below specify projections and note lat/lon; they really seem
+ # to work for any projected coordinates, with _x, _y, in meters, and _xlong,
+ # _ylat, in feet (i.e. they are just northing and easting). The meters and feet
+ # are not just simple conversions, but actually are required for different
+ # modules in the code, and are hence redundant but intentional.
+ hru_columns.append('hru_x double precision') # Easting [m]
+ hru_columns.append('hru_xlong double precision') # Easting [feet]
+ hru_columns.append('hru_y double precision') # Northing [m]
+ hru_columns.append('hru_ylat double precision') # Northing [feet]
+ # Streamflow and lake routing
+ hru_columns.append('K_coef double precision') # Travel time of flood wave to next downstream segment;
+ # this is the Muskingum storage coefficient
+ # 1.0 for reservoirs, diversions, and segments flowing
+ # out of the basin
+ hru_columns.append('x_coef double precision') # Amount of attenuation of flow wave;
+ # this is the Muskingum routing weighting factor
+ # range: 0.0--0.5; default 0.2
+ # 0 for all segments flowing out of the basin
+ hru_columns.append('hru_segment integer') # ID of stream segment to which flow will be routed
+ # this is for non-cascade routing (flow goes directly
+ # from HRU to stream segment)
+ hru_columns.append('obsin_segment integer') # Index of measured streamflow station that replaces
+ # inflow to a segment
+
+ # Create strings
+ hru_columns = ",".join(hru_columns)
+
+ # Add columns to tables
+ v.db_addcolumn(map=HRU, columns=hru_columns, quiet=True)
+
+
+ ###########################
+ # UPDATE DATABASE ENTRIES #
+ ###########################
+
+ colNames = np.array(gscript.vector_db_select(HRU, layer=1)['columns'])
+ colValues = np.array(gscript.vector_db_select(HRU, layer=1)['values'].values())
+ number_of_hrus = colValues.shape[0]
+ cats = colValues[:,colNames == 'cat'].astype(int).squeeze()
+ rnums = colValues[:,colNames == 'rnum'].astype(int).squeeze()
+
+ nhru = np.arange(1, number_of_hrus + 1)
+ nhrut = []
+ for i in range(len(nhru)):
+ nhrut.append( (nhru[i], cats[i]) )
+ # Access the HRUs
+ hru = VectorTopo(HRU)
+ # Open the map with topology:
+ hru.open('rw')
+ # Create a cursor
+ cur = hru.table.conn.cursor()
+ # Use it to loop across the table
+ cur.executemany("update "+HRU+" set id=? where cat=?", nhrut)
+ # Commit changes to the table
+ hru.table.conn.commit()
+ # Close the table
+ hru.close()
+
+ """
+ # Do the same for basins <-------------- DO THIS OR SIMPLY HAVE HRUs OVERLAIN WITH GRID CELLS? IN THIS CASE, RMV AREA ADDITION TO GRAVRES
+ v.db_addcolumn(map=basins, columns='id int', quiet=True)
+ basins = VectorTopo(basins)
+ basins.open('rw')
+ cur = basins.table.conn.cursor()
+ cur.executemany("update basins set id=? where cat=?", nhrut)
+ basins.table.conn.commit()
+ basins.close()
+ """
+
+ # if you want to append to table
+ # cur.executemany("update HRU(id) values(?)", nhrut) # "insert into" will add rows
+
+ #hru_columns.append('hru_area double precision')
+ # Acres b/c USGS
+ v.to_db(map=HRU, option='area', columns='hru_area', units='acres', quiet=True)
+ v.to_db(map=HRU, option='area', columns='hru_area_m2', units='meters', quiet=True)
+
+ # GET MEAN VALUES FOR THESE NEXT ONES, ACROSS THE BASIN
+
+ # SLOPE (and aspect)
+ #####################
+ v.rast_stats(map=HRU, raster=slope, method='average', column_prefix='tmp', flags='c', quiet=True)
+ v.db_update(map=HRU, column='hru_slope', query_column='tmp_average', quiet=True)
+
+ # ASPECT
+ #########
+ v.db_dropcolumn(map=HRU, columns='tmp_average', quiet=True)
+ # Dealing with conversion from degrees (no good average) to something I can
+ # average -- x- and y-vectors
+ # Geographic coordinates, so sin=x, cos=y.... not that it matters so long
+ # as I am consistent in how I return to degrees
+ r.mapcalc('aspect_x = sin(' + aspect + ')', overwrite=gscript.overwrite(), quiet=True)
+ r.mapcalc('aspect_y = cos(' + aspect + ')', overwrite=gscript.overwrite(), quiet=True)
+ #grass.run_command('v.db.addcolumn', map=HRU, columns='aspect_x_sum double precision, aspect_y_sum double precision, ncells_in_hru integer')
+ v.rast_stats(map=HRU, raster='aspect_x', method='sum', column_prefix='aspect_x', flags='c', quiet=True)
+ v.rast_stats(map=HRU, raster='aspect_y', method='sum', column_prefix='aspect_y', flags='c', quiet=True)
+ hru = VectorTopo(HRU)
+ hru.open('rw')
+ cur = hru.table.conn.cursor()
+ cur.execute("SELECT cat,aspect_x_sum,aspect_y_sum FROM %s" %hru.name)
+ _arr = np.array(cur.fetchall()).astype(float)
+ _cat = _arr[:,0]
+ _aspect_x_sum = _arr[:,1]
+ _aspect_y_sum = _arr[:,2]
+ aspect_angle = np.arctan2(_aspect_y_sum, _aspect_x_sum) * 180. / np.pi
+ aspect_angle[aspect_angle < 0] += 360 # all positive
+ aspect_angle_cat = np.vstack((aspect_angle, _cat)).transpose()
+ cur.executemany("update "+ HRU +" set hru_aspect=? where cat=?", aspect_angle_cat)
+ hru.table.conn.commit()
+ hru.close()
+
+ # ELEVATION
+ ############
+ v.rast_stats(map=HRU, raster=elevation, method='average', column_prefix='tmp', flags='c', quiet=True)
+ v.db_update(map=HRU, column='hru_elev', query_column='tmp_average', quiet=True)
+ v.db_dropcolumn(map=HRU, columns='tmp_average', quiet=True)
+
+ # CENTROIDS -- NEED FIXING FOR TRUE CENTERS!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ############
+
+ # get x,y of centroid -- but have areas not in database table, that do have
+ # centroids, and having a hard time finding a good way to get rid of them!
+ # They have duplicate category values!
+ # Perhaps these are little dangles on the edges of the vectorization where
+ # the raster value was the same but pinched out into 1-a few cells?
+ # From looking at map, lots of extra centroids on area boundaries, and removing
+ # small areas (though threshold hard to guess) gets rid of these
+
+ v.db_addcolumn(map=HRU, columns='centroid_x double precision, centroid_y double precision', quiet=True)
+ v.to_db(map=HRU, type='centroid', columns='centroid_x, centroid_y', option='coor', units='meters', quiet=True)
+
+ # hru_columns.append('hru_lat double precision') # Latitude of centroid
+ colNames = np.array(gscript.vector_db_select(HRU, layer=1)['columns'])
+ colValues = np.array(gscript.vector_db_select(HRU, layer=1)['values'].values())
+ xy = colValues[:,(colNames=='centroid_x') + (colNames=='centroid_y')]
+ np.savetxt('_xy.txt', xy, delimiter='|', fmt='%s')
+ m.proj(flags='od', input='_xy.txt', output='_lonlat.txt', overwrite=True, quiet=True)
+ lonlat = np.genfromtxt('_lonlat.txt', delimiter='|',)[:,:2]
+ lonlat_cat = np.concatenate((lonlat, np.expand_dims(_cat, 2)), axis=1)
+
+ # why not just get lon too?
+ v.db_addcolumn(map=HRU, columns='hru_lon double precision')
+
+ hru = VectorTopo(HRU)
+ hru.open('rw')
+ cur = hru.table.conn.cursor()
+ cur.executemany("update "+ HRU +" set hru_lon=?, hru_lat=? where cat=?", lonlat_cat)
+ hru.table.conn.commit()
+ hru.close()
+
+ # FIX!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ # CENTROIDS ARE NOT REALLY CENTERS, AND DO NOT WORK IF THERE ARE MULTIPLE
+ # AREAS
+ # Pygrass to find area edges / cells, and get the center?
+ # https://grass.osgeo.org/grass70/manuals/libpython/pygrass.vector.html
+
+ # Easting and Northing for other columns
+ v.db_update(map=HRU, column='hru_x', query_column='centroid_x', quiet=True)
+ v.db_update(map=HRU, column='hru_xlong', query_column='centroid_x*3.28084', quiet=True) # feet
+ v.db_update(map=HRU, column='hru_y', query_column='centroid_y', quiet=True)
+ v.db_update(map=HRU, column='hru_ylat', query_column='centroid_y*3.28084', quiet=True) # feet
+
+ # ID NUMBER
+ ############
+
+ """
+ hru = VectorTopo(HRU)
+ hru.open('rw')
+ cur = hru.table.conn.cursor()
+ cur.executemany("update "+HRU+" set hru_segment=? where id=?", hru_segmentt)
+ hru.table.conn.commit()
+ hru.close()
+ """
+ # Segment number = HRU ID number
+ v.db_update(map=HRU, column='hru_segment', query_column='id', quiet=True)
+
+
+if __name__ == "__main__":
+ main()
Added: grass-addons/grass7/vector/v.gsflow.reaches/Makefile
===================================================================
--- grass-addons/grass7/vector/v.gsflow.reaches/Makefile (rev 0)
+++ grass-addons/grass7/vector/v.gsflow.reaches/Makefile 2017-09-28 06:08:49 UTC (rev 71511)
@@ -0,0 +1,8 @@
+MODULE_TOPDIR = ../..
+
+PGM = v.gsflow.reaches
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script
+
Added: grass-addons/grass7/vector/v.gsflow.reaches/v.gsflow.reaches.html
===================================================================
--- grass-addons/grass7/vector/v.gsflow.reaches/v.gsflow.reaches.html (rev 0)
+++ grass-addons/grass7/vector/v.gsflow.reaches/v.gsflow.reaches.html 2017-09-28 06:08:49 UTC (rev 71511)
@@ -0,0 +1,25 @@
+<h2>DESCRIPTION</h2>
+
+<em>v.gsflow.reaches</em> produces channel "reaches", which are the intersection between PRMS stream segments and MODFLOW grid cells, for the USGS combined groundwater and surface-water model GSFLOW.
+
+<h2>REFERENCES</h2>
+
+Ng, G.-H. C., A. D. Wickert, R. L. McLaughlin, J. La Frenierre, S. Liess, and L. Sabeeri (2016),
+Modeling the role of groundwater and vegetation in the hydrological response of tropical
+glaciated watersheds to climate change, in AGU Fall Meeting Abstracts, H13L–1590, San
+Francisco, CA.
+
+<h2>SEE ALSO</h2>
+<a href="v.gsflow.export">v.gsflow.export</a>
+<a href="v.gsflow.gravres">v.gsflow.gravres</a>
+<a href="v.gsflow.grid">v.gsflow.grid</a>
+<a href="v.gsflow.hruparams">v.gsflow.hruparams</a>
+<a href="v.gsflow.segments">v.gsflow.segments</a>
+<a href="v.stream.inbasin">v.stream.inbasin</a>
+<a href="v.stream.network">v.stream.network</a>
+
+<h2>AUTHORS</h2>
+
+Andrew D. Wickert<br>
+
+<p><i>Last changed: $Date 2016-09-27$</i>
Added: grass-addons/grass7/vector/v.gsflow.reaches/v.gsflow.reaches.py
===================================================================
--- grass-addons/grass7/vector/v.gsflow.reaches/v.gsflow.reaches.py (rev 0)
+++ grass-addons/grass7/vector/v.gsflow.reaches/v.gsflow.reaches.py 2017-09-28 06:08:49 UTC (rev 71511)
@@ -0,0 +1,354 @@
+#!/usr/bin/env python
+############################################################################
+#
+# MODULE: v.gsflow.reaches
+#
+# AUTHOR(S): Andrew Wickert
+#
+# PURPOSE: Build reaches (intersection of PRMS Segments and MODFLOW
+# grid cells) for the GSFLOW (all USGS models)
+#
+# COPYRIGHT: (c) 2016-2017 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:
+# - uses inputs from r.stream.extract
+
+# More information
+# Started December 2016
+
+#%module
+#% description: Build stream "reaches" that link PRMS segments to MODFLOW cells
+#% keyword: vector
+#% keyword: stream network
+#% keyword: hydrology
+#% keyword: GSFLOW
+#%end
+
+#%option G_OPT_V_INPUT
+#% key: segment_input
+#% label: PRMS stream segments
+#% required: yes
+#% guidependency: layer,column
+#%end
+
+#%option G_OPT_V_INPUT
+#% key: grid_input
+#% label: MODFLOW grid
+#% required: yes
+#% guidependency: layer,column
+#%end
+
+#%option G_OPT_R_INPUT
+#% key: elevation
+#% label: DEM for slope along reaches
+#% required: yes
+#% guidependency: layer,column
+#%end
+
+#%option G_OPT_V_OUTPUT
+#% key: output
+#% label: Reaches for GSFLOW
+#% required: yes
+#% guidependency: layer,column
+#%end
+
+#%option
+#% key: s_min
+#% type: double
+#% description: Minimum reach slope
+#% answer: 0.0001
+#% required: no
+#%end
+
+#%option
+#% key: h_stream
+#% type: double
+#% description: Stream channel depth (bank height) [m]
+#% answer: 1
+#% required: no
+#%end
+
+#%option
+#% key: upstream_easting_column_seg
+#% type: string
+#% description: Upstream easting (or x or lon) column name
+#% answer: x1
+#% required : no
+#%end
+
+#%option
+#% key: upstream_northing_column_seg
+#% type: string
+#% description: Upstream northing (or y or lat) column name
+#% answer: y1
+#% required : no
+#%end
+
+#%option
+#% key: downstream_easting_column_seg
+#% type: string
+#% description: Downstream easting (or x or lon) column name
+#% answer: x2
+#% required : no
+#%end
+
+#%option
+#% key: downstream_northing_column_seg
+#% type: string
+#% description: Downstream northing (or y or lat) column name
+#% answer: y2
+#% required : no
+#%end
+
+#%option
+#% key: tostream_cat_column_seg
+#% type: string
+#% description: Adjacent downstream segment cat (0 = offmap)
+#% answer: tostream
+#% required : no
+#%end
+
+
+##################
+# IMPORT MODULES #
+##################
+# PYTHON
+import numpy as np
+from matplotlib import pyplot as plt
+import sys
+import warnings
+# GRASS
+from grass.pygrass.modules.shortcuts import general as g
+from grass.pygrass.modules.shortcuts import raster as r
+from grass.pygrass.modules.shortcuts import vector as v
+from grass.pygrass.modules.shortcuts import miscellaneous as m
+from grass.pygrass.gis import region
+from grass.pygrass import vector # Change to "v"?
+from grass.script import vector_db_select
+from grass.pygrass.vector import Vector, VectorTopo
+from grass.pygrass.raster import RasterRow
+from grass.pygrass import utils
+from grass import script as gscript
+
+###############
+# MAIN MODULE #
+###############
+
+def main():
+ """
+ Builds river reaches for input to the USGS hydrologic model, GSFLOW.
+ These reaches link the PRMS stream segments to the MODFLOW grid cells.
+ """
+
+ options, flags = gscript.parser()
+ segments = options['segment_input']
+ grid = options['grid_input']
+ reaches = options['output']
+ elevation = options['elevation']
+ Smin = options['s_min']
+ h_stream = options['h_stream']
+ x1 = options['upstream_easting_column_seg']
+ y1 = options['upstream_northing_column_seg']
+ x2 = options['downstream_easting_column_seg']
+ y2 = options['downstream_northing_column_seg']
+ tostream = options['tostream_cat_column_seg']
+
+ # Build reach maps by overlaying segments on grid
+ if len(gscript.find_file(segments, element='vector')['name']) > 0:
+ v.extract(input=segments, output='GSFLOW_TEMP__', type='line', quiet=True, overwrite=True)
+ v.overlay(ainput='GSFLOW_TEMP__', atype='line', binput=grid, output=reaches, operator='and', overwrite=gscript.overwrite(), quiet=True)
+ g.remove(type='vector', name='GSFLOW_TEMP__', quiet=True, flags='f')
+ else:
+ gscript.fatal('No vector file "'+segments+'" found.')
+
+ # Start editing database table
+ reachesTopo = VectorTopo(reaches)
+ reachesTopo.open('rw')
+
+ # Rename a,b columns
+ reachesTopo.table.columns.rename('a_'+x1, 'x1')
+ reachesTopo.table.columns.rename('a_'+x2, 'x2')
+ reachesTopo.table.columns.rename('a_'+y1, 'y1')
+ reachesTopo.table.columns.rename('a_'+y2, 'y2')
+ reachesTopo.table.columns.rename('a_NSEG', 'NSEG')
+ reachesTopo.table.columns.rename('a_ISEG', 'ISEG')
+ reachesTopo.table.columns.rename('a_stream_type', 'stream_type')
+ reachesTopo.table.columns.rename('a_type_code', 'type_code')
+ reachesTopo.table.columns.rename('a_cat', 'rnum_cat')
+ reachesTopo.table.columns.rename('a_'+tostream, 'tostream')
+ reachesTopo.table.columns.rename('a_id', 'segment_id')
+ reachesTopo.table.columns.rename('a_OUTSEG', 'OUTSEG')
+ reachesTopo.table.columns.rename('b_row', 'row')
+ reachesTopo.table.columns.rename('b_col', 'col')
+ reachesTopo.table.columns.rename('b_id', 'cell_id')
+
+ # Drop unnecessary columns
+ cols = reachesTopo.table.columns.names()
+ for col in cols:
+ if (col[:2] == 'a_') or (col[:2] == 'b_'):
+ reachesTopo.table.columns.drop(col)
+
+ # Add new columns to 'reaches'
+ reachesTopo.table.columns.add('KRCH', 'integer')
+ reachesTopo.table.columns.add('IRCH', 'integer')
+ reachesTopo.table.columns.add('JRCH', 'integer')
+ reachesTopo.table.columns.add('IREACH', 'integer')
+ reachesTopo.table.columns.add('RCHLEN', 'integer')
+ reachesTopo.table.columns.add('STRTOP', 'double precision')
+ reachesTopo.table.columns.add('SLOPE', 'double precision')
+ reachesTopo.table.columns.add('STRTHICK', 'double precision')
+ reachesTopo.table.columns.add('xr1', 'double precision')
+ reachesTopo.table.columns.add('xr2', 'double precision')
+ reachesTopo.table.columns.add('yr1', 'double precision')
+ reachesTopo.table.columns.add('yr2', 'double precision')
+
+ # Commit columns before editing (necessary?)
+ reachesTopo.table.conn.commit()
+ reachesTopo.close()
+
+ # Update some columns that can be done now
+ reachesTopo.open('rw')
+ colNames = np.array(gscript.vector_db_select(reaches, layer=1)['columns'])
+ colValues = np.array(gscript.vector_db_select(reaches, layer=1)['values'].values())
+ cats = colValues[:,colNames == 'cat'].astype(int).squeeze()
+ nseg = np.arange(1, len(cats)+1)
+ nseg_cats = []
+ for i in range(len(cats)):
+ nseg_cats.append( (nseg[i], cats[i]) )
+ cur = reachesTopo.table.conn.cursor()
+ cur.execute("update "+reaches+" set KRCH=1") # MAKE A VARIABLE LATER???
+ cur.execute("update "+reaches+" set STRTHICK=0.1") # 10 cm, prescribed -- MAKE A VARIABLE LATER!!!!!
+ cur.executemany("update "+reaches+" set IRCH=? where row=?", nseg_cats)
+ cur.executemany("update "+reaches+" set JRCH=? where col=?", nseg_cats)
+ reachesTopo.table.conn.commit()
+ reachesTopo.close()
+ v.to_db(map=reaches, columns='RCHLEN', option='length', quiet=True)
+
+
+ # Still to go after these:
+ # STRTOP (added with slope)
+ # IREACH (whole next section dedicated to this)
+ # SLOPE (need z_start and z_end)
+
+ # Now, the light stuff is over: time to build the reach order
+ v.to_db(map=reaches, option='start', columns='xr1,yr1')
+ v.to_db(map=reaches, option='end', columns='xr2,yr2')
+
+ # Now just sort by category, find which stream has the same xr1 and yr1 as
+ # x1 and y1 (or a_x1, a_y1) and then find where its endpoint matches another
+ # starting point and move down the line.
+ # v.db.select reaches col=cat,a_id,xr1,xr2 where="a_x1 = xr1"
+
+ # First, get the starting coordinates of each stream segment
+ # and a set of river ID's (ordered from 1...N)
+ colNames = np.array(gscript.vector_db_select(segments, layer=1)['columns'])
+ colValues = np.array(gscript.vector_db_select(segments, layer=1)['values'].values())
+ number_of_segments = colValues.shape[0]
+ segment_x1s = colValues[:,colNames == 'x1'].astype(float).squeeze()
+ segment_y1s = colValues[:,colNames == 'y1'].astype(float).squeeze()
+ segment_ids = colValues[:,colNames == 'id'].astype(float).squeeze()
+
+ # Then move back to the reaches map to produce the ordering
+ colNames = np.array(gscript.vector_db_select(reaches, layer=1)['columns'])
+ colValues = np.array(gscript.vector_db_select(reaches, layer=1)['values'].values())
+ reach_cats = colValues[:,colNames == 'cat'].astype(int).squeeze()
+ reach_x1s = colValues[:,colNames == 'xr1'].astype(float).squeeze()
+ reach_y1s = colValues[:,colNames == 'yr1'].astype(float).squeeze()
+ reach_x2s = colValues[:,colNames == 'xr2'].astype(float).squeeze()
+ reach_y2s = colValues[:,colNames == 'yr2'].astype(float).squeeze()
+ segment_ids__reach = colValues[:,colNames == 'segment_id'].astype(float).squeeze()
+
+ for segment_id in segment_ids:
+ reach_order_cats = []
+ downstream_directed = []
+ ssel = segment_ids == segment_id
+ rsel = segment_ids__reach == segment_id # selector
+ # Find first segment: x1y1 first here, but not necessarily later
+ downstream_directed.append(1)
+ _x_match = reach_x1s[rsel] == segment_x1s[ssel]
+ _y_match = reach_y1s[rsel] == segment_y1s[ssel]
+ _i_match = _x_match * _y_match
+ x1y1 = True # false if x2y2
+ # Find cat
+ _cat = int(reach_cats[rsel][_x_match * _y_match])
+ reach_order_cats.append(_cat)
+ # Get end of reach = start of next one
+ reach_x_end = float(reach_x2s[reach_cats == _cat])
+ reach_y_end = float(reach_y2s[reach_cats == _cat])
+ while _i_match.any():
+ _x_match = reach_x1s[rsel] == reach_x_end
+ _y_match = reach_y1s[rsel] == reach_y_end
+ _i_match = _x_match * _y_match
+ if _i_match.any():
+ _cat = int(reach_cats[rsel][_x_match * _y_match])
+ reach_x_end = float(reach_x2s[reach_cats == _cat])
+ reach_y_end = float(reach_y2s[reach_cats == _cat])
+ reach_order_cats.append(_cat)
+ print len(reach_order_cats), len(reach_cats[rsel])
+
+ # Reach order to database table
+ reach_number__reach_order_cats = []
+ for i in range(len(reach_order_cats)):
+ reach_number__reach_order_cats.append( (i+1, reach_order_cats[i]) )
+ reachesTopo = VectorTopo(reaches)
+ reachesTopo.open('rw')
+ cur = reachesTopo.table.conn.cursor()
+ cur.executemany("update "+reaches+" set IREACH=? where cat=?",
+ reach_number__reach_order_cats)
+ reachesTopo.table.conn.commit()
+ reachesTopo.close()
+
+
+ # TOP AND BOTTOM ARE OUT OF ORDER: SOME SEGS ARE BACKWARDS. UGH!!!!
+ # NEED TO GET THEM IN ORDER TO GET THE Z VALUES AT START AND END
+
+ # Compute slope and starting elevations from the elevations at the start and
+ # end of the reaches and the length of each reach]
+ gscript.message('Obtaining elevation values from raster: may take time.')
+ v.db_addcolumn(map=reaches, columns='zr1 double precision, zr2 double precision')
+ zr1 = []
+ zr2 = []
+ for i in range(len(reach_cats)):
+ _x = reach_x1s[i]
+ _y = reach_y1s[i]
+ _z = float(gscript.parse_command('r.what', map=elevation, coordinates=str(_x)+','+str(_y)).keys()[0].split('|')[-1])
+ zr1.append(_z)
+ _x = reach_x2s[i]
+ _y = reach_y2s[i]
+ _z = float(gscript.parse_command('r.what', map=elevation, coordinates=str(_x)+','+str(_y)).keys()[0].split('|')[-1])
+ zr2.append(_z)
+
+ zr1_cats = []
+ zr2_cats = []
+ for i in range(len(reach_cats)):
+ zr1_cats.append( (zr1[i], reach_cats[i]) )
+ zr2_cats.append( (zr2[i], reach_cats[i]) )
+
+ reachesTopo = VectorTopo(reaches)
+ reachesTopo.open('rw')
+ cur = reachesTopo.table.conn.cursor()
+ cur.executemany("update "+reaches+" set zr1=? where cat=?", zr1_cats)
+ cur.executemany("update "+reaches+" set zr2=? where cat=?", zr2_cats)
+ reachesTopo.table.conn.commit()
+ reachesTopo.close()
+
+ # Use these to create slope -- backwards possible on DEM!
+ v.db_update(map=reaches, column='SLOPE', value='(zr1 - zr2)/RCHLEN')
+ v.db_update(map=reaches, column='SLOPE', value=Smin, where='SLOPE <= '+str(Smin))
+
+ ## srtm_local_filled_grid = srtm_local_filled @ 200m (i.e. current grid)
+ ## resolution
+ ## r.to.vect in=srtm_local_filled_grid out=srtm_local_filled_grid col=z type=area --o#
+ #v.db_addcolumn(map=reaches, columns='z_topo_mean double precision')
+ #v.what_vect(map=reaches, query_map='srtm_local_filled_grid', column='z_topo_mean', query_column='z')
+ #v.db_update(map=reaches, column='STRTOP', value='z_topo_mean -'+str(h_stream))
+
+
+if __name__ == "__main__":
+ main()
Added: grass-addons/grass7/vector/v.gsflow.segments/Makefile
===================================================================
--- grass-addons/grass7/vector/v.gsflow.segments/Makefile (rev 0)
+++ grass-addons/grass7/vector/v.gsflow.segments/Makefile 2017-09-28 06:08:49 UTC (rev 71511)
@@ -0,0 +1,8 @@
+MODULE_TOPDIR = ../..
+
+PGM = v.gsflow.segments
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script
+
Added: grass-addons/grass7/vector/v.gsflow.segments/v.gsflow.segments.html
===================================================================
--- grass-addons/grass7/vector/v.gsflow.segments/v.gsflow.segments.html (rev 0)
+++ grass-addons/grass7/vector/v.gsflow.segments/v.gsflow.segments.html 2017-09-28 06:08:49 UTC (rev 71511)
@@ -0,0 +1,26 @@
+<h2>DESCRIPTION</h2>
+
+<em>v.gsflow.segments</em> produces channel segments that are each one link in the full drainage network. Each segment also corresponds to a hydrologic response unit (HRU) in this formulation. These are used for the USGS combined groundwater (MODFLOW) and surface-water (PRMS) model GSFLOW.
+
+<h2>REFERENCES</h2>
+
+Ng, G.-H. C., A. D. Wickert, R. L. McLaughlin, J. La Frenierre, S. Liess, and L. Sabeeri (2016),
+Modeling the role of groundwater and vegetation in the hydrological response of tropical
+glaciated watersheds to climate change, in AGU Fall Meeting Abstracts, H13L–1590, San
+Francisco, CA.
+
+<h2>SEE ALSO</h2>
+<a href="v.gsflow.export">v.gsflow.export</a>
+<a href="v.gsflow.gravres">v.gsflow.gravres</a>
+<a href="v.gsflow.grid">v.gsflow.grid</a>
+<a href="v.gsflow.hruparams">v.gsflow.hruparams</a>
+<a href="v.gsflow.reaches">v.gsflow.reaches</a>
+<a href="v.gsflow.segments">v.gsflow.segments</a>
+<a href="v.stream.inbasin">v.stream.inbasin</a>
+<a href="v.stream.network">v.stream.network</a>
+
+<h2>AUTHORS</h2>
+
+Andrew D. Wickert<br>
+
+<p><i>Last changed: $Date 2016-09-27$</i>
Added: grass-addons/grass7/vector/v.gsflow.segments/v.gsflow.segments.py
===================================================================
--- grass-addons/grass7/vector/v.gsflow.segments/v.gsflow.segments.py (rev 0)
+++ grass-addons/grass7/vector/v.gsflow.segments/v.gsflow.segments.py 2017-09-28 06:08:49 UTC (rev 71511)
@@ -0,0 +1,327 @@
+#!/usr/bin/env python
+############################################################################
+#
+# MODULE: v.gsflow.segments
+#
+# AUTHOR(S): Andrew Wickert
+#
+# PURPOSE: Build stream segments for the USGS models MODFLOW or GSFLOW.
+#
+# COPYRIGHT: (c) 2016-2017 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:
+# - uses inputs from v.stream.network
+
+# More information
+# Started December 2016
+
+#%module
+#% description: Prepares stream segments for PRMS and GSFLOW
+#% keyword: vector
+#% keyword: stream network
+#% keyword: hydrology
+#% keyword: GSFLOW
+#%end
+
+#%option G_OPT_V_INPUT
+#% key: input
+#% label: Vector stream network from r.stream.extract
+#% required: yes
+#% guidependency: layer,column
+#%end
+
+#%option G_OPT_V_OUTPUT
+#% key: output
+#% label: Segments: stream segments for GSFLOW / PRMS
+#% required: yes
+#% guidependency: layer,column
+#%end
+
+#%option
+#% key: ICALC
+#% type: integer
+#% description: Stream depth option: 0-const; 1,2-Manning; 3-aQ^b
+#% answer: 3
+#% required: no
+#%end
+
+#%option
+#% key: CDPTH
+#% type: double
+#% description: Flow depth coefficient; used if ICALC=3
+#% answer: 0.25
+#% required: no
+#%end
+
+#%option
+#% key: FDPTH
+#% type: double
+#% description: Flow depth exponent; used if ICALC=3
+#% answer: 0.4
+#% required: no
+#%end
+
+#%option
+#% key: AWDTH
+#% type: double
+#% description: Flow width coefficient; used if ICALC=3
+#% answer: 8
+#% required: no
+#%end
+
+#%option
+#% key: BWDTH
+#% type: double
+#% description: Flow width exponent; used if ICALC=3
+#% answer: 0.5
+#% required: no
+#%end
+
+#%option
+#% key: IUPSEG
+#% type: string
+#% description: Category of upstream diversion segment (from_cat,to_cat,...)
+#% answer: 0,0
+#% required: no
+#%end
+
+#%option
+#% key: FLOW
+#% type: string
+#% description: Streamflow entering the upstream-most segments (cat,Q,cat,Q,...)
+#% answer: 0,0
+#% required: no
+#%end
+
+#%option
+#% key: RUNOFF
+#% type: string
+#% description: Diffuse runoff entering each segment (cat,Q,cat,Q,...)
+#% answer: 0,0
+#% required: no
+#%end
+
+#%option
+#% key: ETSW
+#% type: string
+#% description: Direct removal of in-channel water by ET (cat,Q,cat,Q)
+#% answer: 0,0
+#% required: no
+#%end
+
+#%option
+#% key: PPTSW
+#% type: string
+#% description: Direct precipitation on the stream
+#% answer: 0,0
+#% required: no
+#%end
+
+#%option
+#% key: PPTSW
+#% type: string
+#% description: Direct precipitation on the stream
+#% answer: 0,0
+#% required: no
+#%end
+
+#%option
+#% key: ROUGHCH
+#% type: double
+#% description: In-channel Manning's n for ICALC=1,2
+#% answer: 0.035
+#% required: no
+#%end
+
+#%option
+#% key: ROUGHBK
+#% type: double
+#% description: Overbank Manning's n for ICALC=2
+#% answer: 0.06
+#% required: no
+#%end
+
+#%option
+#% key: WIDTH1
+#% type: double
+#% description: Upstream width in segment, assumed constant through watershed
+#% answer: 5
+#% required: no
+#%end
+
+#%option
+#% key: WIDTH2
+#% type: double
+#% description: Downstream width in segment, assumed constant through watershed
+#% answer: 5
+#% required: no
+#%end
+
+##################
+# IMPORT MODULES #
+##################
+# PYTHON
+import numpy as np
+from matplotlib import pyplot as plt
+import sys
+# GRASS
+from grass.pygrass.modules.shortcuts import general as g
+from grass.pygrass.modules.shortcuts import raster as r
+from grass.pygrass.modules.shortcuts import vector as v
+from grass.pygrass.modules.shortcuts import miscellaneous as m
+from grass.pygrass.gis import region
+from grass.pygrass import vector # Change to "v"?
+from grass.script import vector_db_select
+from grass.pygrass.vector import Vector, VectorTopo
+from grass.pygrass.raster import RasterRow
+from grass.pygrass import utils
+from grass import script as gscript
+
+###############
+# MAIN MODULE #
+###############
+
+def main():
+ """
+ Builds river segments for input to the USGS hydrologic models
+ PRMS and GSFLOW.
+ """
+
+ ##################
+ # OPTION PARSING #
+ ##################
+
+ options, flags = gscript.parser()
+
+ # I/O
+ streams = options['input']
+ segments = options['output']
+
+ # Hydraulic geometry
+ ICALC = options['ICALC']
+
+ # ICALC=0: Constant depth
+ WIDTH1 = options['WIDTH1']
+ WIDTH2 = options['WIDTH2']
+
+ # ICALC=1: Manning
+ ROUGHCH = options['ROUGHCH']
+
+ # ICALC=2: Manning
+ ROUGHBK = options['ROUGHBK']
+
+ # ICALC=3: Power-law relationships (following Leopold and others)
+ CDPTH = options['CDPTH']
+ FDPTH = options['FDPTH']
+ AWDTH = options['AWDTH']
+ BWDTH = options['BWDTH']
+
+ ##################################################
+ # CHECKING DEPENDENCIES WITH OPTIONAL PARAMETERS #
+ ##################################################
+
+ if ICALC == 3:
+ if CDPTH and FDPTH and AWDTH and BWDTH:
+ pass
+ else:
+ grass.fatal('Missing CDPTH, FDPTH, AWDTH, and/or BWDTH. \
+ These are required when ICALC = 3.')
+
+ ###########
+ # RUNNING #
+ ###########
+
+ # New Columns for Segments
+ segment_columns = []
+ # Self ID
+ segment_columns.append('id integer') # segment number
+ segment_columns.append('ISEG integer') # segment number
+ segment_columns.append('NSEG integer') # segment number
+ # for GSFLOW
+ segment_columns.append('ICALC integer') # 3 for power function
+ segment_columns.append('OUTSEG integer') # downstream segment -- tostream, renumbered
+ segment_columns.append('ROUGHCH double precision') # overbank roughness
+ segment_columns.append('ROUGHBK double precision') # in-channel roughness
+ segment_columns.append('WIDTH1 double precision') # overbank roughness
+ segment_columns.append('WIDTH2 double precision') # in-channel roughness
+ segment_columns.append('CDPTH double precision') # depth coeff
+ segment_columns.append('FDPTH double precision') # depth exp
+ segment_columns.append('AWDTH double precision') # width coeff
+ segment_columns.append('BWDTH double precision') # width exp
+ # The below will be all 0
+ segment_columns.append('IUPSEG integer') # upstream segment ID number, for diversions
+ segment_columns.append('FLOW integer')
+ segment_columns.append('RUNOFF integer')
+ segment_columns.append('ETSW integer')
+ segment_columns.append('PPTSW integer')
+
+ segment_columns = ",".join(segment_columns)
+
+ # CONSIDER THE EFFECT OF OVERWRITING COLUMNS -- WARN FOR THIS
+ # IF MAP EXISTS ALREADY?
+
+ # Create a map to work with
+ g.copy(vector=(streams, segments), overwrite=gscript.overwrite())
+ # and add its columns
+ v.db_addcolumn(map=segments, columns=segment_columns)
+
+ # Produce the data table entries
+ ##################################
+ colNames = np.array(gscript.vector_db_select(segments, layer=1)['columns'])
+ colValues = np.array(gscript.vector_db_select(segments, layer=1)['values'].values())
+ number_of_segments = colValues.shape[0]
+ cats = colValues[:,colNames == 'cat'].astype(int).squeeze()
+
+ nseg = np.arange(1, len(cats)+1)
+ nseg_cats = []
+ for i in range(len(cats)):
+ nseg_cats.append( (nseg[i], cats[i]) )
+
+ segmentsTopo = VectorTopo(segments)
+ segmentsTopo.open('rw')
+ cur = segmentsTopo.table.conn.cursor()
+
+ # id = cat (as does ISEG and NSEG)
+ cur.executemany("update "+segments+" set id=? where cat=?", nseg_cats)
+ cur.executemany("update "+segments+" set ISEG=? where cat=?", nseg_cats)
+ cur.executemany("update "+segments+" set NSEG=? where cat=?", nseg_cats)
+
+ # outseg = tostream
+ cur.executemany("update "+segments+" set OUTSEG=? where tostream=?", nseg_cats)
+
+ # Discharge and hydraulic geometry
+ cur.execute("update "+segments+" set WIDTH1="+str(WIDTH1))
+ cur.execute("update "+segments+" set WIDTH2="+str(WIDTH2))
+ cur.execute("update "+segments+" set ROUGHCH="+str(ROUGHCH))
+ cur.execute("update "+segments+" set ROUGHBK="+str(ROUGHBK))
+ cur.execute("update "+segments+" set ICALC="+str(ICALC))
+ cur.execute("update "+segments+" set CDPTH="+str(CDPTH))
+ cur.execute("update "+segments+" set FDPTH="+str(FDPTH))
+ cur.execute("update "+segments+" set AWDTH="+str(AWDTH))
+ cur.execute("update "+segments+" set BWDTH="+str(BWDTH))
+
+ gscript.message('')
+ gscript.message('NOTICE: not currently used:')
+ gscript.message('IUPSEG, FLOW, RUNOFF, ETSW, and PPTSW.')
+ gscript.message('All set to 0.')
+ gscript.message('')
+
+ # values that are 0
+ cur.execute("update "+segments+" set IUPSEG="+str(0))
+ cur.execute("update "+segments+" set FLOW="+str(0))
+ cur.execute("update "+segments+" set RUNOFF="+str(0))
+ cur.execute("update "+segments+" set ETSW="+str(0))
+ cur.execute("update "+segments+" set PPTSW="+str(0))
+
+ segmentsTopo.table.conn.commit()
+ segmentsTopo.close()
+
+
+if __name__ == "__main__":
+ main()
Added: grass-addons/grass7/vector/v.stream.inbasin/Makefile
===================================================================
--- grass-addons/grass7/vector/v.stream.inbasin/Makefile (rev 0)
+++ grass-addons/grass7/vector/v.stream.inbasin/Makefile 2017-09-28 06:08:49 UTC (rev 71511)
@@ -0,0 +1,8 @@
+MODULE_TOPDIR = ../..
+
+PGM = v.stream.inbasin
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script
+
Added: grass-addons/grass7/vector/v.stream.inbasin/v.stream.inbasin.html
===================================================================
--- grass-addons/grass7/vector/v.stream.inbasin/v.stream.inbasin.html (rev 0)
+++ grass-addons/grass7/vector/v.stream.inbasin/v.stream.inbasin.html 2017-09-28 06:08:49 UTC (rev 71511)
@@ -0,0 +1,21 @@
+<h2>DESCRIPTION</h2>
+
+<em>v.stream.inbasin</em> uses the output of v.stream.network to select only those streams (and sub-basins) that are upstream of (and inclusive of) a selected link in the network. It is used as a step to develop GSFLOW model inputs for a watershed, but need not be exclusively used for that purpose.
+
+<h2>REFERENCES</h2>
+
+Ng, G.-H. C., A. D. Wickert, R. L. McLaughlin, J. La Frenierre, S. Liess, and L. Sabeeri (2016),
+Modeling the role of groundwater and vegetation in the hydrological response of tropical
+glaciated watersheds to climate change, in AGU Fall Meeting Abstracts, H13L–1590, San
+Francisco, CA.
+
+<h2>SEE ALSO</h2>
+<a href="v.gsflow.hruparams">v.gsflow.hruparams</a>
+<a href="v.gsflow.segments">v.gsflow.segments</a>
+<a href="v.stream.network">v.stream.network</a>
+
+<h2>AUTHORS</h2>
+
+Andrew D. Wickert<br>
+
+<p><i>Last changed: $Date 2016-09-27$</i>
Added: grass-addons/grass7/vector/v.stream.inbasin/v.stream.inbasin.py
===================================================================
--- grass-addons/grass7/vector/v.stream.inbasin/v.stream.inbasin.py (rev 0)
+++ grass-addons/grass7/vector/v.stream.inbasin/v.stream.inbasin.py 2017-09-28 06:08:49 UTC (rev 71511)
@@ -0,0 +1,173 @@
+#!/usr/bin/env python
+############################################################################
+#
+# MODULE: v.stream.inbasin
+#
+# AUTHOR(S): Andrew Wickert
+#
+# PURPOSE: Build a drainage basin from the subwatersheds of a river
+# network, based on the structure of the network.
+#
+# COPYRIGHT: (c) 2016 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:
+# - uses inputs from v.stream.network
+
+# More information
+# Started 14 October 2016
+
+#%module
+#% description: Subset a stream network into just one of its basins
+#% keyword: vector
+#% keyword: stream network
+#% keyword: basins
+#% keyword: hydrology
+#% keyword: geomorphology
+#%end
+
+#%option G_OPT_V_INPUT
+#% key: input_streams
+#% label: Stream network
+#% required: yes
+#%end
+
+#%option G_OPT_V_INPUT
+#% key: input_basins
+#% label: Subbasins built alongside stream network
+#% required: no
+#%end
+
+#%option
+#% key: cat
+#% label: Farthest downstream segment category
+#% required: yes
+#% guidependency: layer,column
+#%end
+
+#%option G_OPT_V_OUTPUT
+#% key: output_basin
+#% label: Vector output drainage basin
+#% required: no
+#%end
+
+#%option G_OPT_V_OUTPUT
+#% key: output_streams
+#% label: Streams within vector output drainage basin
+#% required: no
+#%end
+
+#%option G_OPT_V_OUTPUT
+#% key: output_pour_point
+#% label: Basin outlet
+#% required: no
+#%end
+
+##################
+# IMPORT MODULES #
+##################
+# PYTHON
+import numpy as np
+from matplotlib import pyplot as plt
+import sys
+# GRASS
+from grass.pygrass.modules.shortcuts import general as g
+from grass.pygrass.modules.shortcuts import raster as r
+from grass.pygrass.modules.shortcuts import vector as v
+from grass.pygrass.gis import region
+from grass.pygrass import vector # Change to "v"?
+from grass.script import vector_db_select
+from grass.pygrass.vector import Vector, VectorTopo
+from grass.pygrass.raster import RasterRow
+from grass.pygrass import utils
+from grass import script as gscript
+from grass.pygrass.vector.geometry import Point
+
+###############
+# MAIN MODULE #
+###############
+
+def find_upstream_segments(cats, tostream, cat):
+ """
+ Find segments immediately upstream of the given one
+ """
+ pass
+
+
+def main():
+ """
+ Links each river segment to the next downstream segment in a tributary
+ network by referencing its category (cat) number in a new column. "0"
+ means that the river exits the map.
+ """
+
+ options, flags = gscript.parser()
+
+ streams = options['input_streams']
+ basins = options['input_basins']
+ cat = options['cat']
+ output_basins = options['output_basin']
+ output_streams = options['output_streams']
+ output_pour_point = options['output_pour_point']
+
+ #print options
+ #print flags
+
+ # Attributes of streams
+ colNames = np.array(vector_db_select(streams)['columns'])
+ colValues = np.array(vector_db_select(streams)['values'].values())
+ tostream = colValues[:,colNames == 'tostream'].astype(int).squeeze()
+ cats = colValues[:,colNames == 'cat'].astype(int).squeeze() # = "fromstream"
+
+ # Find network
+ basincats = [cat] # start here
+ most_upstream_cats = [cat] # all of those for which new cats must be sought
+ while True:
+ if len(most_upstream_cats) == 0:
+ break
+ tmp = list(most_upstream_cats) # copy to a temp file: old values
+ most_upstream_cats = [] # Ready to accept new values
+ for ucat in tmp:
+ most_upstream_cats += list(cats[tostream == int(ucat)])
+ basincats += most_upstream_cats
+
+ basincats = list(set(list(basincats)))
+
+ basincats_str = ','.join(map(str, basincats))
+
+ # Many basins out -- need to use overwrite flag in future!
+ #SQL_OR = 'rnum = ' + ' OR rnum = '.join(map(str, basincats))
+ SQL_OR = 'cat = ' + ' OR cat = '.join(map(str, basincats))
+ if len(basins) > 0:
+ v.extract(input=basins, output=output_basins, where=SQL_OR, overwrite=gscript.overwrite(), quiet=True)
+ if len(streams) > 0:
+ v.extract(input=streams, output=output_streams, cats=basincats_str, overwrite=gscript.overwrite(), quiet=True)
+ if len(output_pour_point) > 0:
+ _pp = gscript.vector_db_select(map=streams, columns='x2,y2', where='cat='+str(cat))
+ _xy = np.squeeze(_pp['values'].values())
+ _x = float(_xy[0])
+ _y = float(_xy[1])
+ # NEED TO ADD IF-STATEMENT HERE TO AVOID AUTOMATIC OVERWRITING!!!!!!!!!!!
+ try:
+ v.db_droptable(table=output_pour_point, flags='f')
+ except:
+ pass
+ pptmp = vector.Vector(output_pour_point)
+ _cols = [(u'cat', 'INTEGER PRIMARY KEY'),
+ (u'x', 'DOUBLE PRECISION'),
+ (u'y', 'DOUBLE PRECISION')]
+ pptmp.open('w', tab_name=output_pour_point, tab_cols=_cols)
+ point0 = Point(_x,_y)
+ pptmp.write(point0, cat=1, attrs=(str(_x), str(_y)), )
+ pptmp.table.conn.commit()
+ pptmp.build()
+ pptmp.close()
+
+
+if __name__ == "__main__":
+ main()
Added: grass-addons/grass7/vector/v.stream.network/Makefile
===================================================================
--- grass-addons/grass7/vector/v.stream.network/Makefile (rev 0)
+++ grass-addons/grass7/vector/v.stream.network/Makefile 2017-09-28 06:08:49 UTC (rev 71511)
@@ -0,0 +1,8 @@
+MODULE_TOPDIR = ../..
+
+PGM = v.stream.network
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script
+
Added: grass-addons/grass7/vector/v.stream.network/v.stream.network.html
===================================================================
--- grass-addons/grass7/vector/v.stream.network/v.stream.network.html (rev 0)
+++ grass-addons/grass7/vector/v.stream.network/v.stream.network.html 2017-09-28 06:08:49 UTC (rev 71511)
@@ -0,0 +1,49 @@
+<h2>DESCRIPTION</h2>
+
+<em>v.stream.network</em> populates the attribute table of a river network built by <a href="r.stream.extract"><em>r.stream.extract</em></a> with five new columns that define the nodes and links within a river network. These are (default names given):
+<ul>
+ <li><em>x1</em>, <em>y1</em>: upstream coordinates</li>
+ <li><em>x2</em>, <em>y2</em>: downstream coordinates</li>
+ <li><em>tostream</em>: category number of the next stream; is equal to 0 if the stream flows off of the map</li>
+</ul>
+
+<h2>NOTES</h2>
+
+<b>streams</b> is a set of vector lines that is generated by r.stream.extract. It is recommended to be built as follows (Python code):
+
+<div class="code"><pre>
+# Choose elevation map and threshold
+elevation = srtm
+threshold = 1E6 # [m**2], minimum catchment size = 1 km**2
+# Build a map of cell areas: pick one of these
+# km2 is more intuitive, but 1 km2 is the threshold for the smallest
+# threshold catchment size
+r.mapcalc('cellArea_meters2 = '+str(reg.nsres)+' * '+str(reg.ewres), overwrite=True)
+r.mapcalc("cellArea_km2 = cellArea_meters2 / 10^6", overwrite=True)
+
+# Use r.watershed to build a drainage basin with units of length**2
+# instead of arbitrary units (cells)
+r.watershed(elevation=elevation, flow='cellArea_meters2', accumulation='drainageArea_m2', drainage='drainageDirection', stream='streams', threshold=thresh, flags='s', overwrite=True)
+# Note that this will include areas of negative (i.e. offmap) flow accumulation
+
+# Build watershed network using r.stream.extract: single-flow-direction (SFD)
+print "Building drainage network"
+r.stream_extract(elevation=elevation, accumulation='drainageArea_m2', threshold=thresh, d8cut=0, mexp=0, stream_raster='streams', stream_vector='streams', direction='draindir', overwrite=True)
+</pre></div>
+
+<h2>SEE ALSO</h2>
+
+<em>
+<a href="r.stream.extract">r.stream.extract</a>
+<a href="v.stream.order">v.stream.order</a>
+</em>
+
+<h2>REFERENCES</h2>
+
+None (yet)
+
+<h2>AUTHORS</h2>
+
+Andrew D. Wickert<br>
+
+<p><i>Last changed: $Date 2016-09-27$</i>
Added: grass-addons/grass7/vector/v.stream.network/v.stream.network.py
===================================================================
--- grass-addons/grass7/vector/v.stream.network/v.stream.network.py (rev 0)
+++ grass-addons/grass7/vector/v.stream.network/v.stream.network.py 2017-09-28 06:08:49 UTC (rev 71511)
@@ -0,0 +1,210 @@
+#!/usr/bin/env python
+############################################################################
+#
+# MODULE: v.stream.network
+#
+# AUTHOR(S): Andrew Wickert
+#
+# PURPOSE: Attach IDs of upstream and downstream nodes as well as the
+# category value of the next downstream stream segment
+# (0 if the stream exits the map)
+#
+# COPYRIGHT: (c) 2016-2017 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:
+# - uses inputs from r.stream.extract
+
+# More information
+# Started 14 October 2016
+
+#%module
+#% description: Build a linked stream network: each link knows its downstream link
+#% keyword: vector
+#% keyword: stream network
+#% keyword: hydrology
+#% keyword: geomorphology
+#%end
+
+#%option G_OPT_V_INPUT
+#% key: map
+#% label: Vector stream network from r.stream.extract
+#% required: yes
+#% guidependency: layer,column
+#%end
+
+#%option
+#% key: upstream_easting_column
+#% type: string
+#% description: Upstream easting (or x or lon) column name
+#% answer: x1
+#% required : no
+#%end
+
+#%option
+#% key: upstream_northing_column
+#% type: string
+#% description: Upstream northing (or y or lat) column name
+#% answer: y1
+#% required : no
+#%end
+
+#%option
+#% key: downstream_easting_column
+#% type: string
+#% description: Downstream easting (or x or lon) column name
+#% answer: x2
+#% required : no
+#%end
+
+#%option
+#% key: downstream_northing_column
+#% type: string
+#% description: Downstream northing (or y or lat) column name
+#% answer: y2
+#% required : no
+#%end
+
+#%option
+#% key: tostream_cat_column
+#% type: string
+#% description: Adjacent downstream segment cat (0 = offmap flow)
+#% answer: tostream
+#% required : no
+#%end
+
+##################
+# IMPORT MODULES #
+##################
+# PYTHON
+import numpy as np
+from matplotlib import pyplot as plt
+import sys
+# GRASS
+from grass.pygrass.modules.shortcuts import general as g
+from grass.pygrass.modules.shortcuts import raster as r
+from grass.pygrass.modules.shortcuts import vector as v
+from grass.pygrass.gis import region
+from grass.pygrass import vector # Change to "v"?
+from grass.script import vector_db_select
+from grass.pygrass.vector import Vector, VectorTopo
+from grass.pygrass.raster import RasterRow
+from grass.pygrass import utils
+from grass import script as gscript
+
+###############
+# MAIN MODULE #
+###############
+
+def main():
+ """
+ Links each river segment to the next downstream segment in a tributary
+ network by referencing its category (cat) number in a new column. "0"
+ means that the river exits the map.
+ """
+
+ options, flags = gscript.parser()
+ streams = options['map']
+ x1 = options['upstream_easting_column']
+ y1 = options['upstream_northing_column']
+ x2 = options['downstream_easting_column']
+ y2 = options['downstream_northing_column']
+
+ streamsTopo = VectorTopo(streams)
+ #streamsTopo.build()
+
+ # 1. Get vectorTopo
+ streamsTopo.open(mode='rw')
+ """
+ points_in_streams = []
+ cat_of_line_segment = []
+
+ # 2. Get coordinates
+ for row in streamsTopo:
+ cat_of_line_segment.append(row.cat)
+ if type(row) == vector.geometry.Line:
+ points_in_streams.append(row)
+ """
+
+ # 3. Coordinates of points: 1 = start, 2 = end
+ try:
+ streamsTopo.table.columns.add(x1,'double precision')
+ except:
+ pass
+ try:
+ streamsTopo.table.columns.add(y1,'double precision')
+ except:
+ pass
+ try:
+ streamsTopo.table.columns.add(x2,'double precision')
+ except:
+ pass
+ try:
+ streamsTopo.table.columns.add(y2,'double precision')
+ except:
+ pass
+ try:
+ streamsTopo.table.columns.add('tostream','int')
+ except:
+ pass
+ streamsTopo.table.conn.commit()
+
+ # Is this faster than v.to.db?
+ """
+ cur = streamsTopo.table.conn.cursor()
+ for i in range(len(points_in_streams)):
+ cur.execute("update streams set x1="+str(points_in_streams[i][0].x)+" where cat="+str(cat_of_line_segment[i]))
+ cur.execute("update streams set y1="+str(points_in_streams[i][0].y)+" where cat="+str(cat_of_line_segment[i]))
+ cur.execute("update streams set x2="+str(points_in_streams[i][-1].x)+" where cat="+str(cat_of_line_segment[i]))
+ cur.execute("update streams set y2="+str(points_in_streams[i][-1].y)+" where cat="+str(cat_of_line_segment[i]))
+ streamsTopo.table.conn.commit()
+ streamsTopo.build()
+ """
+ # v.to.db Works more consistently, at least
+ streamsTopo.close()
+ v.to_db(map=streams, option='start', columns=x1+','+y1)
+ v.to_db(map=streams, option='end', columns=x2+','+y2)
+
+ # 4. Read in and save the start and end coordinate points
+ colNames = np.array(vector_db_select(streams)['columns'])
+ colValues = np.array(vector_db_select(streams)['values'].values())
+ cats = colValues[:,colNames == 'cat'].astype(int).squeeze() # river number
+ xy1 = colValues[:,(colNames == 'x1') + (colNames == 'y1')].astype(float) # upstream
+ xy2 = colValues[:,(colNames == 'x2') + (colNames == 'y2')].astype(float) # downstream
+
+ # 5. Build river network
+ tocat = []
+ for i in range(len(cats)):
+ tosegment_mask = np.prod(xy1 == xy2[i], axis=1)
+ if np.sum(tosegment_mask) == 0:
+ tocat.append(0)
+ else:
+ tocat.append(cats[tosegment_mask.nonzero()[0][0]])
+ tocat = np.asarray(tocat).astype(int)
+
+ # This gives us a set of downstream-facing adjacencies.
+ # We will update the database with it.
+ streamsTopo.build()
+ streamsTopo.open('rw')
+ cur = streamsTopo.table.conn.cursor()\
+ # use "executemany" ????????????????????????????????????????????????????
+ # see v.gsflow.segments.py
+ for i in range(len(tocat)):
+ cur.execute("update "+streams+" set tostream="+str(tocat[i])+" where cat="+str(cats[i]))
+ streamsTopo.table.conn.commit()
+ #streamsTopo.build()
+ streamsTopo.close()
+
+ gscript.message('')
+ gscript.message('Drainage topology built. Check "tostream" column for the downstream cat.')
+ gscript.message('A cat value of 0 indicates the downstream-most segment.')
+ gscript.message('')
+
+
+if __name__ == "__main__":
+ main()
More information about the grass-commit
mailing list