[GRASS-SVN] r71526 - in grass-addons/grass7/vector: v.gsflow.gravres v.gsflow.grid v.gsflow.hruparams v.gsflow.segments v.stream.network
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon Oct 2 15:41:40 PDT 2017
Author: awickert
Date: 2017-10-02 15:41:40 -0700 (Mon, 02 Oct 2017)
New Revision: 71526
Modified:
grass-addons/grass7/vector/v.gsflow.gravres/v.gsflow.gravres.py
grass-addons/grass7/vector/v.gsflow.grid/v.gsflow.grid.py
grass-addons/grass7/vector/v.gsflow.hruparams/v.gsflow.hruparams.py
grass-addons/grass7/vector/v.gsflow.segments/v.gsflow.segments.py
grass-addons/grass7/vector/v.stream.network/v.stream.network.py
Log:
multiple fixes from https://github.com/awickert/GRASS-fluvial-profiler
Modified: grass-addons/grass7/vector/v.gsflow.gravres/v.gsflow.gravres.py
===================================================================
--- grass-addons/grass7/vector/v.gsflow.gravres/v.gsflow.gravres.py 2017-10-02 20:27:52 UTC (rev 71525)
+++ grass-addons/grass7/vector/v.gsflow.gravres/v.gsflow.gravres.py 2017-10-02 22:41:40 UTC (rev 71526)
@@ -107,7 +107,7 @@
"""
# Create gravity reservoirs -- overlay cells=grid and HRUs
- v.overlay(ainput=HRUs, binput=grid, operator='and', output=gravity_reservoirs, overwrite=gscript.overwrite())
+ v.overlay(ainput=HRUs, binput=grid, atype='area', btype='area', 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)
@@ -120,6 +120,8 @@
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)
+ v.extract(input=gravity_reservoirs, output='_tmp', where="gvr_cell_pct > 0.001", overwrite=True, quiet=True)
+ g.rename(vector='_tmp,'+gravity_reservoirs, overwrite=True, quiet=True)
if __name__ == "__main__":
Modified: grass-addons/grass7/vector/v.gsflow.grid/v.gsflow.grid.py
===================================================================
--- grass-addons/grass7/vector/v.gsflow.grid/v.gsflow.grid.py 2017-10-02 20:27:52 UTC (rev 71525)
+++ grass-addons/grass7/vector/v.gsflow.grid/v.gsflow.grid.py 2017-10-02 22:41:40 UTC (rev 71526)
@@ -43,13 +43,13 @@
#%option
#% key: dx
-#% label: Cell size (x / E / zonal), in map units
+#% label: Cell size suggestion (x / E / zonal), map units: rounds to DEM
#% required: yes
#%end
#%option
#% key: dy
-#% label: Cell size (y / N / meridional), in map units
+#% label: Cell size suggestion (y / N / meridional), map units: rounds to DEM
#% required: yes
#%end
@@ -104,9 +104,35 @@
mask = options['mask_output']
pp = options['pour_point']
- # Create grid
+ # Create grid -- overlaps DEM, one cell of padding
gscript.use_temp_region()
+ reg = gscript.region()
+ reg_grid_edges_sn = np.linspace(reg['s'], reg['n'], reg['rows'])
+ reg_grid_edges_we = np.linspace(reg['w'], reg['e'], reg['cols'])
g.region(vector=basin, ewres=dx, nsres=dy)
+ regnew = gscript.region()
+ # Use a grid ratio -- don't match exactly the desired MODFLOW resolution
+ grid_ratio_ns = np.round(regnew['nsres']/reg['nsres'])
+ grid_ratio_ew = np.round(regnew['ewres']/reg['ewres'])
+ # Get S, W, and then move the unit number of grid cells over to get N and E
+ # and include 1 (new) cell of padding around the whole watershed
+ _s_dist = np.abs(reg_grid_edges_sn - (regnew['s'] - 2.*regnew['nsres']) )
+ _s_idx = np.where(_s_dist == np.min(_s_dist))[0][0]
+ _s = float(reg_grid_edges_sn[_s_idx])
+ _n_grid = np.arange(_s, reg['n'] + 2*grid_ratio_ns*reg['nsres'], grid_ratio_ns*reg['nsres'])
+ _n_dist = np.abs(_n_grid - (regnew['n'] + 2.*regnew['nsres']))
+ _n_idx = np.where(_n_dist == np.min(_n_dist))[0][0]
+ _n = float(_n_grid[_n_idx])
+ _w_dist = np.abs(reg_grid_edges_we - (regnew['w'] - 2.*regnew['ewres']))
+ _w_idx = np.where(_w_dist == np.min(_w_dist))[0][0]
+ _w = float(reg_grid_edges_we[_w_idx])
+ _e_grid = np.arange(_w, reg['e'] + 2*grid_ratio_ew*reg['ewres'], grid_ratio_ew*reg['ewres'])
+ _e_dist = np.abs(_e_grid - (regnew['e'] + 2.*regnew['ewres']))
+ _e_idx = np.where(_e_dist == np.min(_e_dist))[0][0]
+ _e = float(_e_grid[_e_idx])
+ # Finally make the region
+ g.region(w=str(_w), e=str(_e), s=str(_s), n=str(_n), nsres=str(grid_ratio_ns*reg['nsres']), ewres=str(grid_ratio_ew*reg['ewres']))
+ # And then make the grid
v.mkgrid(map=grid, overwrite=gscript.overwrite())
# Cell numbers (row, column, continuous ID)
@@ -145,6 +171,8 @@
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)
+ g.region(n=reg['n'], s=reg['s'], w=reg['w'], e=reg['e'], nsres=reg['nsres'], ewres=reg['ewres'])
+
if __name__ == "__main__":
main()
Modified: grass-addons/grass7/vector/v.gsflow.hruparams/v.gsflow.hruparams.py
===================================================================
--- grass-addons/grass7/vector/v.gsflow.hruparams/v.gsflow.hruparams.py 2017-10-02 20:27:52 UTC (rev 71525)
+++ grass-addons/grass7/vector/v.gsflow.hruparams/v.gsflow.hruparams.py 2017-10-02 22:41:40 UTC (rev 71526)
@@ -81,14 +81,26 @@
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.vector.geometry import Point
from grass.pygrass.raster import RasterRow
from grass.pygrass import utils
from grass import script as gscript
-###############
-# MAIN MODULE #
-###############
+################
+# MAIN MODULES #
+################
+def create_iterator(vect):
+ colNames = np.array(gscript.vector_db_select(vect, layer=1)['columns'])
+ colValues = np.array(gscript.vector_db_select(vect, layer=1)['values'].values())
+ cats = colValues[:,colNames == 'cat'].astype(int).squeeze()
+ _n = np.arange(1, len(cats)+1)
+ _n_cats = []
+ for i in range(len(cats)):
+ _n_cats.append( (_n[i], cats[i]) )
+ return _n_cats
+
+
def main():
"""
Adds GSFLOW parameters to a set of HRU sub-basins
@@ -125,6 +137,8 @@
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_lon double precision') # Longitude of centroid
+ # unnecessary but why not?
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
@@ -249,7 +263,7 @@
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!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ # CENTROIDS
############
# get x,y of centroid -- but have areas not in database table, that do have
@@ -260,51 +274,99 @@
# 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 = VectorTopo(HRU)
+ hru.open('rw')
+ hru_cats = []
+ hru_coords = []
+ for hru_i in hru:
+ if type(hru_i) is vector.geometry.Centroid:
+ hru_cats.append(hru_i.cat)
+ hru_coords.append(hru_i.coords())
+ hru_cats = np.array(hru_cats)
+ hru_coords = np.array(hru_coords)
+ hru.rewind()
+
+ hru_area_ids = []
+ for coor in hru_coords:
+ _area = hru.find_by_point.area(Point(coor[0], coor[1]))
+ hru_area_ids.append(_area)
+ hru_area_ids = np.array(hru_area_ids)
+ hru.rewind()
- # 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)
+ hru_areas = []
+ for _area_id in hru_area_ids:
+ hru_areas.append(_area_id.area())
+ hru_areas = np.array(hru_areas)
+ hru.rewind()
+
+ allcats = sorted(list(set(list(hru_cats))))
+
+ # Now create weighted mean
+ hru_centroid_locations = []
+ for cat in allcats:
+ hrus_with_cat = hru_cats[hru_cats == cat]
+ if len(hrus_with_cat) == 1:
+ hru_centroid_locations.append((hru_coords[hru_cats == cat]).squeeze())
+ else:
+ _centroids = hru_coords[hru_cats == cat]
+ #print _centroids
+ _areas = hru_areas[hru_cats == cat]
+ #print _areas
+ _x = np.average(_centroids[:,0], weights=_areas)
+ _y = np.average(_centroids[:,1], weights=_areas)
+ #print _x, _y
+ hru_centroid_locations.append(np.array([_x, _y]))
+
+ # Now upload weighted mean to database table
+ # allcats and hru_centroid_locations are co-indexed
+ index__cats = create_iterator(HRU)
+ cur = hru.table.conn.cursor()
+ for i in range(len(allcats)):
+ # meters
+ cur.execute('update '+HRU
+ +' set hru_x='+str(hru_centroid_locations[i][0])
+ +' where cat='+str(allcats[i]))
+ cur.execute('update '+HRU
+ +' set hru_y='+str(hru_centroid_locations[i][1])
+ +' where cat='+str(allcats[i]))
+ # feet
+ cur.execute('update '+HRU
+ +' set hru_xlong='+str(hru_centroid_locations[i][0]*3.28084)
+ +' where cat='+str(allcats[i]))
+ cur.execute('update '+HRU
+ +' set hru_ylat='+str(hru_centroid_locations[i][1]*3.28084)
+ +' where cat='+str(allcats[i]))
+ # (un)Project to lat/lon
+ _centroid_ll = gscript.parse_command('m.proj',
+ coordinates=
+ list(hru_centroid_locations[i]),
+ flags='od').keys()[0]
+ _lon, _lat, _z = _centroid_ll.split('|')
+ cur.execute('update '+HRU
+ +' set hru_lon='+_lon
+ +' where cat='+str(allcats[i]))
+ cur.execute('update '+HRU
+ +' set hru_lat='+_lat
+ +' where cat='+str(allcats[i]))
- # why not just get lon too?
- v.db_addcolumn(map=HRU, columns='hru_lon double precision')
+ # feet -- not working.
+ # Probably an issue with index__cats -- maybe fix later, if needed
+ # But currently not a major speed issue
+ """
+ cur.executemany("update "+HRU+" set hru_xlong=?*3.28084 where hru_x=?",
+ index__cats)
+ cur.executemany("update "+HRU+" set hru_ylat=?*3.28084 where hru_y=?",
+ index__cats)
+ """
- hru = VectorTopo(HRU)
- hru.open('rw')
- cur = hru.table.conn.cursor()
- cur.executemany("update "+ HRU +" set hru_lon=?, hru_lat=? where cat=?", lonlat_cat)
+ cur.close()
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()
- """
+ #cur.executemany("update "+HRU+" set hru_segment=? where id=?",
+ # index__cats)
# Segment number = HRU ID number
v.db_update(map=HRU, column='hru_segment', query_column='id', quiet=True)
Modified: grass-addons/grass7/vector/v.gsflow.segments/v.gsflow.segments.py
===================================================================
--- grass-addons/grass7/vector/v.gsflow.segments/v.gsflow.segments.py 2017-10-02 20:27:52 UTC (rev 71525)
+++ grass-addons/grass7/vector/v.gsflow.segments/v.gsflow.segments.py 2017-10-02 22:41:40 UTC (rev 71526)
@@ -54,8 +54,8 @@
#%option
#% key: cdpth
#% type: double
-#% description: Flow depth coefficient; used if ICALC=3
-#% answer: 0.25
+#% description: Flow depth coefficient [meters]; used if ICALC=3
+#% answer: 1
#% required: no
#%end
@@ -63,15 +63,15 @@
#% key: fdpth
#% type: double
#% description: Flow depth exponent; used if ICALC=3
-#% answer: 0.4
+#% answer: 0.42
#% required: no
#%end
#%option
#% key: awdth
#% type: double
-#% description: Flow width coefficient; used if ICALC=3
-#% answer: 8
+#% description: Flow width coefficient [meters]; used if ICALC=3
+#% answer: 20
#% required: no
#%end
@@ -79,7 +79,7 @@
#% key: bwdth
#% type: double
#% description: Flow width exponent; used if ICALC=3
-#% answer: 0.5
+#% answer: 0.23
#% required: no
#%end
@@ -209,9 +209,10 @@
ROUGHBK = options['roughbk']
# ICALC=3: Power-law relationships (following Leopold and others)
- CDPTH = options['cdpth']
+ # The at-a-station default exponents are from Rhodes (1977)
+ CDPTH = str(float(options['cdpth']) / 35.3146667) # cfs to m^3/s
FDPTH = options['fdpth']
- AWDTH = options['awdth']
+ AWDTH = str(float(options['awdth']) / 35.3146667) # cfs to m^3/s
BWDTH = options['bwdth']
##################################################
Modified: grass-addons/grass7/vector/v.stream.network/v.stream.network.py
===================================================================
--- grass-addons/grass7/vector/v.stream.network/v.stream.network.py 2017-10-02 20:27:52 UTC (rev 71525)
+++ grass-addons/grass7/vector/v.stream.network/v.stream.network.py 2017-10-02 22:41:40 UTC (rev 71526)
@@ -194,8 +194,6 @@
cur = streamsTopo.table.conn.cursor()
# Default to 0 if no stream flows to it
cur.execute("update "+streams+" set tostream=0")
- # 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()
More information about the grass-commit
mailing list