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 <[email protected]> `@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