Skip to content

Commit

Permalink
v.stream.network: Fix tostream name, fix for v8 (#936)
Browse files Browse the repository at this point in the history 
- Fixes use of custom tostream which was ignored.
- Updates interfacing with numpy for Python 3.
- Updates v.to.db use for v8 where columns are created automatically (saves >10 lines here).
- Simplify code and remove old code and comments.
- Rename variables and comply more with Flake8 and Pylint.
- Explain input more.
- Add reference.
  • Loading branch information
@wenzeslaus
wenzeslaus committed Jul 29, 2023
1 parent 5633a4b commit d683353
Show file tree
Hide file tree
Showing 2 changed files with 62 additions and 118 deletions.
16 changes: 15 additions & 1 deletion src/vector/v.stream.network/v.stream.network.html
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,10 @@ <h2>DESCRIPTION</h2>
<li><em>tostream</em>: category number of the next stream; is equal to 0 if the stream flows off of the map</li>
</ul>

Any stream-like network will work, but the lines need to be connected.
In terms of graph theory, a tree and forest are supported. Behavior for
a cyclic graph is undefined.

<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):
Expand All @@ -33,7 +37,17 @@ <h2>NOTES</h2>

<h2>REFERENCES</h2>

None (yet)
<ul>
<li>
Ng, G-H. Crystal, Andrew D. Wickert, Lauren D. Somers, Leila Saberi,
Collin Cronkite-Ratcliff, Richard G. Niswonger,
and Jeffrey M. McKenzie.
"GSFLOW–GRASS v1. 0.0: GIS-enabled hydrologic modeling of coupled
groundwater–surface-water systems."
<em>Geoscientific Model Development</em> 11 (2018): 4755-4777.
<a href="https://doi.org/10.5194/gmd-11-4755-2018">DOI 10.5194/gmd-11-4755-2018</a>
</li>
</ul>

<h2>SEE ALSO</h2>

Expand Down
164 changes: 47 additions & 117 deletions src/vector/v.stream.network/v.stream.network.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,12 +4,13 @@
# MODULE: v.stream.network
#
# AUTHOR(S): Andrew Wickert
# Vaclav Petras (v8 fixes and interface improvements)
#
# 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
# COPYRIGHT: (c) 2016-2023 Andrew Wickert and the GRASS Development Team
#
# This program is free software under the GNU General Public
# License (>=v2). Read the file COPYING that comes with GRASS
Expand All @@ -20,9 +21,6 @@
# 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
Expand Down Expand Up @@ -73,150 +71,82 @@
# %option
# % key: tostream_cat_column
# % type: string
# % description: Adjacent downstream segment cat (0 = offmap flow)
# % label: Adjacent downstream segment category
# % description: Zero (0) indicates off-map flow
# % answer: tostream
# % required : no
# %end

##################
# IMPORT MODULES #
##################
# PYTHON
"""Add stream network links to attribute table"""

import numpy as np

# GRASS
from grass.pygrass.modules.shortcuts import general as g
from grass.pygrass.modules.shortcuts import raster as r
from grass import script as gs

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.pygrass.vector import VectorTopo
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 segment to its downstream segment
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()
options, unused_flags = gs.parser()
streams = options["map"]
x1 = options["upstream_easting_column"]
y1 = options["upstream_northing_column"]
x2 = options["downstream_easting_column"]
y2 = options["downstream_northing_column"]
to_stream = options["tostream_cat_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()
# Get coordinates of points: 1 = start, 2 = end
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 = []
# Read in and save the start and end coordinate points
col_names = np.array(vector_db_select(streams)["columns"])
col_values = np.array(list(vector_db_select(streams)["values"].values()))
# river number
cats = col_values[:, col_names == "cat"].astype(int).squeeze()
# upstream
xy1 = col_values[:, (col_names == "x1") + (col_names == "y1")].astype(float)
# downstream
xy2 = col_values[:, (col_names == "x2") + (col_names == "y2")].astype(float)

# Build river network
to_cats = []
for i in range(len(cats)):
tosegment_mask = np.prod(xy1 == xy2[i], axis=1)
if np.sum(tosegment_mask) == 0:
tocat.append(0)
to_segment_mask = np.prod(xy1 == xy2[i], axis=1)
if np.sum(to_segment_mask) == 0:
to_cats.append(0)
else:
tocat.append(cats[tosegment_mask.nonzero()[0][0]])
tocat = np.asarray(tocat).astype(int)
to_cats.append(cats[to_segment_mask.nonzero()[0][0]])
to_cats = np.asarray(to_cats).astype(int)

streams_vector = VectorTopo(streams)
streams_vector.open("rw")
streams_vector.table.columns.add(f"{to_stream}", "int")
streams_vector.table.conn.commit()
# 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()
cur = streams_vector.table.conn.cursor()
# Default to 0 if no stream flows to it
cur.execute("update " + streams + " set tostream=0")
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.'
cur.execute(f"update {streams} set {to_stream}=0")
for to_cat, where_cat in zip(to_cats, cats):
cur.execute(f"update {streams} set {to_stream}={to_cat} where cat={where_cat}")
streams_vector.table.conn.commit()
streams_vector.close()

gs.message(
_(
'Drainage topology built. See "{to_stream}" column for the downstream cat.'
).format(to_stream=to_stream)
)
gscript.message("A cat value of 0 indicates the downstream-most segment.")
gscript.message("")
gs.message(_("A cat value of 0 indicates the downstream-most segment."))


if __name__ == "__main__":
Expand Down

0 comments on commit d683353

Please sign in to comment.