Update routing_product.py

.zenodo.json

@@ -24,6 +24,11 @@
  "name": "Arsenault, Richard",
  "affiliation": "École de technologie supérieure",
  "orcid": "0000-0003-2834-2750"
+ },
+ {
+ "name": "Arnal, Louise",
+ "affiliation": "Ouranos",
+ "orcid": "0000-0002-0208-2324"
  }
  ],
  "keywords": [

AUTHORS.rst

@@ -18,4 +18,4 @@ Co-Developers
 Contributors
 ------------
 
-None yet. Why not be the first?
+* Louise Arnal <arnal.louise@ouranos.ca> `@lou-a <https://github.com/lou-a>`_

CHANGELOG.rst

@@ -18,6 +18,7 @@ Internal changes
  * `ravenpy` now has a `CODE_OF_CONDUCT.md` file.
  * Many `numpydoc`-style docstrings have been adjusted for consistency.
 * Added `setuptools` to the `gis` build recipe to ensure that the `gdal` bindings are built successfully. (PR #400)
+* Modified the sub-basin and channel profile extraction functions to correctly set the river length to zero and set default values for reach attributes in sub-basins with no channel routing (i.e., sub-basins with lakes or headwater basins). (issue #354, PR #401)
 
 v0.15.0 (2024-06-20)
 --------------------

src/ravenpy/extractors/routing_product.py

@@ -148,7 +148,11 @@ def _extract_subbasin(self, row, is_lake, subbasin_ids) -> dict:
  riv_length_field = (
  "Rivlen" if self.routing_product_version == "1.0" else "RivLength"
  )
- river_length_in_kms = 0 if is_lake else round(row[riv_length_field] / 1000, 5)
+ # Correctly setting the river length to zero for sub-basins with no channel routing, such as sub-basins with lakes or headwater basins
+ if is_lake or row[riv_length_field] <= 1.0:
+ river_length_in_kms = 0
+ else:
+ river_length_in_kms = round(row[riv_length_field] / 1000, 5)
  # river_slope = max(
  # row["RivSlope"], RoutingProductShapefileExtractor.MAX_RIVER_SLOPE
  # )
@@ -192,15 +196,23 @@ def _extract_reservoir(row) -> dict:
  @staticmethod
  def _extract_channel_profile(row) -> dict:
  subbasin_id = int(row["SubId"])
- slope = max(row["RivSlope"], BasinMakerExtractor.MAX_RIVER_SLOPE)
 
- # SWAT: top width of channel when filled with water; bankfull width W_bnkfull
- channel_width = max(row["BkfWidth"], 1)
- # SWAT: depth of water in channel when filled to top of bank
- channel_depth = max(row["BkfDepth"], 1)
+ # Correctly defining the reach attributes for sub-basins with no channel routing (i.e., lakes or headwater basins)
  channel_elev = row["MeanElev"]
- floodn = row["FloodP_n"]
- channeln = row["Ch_n"]
+ if row["RivLength"] > 1.0:
+ slope = max(row["RivSlope"], BasinMakerExtractor.MAX_RIVER_SLOPE)
+ # SWAT: top width of channel when filled with water; bankfull width W_bnkfull
+ channel_width = max(row["BkfWidth"], 1)
+ # SWAT: depth of water in channel when filled to top of bank
+ channel_depth = max(row["BkfDepth"], 1)
+ floodn = row["FloodP_n"]
+ channeln = row["Ch_n"]
+ else:
+ slope = 0.12345
+ channeln = 0.12345
+ floodn = 0.12345
+ channel_width = 0.12345
+ channel_depth = 0.12345
 
  # channel profile calculations are based on theory SWAT model is based on
  # see: https://swat.tamu.edu/media/99192/swat2009-theory.pdf

tests/test_extractor.py

@@ -15,7 +15,25 @@ def test_basinmaker_extractor(get_local_testdata, tmp_path):
  routing_product_version="2.1",
  )
  rvh_config = rvh_extractor.extract(hru_from_sb=True)
+
+ # Create lists of values to check
+ bedslope_list = [item["bed_slope"] for item in rvh_config["channel_profile"]]
+ mannings_list = [
+ value
+ for d in rvh_config["channel_profile"]
+ for value in [t[1] for t in d["roughness_zones"]]
+ ]
+ reach_length_list = [item["reach_length"] for item in rvh_config["sub_basins"]]
+
  rvh_config.pop("channel_profile")
 
  config = BasicRoute(**rvh_config)
  config.write_rv(tmp_path, modelname="routing")
+
+ # Checks that the bedslope, Manning and reach length values are non negative numbers
+ assert all(isinstance(x, (int, float)) for x in bedslope_list) is True
+ assert any(x < 0 for x in bedslope_list) is False
+ assert all(isinstance(y, (int, float)) for y in mannings_list) is True
+ assert any(y < 0 for y in mannings_list) is False
+ assert all(isinstance(z, (int, float)) for z in bedslope_list) is True
+ assert any(z < 0 for z in reach_length_list) is False