Skip to content

Commit

Permalink
v.stream.inbasin: Modernize interface, fix for v8 (#935)
Browse files Browse the repository at this point in the history
- Changes x_outlet and y_outlet to single coordinates standard option which is recognized by GUI and user can click to get coordinates there. Call in Python has less parameters and same readability (x,y).
- Makes coordinates optional which seems to be the original intention.
- Uses rules to decide about valid options instead of custom if statements.
- Updates interfacing with numpy to Python 3.
- Adds more recent reference.
- Describes stream input more.
- Adds tests.
  • Loading branch information
wenzeslaus committed Jul 29, 2023
1 parent d683353 commit bfd8226
Show file tree
Hide file tree
Showing 3 changed files with 200 additions and 38 deletions.
153 changes: 153 additions & 0 deletions src/vector/v.stream.inbasin/testsuite/test_v_stream_inbasin.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,153 @@
##############################################################################
# AUTHOR(S): Vaclav Petras <wenzeslaus gmail com>
#
# COPYRIGHT: (C) 2023 Vaclav Petras 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
# for details.
##############################################################################

"""Test v.stream.inbasin tool using grass.gunittest with NC SPM dataset"""

import json

import grass.script as gs
from grass.gunittest.case import TestCase
from grass.gunittest.main import test


class FunctionalityVStreamNetwork(TestCase):
"""Test case for v.stream.inbasin focused on functionality
The test needs v.stream.network to run.
"""

original = "streams@PERMANENT"
ours = "functionality_streams"
input_as_raster = "input_as_raster"
output_as_raster = "output_as_raster"
bad_difference = "bad_difference"
subset_output = "selected_subset"
point_output = "pour_point"
coordinates = (639089, 223304)

@classmethod
def setUpClass(cls):
cls.runModule("g.copy", vector=(cls.original, cls.ours))
cls.runModule("v.stream.network", map=cls.ours)
cls.runModule("g.region", vector=cls.ours, res=5)
cls.runModule(
"v.to.rast", input=cls.ours, output=cls.input_as_raster, use="val", value=1
)
cls.runModule(
"v.stream.inbasin",
input_streams=cls.ours,
coordinates=cls.coordinates,
output_streams=cls.subset_output,
output_pour_point=cls.point_output,
flags="s",
)
cls.runModule(
"v.to.rast",
input=cls.subset_output,
output=cls.output_as_raster,
use="val",
value=1,
)
cls.runModule(
"r.mapcalc",
expression=f"{cls.bad_difference} = if(isnull({cls.input_as_raster}) && not(isnull({cls.input_as_raster})), 1, 0)",
)

@classmethod
def tearDownClass(cls):
cls.runModule(
"g.remove",
type="vector",
name=[cls.ours, cls.subset_output, cls.point_output],
flags="f",
)
cls.runModule(
"g.remove",
type="raster",
name=[cls.input_as_raster, cls.output_as_raster, cls.bad_difference],
flags="f",
)

def test_geometry_is_subset(self):
"""Check that geometry is subset of the input by comparing raster versions"""
self.assertRasterFitsInfo(self.bad_difference, {"max": 0, "min": 0})

def test_geometry_statistics(self):
"""Check against values determined by running the tool"""
info = gs.vector_info(self.subset_output)
self.assertEqual(info["level"], 2)
self.assertEqual(info["lines"], 918)
self.assertEqual(info["nodes"], 915)
self.assertEqual(info["primitives"], 918)

def test_attribute_info(self):
"""Check attribute info against expected values"""
info = gs.vector_info(self.subset_output)
self.assertEqual(info["num_dblinks"], 1)
self.assertEqual(info["attribute_layer_number"], "1")
self.assertEqual(info["attribute_table"], self.subset_output)
self.assertEqual(info["attribute_primary_key"], "cat")

def compare_columns(self, old_vector, new_vector):
original_columns = gs.vector_columns(old_vector)
our_columns = gs.vector_columns(new_vector)
for original_column in original_columns:
self.assertIn(original_column, our_columns)
for name, original_column in original_columns.items():
for key, value in original_column.items():
self.assertEqual(value, our_columns[name][key])

def test_columns_preserved_from_original_streams(self):
"""Check that columns are preserved from original unchanged input"""
self.compare_columns(self.original, self.ours)

def test_columns_preserved_from_input(self):
"""Check that columns are preserved from input"""
self.compare_columns(self.ours, self.subset_output)

def test_point_geometry_statistics(self):
"""Check against values determined by running the tool"""
self.assertVectorExists(self.point_output)
info = gs.vector_info(self.point_output)
self.assertEqual(info["level"], 2)
self.assertEqual(info["points"], 1)
self.assertEqual(info["nodes"], 0)
self.assertEqual(info["primitives"], 1)

def test_point_attribute_info(self):
"""Check attribute info against expected values"""
self.assertVectorExists(self.point_output)
info = gs.vector_info(self.point_output)
self.assertEqual(info["num_dblinks"], 1)
self.assertEqual(info["attribute_layer_number"], "1")
self.assertEqual(info["attribute_table"], self.point_output)
self.assertEqual(info["attribute_primary_key"], "cat")

def test_pour_point_geometry(self):
"""Check pour point against coordinates determined by running the tool"""
data = (
gs.read_command("v.out.ascii", input=self.point_output).strip().split("|")
)
self.assertEqual(len(data), 3)
self.assertAlmostEqual(float(data[0]), 639258.424, places=3)
self.assertAlmostEqual(float(data[1]), 223266.931, places=3)

def test_pour_point_attributes(self):
"""Check pour point against coordinates determined by running the tool"""
records = json.loads(
gs.read_command("v.db.select", map=self.point_output, format="json")
)["records"]
self.assertEqual(len(records), 1)
self.assertAlmostEqual(records[0]["x"], 639258.424, places=3)
self.assertAlmostEqual(records[0]["y"], 223266.931, places=3)


if __name__ == "__main__":
test()
33 changes: 28 additions & 5 deletions src/vector/v.stream.inbasin/v.stream.inbasin.html
Original file line number Diff line number Diff line change
@@ -1,13 +1,36 @@
<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.
<em>v.stream.inbasin</em> uses the output of <em>v.stream.network</em> 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.

<em>v.stream.inbasin</em> expects the stream network attributes created by
<em>v.stream.network</em> and named using the names <em>v.stream.network</em>
uses by default. In other words, <em>v.stream.inbasin</em> will work on the
output <em>v.stream.network</em> with default attribute names.

<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.
<ul>
<li>
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.
</li>
<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
52 changes: 19 additions & 33 deletions src/vector/v.stream.inbasin/v.stream.inbasin.py
Original file line number Diff line number Diff line change
Expand Up @@ -4,11 +4,12 @@
# MODULE: v.stream.inbasin
#
# AUTHOR(S): Andrew Wickert
# Vaclav Petras (v8 fixes and interface improvements)
#
# PURPOSE: Build a drainage basin from the subwatersheds of a river
# network, based on the structure of the network.
#
# COPYRIGHT: (c) 2016 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 @@ -19,9 +20,6 @@
# 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
Expand Down Expand Up @@ -53,21 +51,12 @@
# % key: cat
# % label: Farthest downstream segment category
# % required: no
# % guidependency: layer,column
# %end

# %option
# % key: x_outlet
# % label: Approx. pour point x/Easting: will find closest segment
# % required: no
# % guidependency: layer,column
# %end

# %option
# % key: y_outlet
# % label: Approx. pour point y/Northing: will find closest segment
# %option G_OPT_M_COORDS
# % label: Pour point coordinates
# % description: The alogorithm will find the closest stream segment
# % required: no
# % guidependency: layer,column
# %end

# %option G_OPT_V_OUTPUT
Expand All @@ -94,22 +83,22 @@
# % guisection: Settings
# %end

# %rules
# % required: coordinates,cat
# % exclusive: coordinates,cat
# %end

##################
# IMPORT MODULES #
##################
# PYTHON
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.pygrass.modules.shortcuts import vector as v
from grass.pygrass.gis import region
from grass.pygrass import vector # Change to "v"?
from grass.pygrass import vector
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

Expand Down Expand Up @@ -137,8 +126,11 @@ def main():
streams = options["input_streams"]
basins = options["input_basins"]
downstream_cat = options["cat"]
x_outlet = float(options["x_outlet"])
y_outlet = float(options["y_outlet"])
if options["coordinates"]:
x_outlet, y_outlet = options["coordinates"].split(",")
x_outlet, y_outlet = float(x_outlet), float(y_outlet)
else:
x_outlet, y_outlet = None, None
output_basins = options["output_basin"]
output_streams = options["output_streams"]
output_pour_point = options["output_pour_point"]
Expand All @@ -148,12 +140,6 @@ def main():
# print options
# print flags

# Check that either x,y or cat are set
if (downstream_cat != "") or ((x_outlet != "") and (y_outlet != "")):
pass
else:
gscript.fatal('You must set either "cat" or "x_outlet" and "y_outlet".')

# NEED TO ADD IF-STATEMENT HERE TO AVOID AUTOMATIC OVERWRITING!!!!!!!!!!!
if snapflag or (downstream_cat != ""):
if downstream_cat == "":
Expand Down Expand Up @@ -186,11 +172,11 @@ def main():
)
# v.distance(_from_='tmp', to=streams, upload='cat', column='strcat')
downstream_cat = gscript.vector_db_select(map="tmp", columns="strcat")
downstream_cat = int(downstream_cat["values"].values()[0][0])
downstream_cat = int(downstream_cat["values"][1][0])

# Attributes of streams
colNames = np.array(vector_db_select(streams)["columns"])
colValues = np.array(vector_db_select(streams)["values"].values())
colValues = np.array(list(vector_db_select(streams)["values"].values()))
tostream = colValues[:, colNames == "tostream"].astype(int).squeeze()
cats = colValues[:, colNames == "cat"].astype(int).squeeze() # = "fromstream"

Expand Down Expand Up @@ -265,7 +251,7 @@ def main():
_pp = gscript.vector_db_select(
map=streams, columns="x2,y2", where="cat=" + str(downstream_cat)
)
_xy = np.squeeze(_pp["values"].values())
_xy = np.squeeze(list(_pp["values"].values()))
_x = float(_xy[0])
_y = float(_xy[1])
else:
Expand Down

0 comments on commit bfd8226

Please sign in to comment.