From 97f4c97ad05381e15155fa8e444604f640600787 Mon Sep 17 00:00:00 2001 From: Vaclav Petras Date: Wed, 26 Jul 2023 21:04:24 -0400 Subject: [PATCH] v.stream.inbasin: Modernize interface, fix for v8 - 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. --- .../testsuite/test_v_stream_inbasin.py | 153 ++++++++++++++++++ .../v.stream.inbasin/v.stream.inbasin.html | 33 +++- .../v.stream.inbasin/v.stream.inbasin.py | 52 +++--- 3 files changed, 200 insertions(+), 38 deletions(-) create mode 100644 src/vector/v.stream.inbasin/testsuite/test_v_stream_inbasin.py diff --git a/src/vector/v.stream.inbasin/testsuite/test_v_stream_inbasin.py b/src/vector/v.stream.inbasin/testsuite/test_v_stream_inbasin.py new file mode 100644 index 0000000000..5b4f71af14 --- /dev/null +++ b/src/vector/v.stream.inbasin/testsuite/test_v_stream_inbasin.py @@ -0,0 +1,153 @@ +############################################################################## +# AUTHOR(S): Vaclav Petras +# +# 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() diff --git a/src/vector/v.stream.inbasin/v.stream.inbasin.html b/src/vector/v.stream.inbasin/v.stream.inbasin.html index 88e6b372b5..7728a5fe03 100644 --- a/src/vector/v.stream.inbasin/v.stream.inbasin.html +++ b/src/vector/v.stream.inbasin/v.stream.inbasin.html @@ -1,13 +1,36 @@

DESCRIPTION

-v.stream.inbasin 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. +v.stream.inbasin 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. + +v.stream.inbasin expects the stream network attributes created by +v.stream.network and named using the names v.stream.network +uses by default. In other words, v.stream.inbasin will work on the +output v.stream.network with default attribute names.

REFERENCES

-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. +

SEE ALSO

diff --git a/src/vector/v.stream.inbasin/v.stream.inbasin.py b/src/vector/v.stream.inbasin/v.stream.inbasin.py index bfacf9ac0d..8fc9864055 100644 --- a/src/vector/v.stream.inbasin/v.stream.inbasin.py +++ b/src/vector/v.stream.inbasin/v.stream.inbasin.py @@ -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 @@ -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 @@ -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 @@ -94,6 +83,11 @@ # % guisection: Settings # %end +# %rules +# % required: coordinates,cat +# % exclusive: coordinates,cat +# %end + ################## # IMPORT MODULES # ################## @@ -101,15 +95,10 @@ 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 @@ -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"] @@ -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 == "": @@ -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" @@ -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: