fix(SFRData.from_yaml(); SFRData.set_streambed_top_elevations_from_de…

…m()): refactor get_slopes() method to use general function in elevations.py module; rename SFRData.get_slopes to update_slopes. Updates slopes in from_yaml and set_streambed_top_elevations_from_dem after streambed top elevations are populated, so that slopes are consistent with streambed tops. Add tests for scripting and from_yaml contexts.

sfrmaker/elevations.py

@@ -6,6 +6,57 @@
 from sfrmaker.routing import get_nextupsegs, get_upsegs, make_graph
 
 
+def get_slopes(streambed_tops, reach_lengths, reach_numbers, outreach_numbers, 
+               default_slope=0.001, minimum_slope=0.0001,
+                maximum_slope=1.):
+    """Compute slopes by reach using values in strtop (streambed top) and rchlen (reach length)
+    columns of reach_data. The slope for a reach n is computed as strtop(n) - strtop(n+1) / rchlen(n).
+    Slopes for outlet reaches are set equal to a default value (default_slope).
+    Populates the slope column in reach_data.
+
+    Parameters
+    ----------
+    streambed_top : sequence
+        Streambed top elevations
+    reach_lengths : sequence
+        Reach lengths
+    reach_numbers : sequence
+        Unique identifiers for each reach.
+    outreach_numbers : sequence
+        Unique identifier of next downtream reach.
+    default_slope : float
+        Slope value applied to outlet reaches (where water leaves the model).
+        Default value is 0.001
+    minimum_slope : float
+        Assigned to reaches with computed slopes less than this value.
+        This ensures that the Manning's equation won't produce unreasonable values of stage
+        (in other words, that stage is consistent with assumption that
+        streamflow is primarily drive by the streambed gradient).
+        Default value is 0.0001.
+    maximum_slope : float
+        Assigned to reaches with computed slopes more than this value.
+        Default value is 1.
+    """
+    # cast everything to lists to avoid confusion with numpy vs. pandas indexers
+    streambed_tops = list(streambed_tops)
+    reach_lengths = list(reach_lengths)
+    reach_numbers = list(reach_numbers)
+    outreach_numbers = list(outreach_numbers)
+    assert np.sum(outreach_numbers) > 0, \
+        ("outreach_numbers appear to be invalid; make sure outreaches are popluated, "
+         "for example by running SFRData.set_outreaches()")
+    elev = dict(zip(reach_numbers, streambed_tops))
+    dist = dict(zip(reach_numbers, reach_lengths))
+    dnelev = {rno: elev[outreach_numbers[i]] if outreach_numbers[i] != 0
+              else -9999 for i, rno in enumerate(reach_numbers)}
+    slopes = np.array(
+        [(elev[rno] - dnelev[rno]) / dist[rno] if dnelev[rno] != -9999 and dist[rno] > 0
+            else default_slope for rno in reach_numbers])
+    slopes[slopes < minimum_slope] = minimum_slope
+    slopes[slopes > maximum_slope] = maximum_slope
+    return slopes
+
+
 def smooth_elevations(fromids, toids, elevations, start_elevations=None):  # elevup, elevdn):
     """
 

sfrmaker/sfrdata.py

@@ -12,7 +12,7 @@
 from gisutils import df2shp, get_authority_crs
 from sfrmaker.routing import find_path, renumber_segments
 from sfrmaker.checks import valid_rnos, valid_nsegs, rno_nseg_routing_consistent
-from sfrmaker.elevations import smooth_elevations
+from sfrmaker.elevations import get_slopes, smooth_elevations
 from sfrmaker.flows import add_to_perioddata, add_to_segment_data
 from sfrmaker.gis import export_reach_data, project
 from sfrmaker.observations import write_gage_package, write_mf6_sfr_obsfile, add_observations
@@ -149,6 +149,9 @@ def __init__(self, reach_data=None,
         self._period_data = None
         self._observations = None
         self._observations_filename = None
+        self._minimum_slope = minimum_slope
+        self._default_slope = default_slope
+        self._maximum_slope = maximum_slope
 
         # convert any modflow6 kwargs to modflow5
         kwargs = {SFRData.mf5names[k] if k in SFRData.mf6names else k:
@@ -175,7 +178,7 @@ def __init__(self, reach_data=None,
         # not the ideal logic for MODFLOW 6 case where only
         # rno and connections are supplied
         self.set_outreaches()
-        self.get_slopes(default_slope=default_slope, minimum_slope=minimum_slope,
+        self.update_slopes(default_slope=default_slope, minimum_slope=minimum_slope,
                         maximum_slope=maximum_slope)
 
         # have to set the model last, because it also sets up a flopy sfr package instance
@@ -1018,40 +1021,47 @@ def set_streambed_top_elevations_from_dem(self, filename, elevation_units=None,
         self.reach_data['strtop'] = [sampled_elevs[rno]
                                      for rno in self.reach_data['rno'].values]
         self.reach_data['strtop'] *= mult
+        # update the slopes as well
+        self.update_slopes()
 
-    def get_slopes(self, default_slope=0.001, minimum_slope=0.0001,
-                   maximum_slope=1.):
+    def update_slopes(self, default_slope=None, minimum_slope=None,
+                   maximum_slope=None):
         """Compute slopes by reach using values in strtop (streambed top) and rchlen (reach length)
-        columns of reach_data. The slope for a reach n is computed as strtop(n+1) - strtop(n) / rchlen(n).
+        columns of reach_data. The slope for a reach n is computed as strtop(n) - strtop(n+1) / rchlen(n).
         Slopes for outlet reaches are set equal to a default value (default_slope).
         Populates the slope column in reach_data.
 
         Parameters
         ----------
         default_slope : float
             Slope value applied to outlet reaches (where water leaves the model).
-            Default value is 0.001
+            By default None, in which case '_default_slope' attribute is used.
         minimum_slope : float
             Assigned to reaches with computed slopes less than this value.
             This ensures that the Manning's equation won't produce unreasonable values of stage
             (in other words, that stage is consistent with assumption that
             streamflow is primarily drive by the streambed gradient).
-            Default value is 0.0001.
+            By default None, in which case '_minimum_slope' attribute is used.
         maximum_slope : float
             Assigned to reaches with computed slopes more than this value.
-            Default value is 1.
+            By default None, in which case '_maximum_slope' attribute is used.
+            
+        Notes
+        -----
+        If default_slope, minimum_slope or maximum_slope are not None, the
+        respective attribute (for example '_default_slope') will be set with the supplied value.
         """
+        if default_slope is not None:
+            self._default_slope = default_slope
+        if minimum_slope is not None:
+            self._minimum_slope = minimum_slope
+        if maximum_slope is not None:
+            self._maximum_slope = maximum_slope
         rd = self.reach_data
-        assert rd.outreach.sum() > 0, "requires reach routing, must be called after set_outreaches()"
-        elev = dict(zip(rd.rno, rd.strtop))
-        dist = dict(zip(rd.rno, rd.rchlen))
-        dnelev = {rid: elev[rd.outreach.values[i]] if rd.outreach.values[i] != 0
-        else -9999 for i, rid in enumerate(rd.rno)}
-        slopes = np.array(
-            [(elev[i] - dnelev[i]) / dist[i] if dnelev[i] != -9999 and dist[i] > 0
-             else default_slope for i in rd.rno])
-        slopes[slopes < minimum_slope] = minimum_slope
-        slopes[slopes > maximum_slope] = maximum_slope
+        slopes = get_slopes(rd['strtop'], rd['rchlen'], rd['rno'], rd['outreach'], 
+                            default_slope=self._default_slope, 
+                            minimum_slope=self._minimum_slope,
+                            maximum_slope=self._maximum_slope)
         self.reach_data['slope'] = slopes
 
     @classmethod
@@ -1293,6 +1303,7 @@ def from_yaml(cls, config_file, package_name=None, output_path=None,
             sfrdata.set_streambed_top_elevations_from_dem(**dem_kwargs)
         else:
             sfrdata.reach_data['strtop'] = sfrdata.interpolate_to_reaches('elevup', 'elevdn')
+            sfrdata.update_slopes()
 
         # assign layers to the sfr reaches
         if model is not None:

sfrmaker/test/test_from_yaml.py

@@ -181,6 +181,19 @@ def test_tylerforks_from_config(config_file, dem):
     assert np.allclose(rd.strtop.values.min(), 1125, rtol=0.01)
     assert sfrdata.grid.crs.srs == 'EPSG:26715'
 
+    # check that slopes were updated after sampling streambed tops from the DEM
+    reach_data = sfrdata.reach_data.copy()
+    reach_data.index = reach_data['rno']
+    path = sfrmaker.routing.find_path(sfrdata.rno_routing, 1)[:-1]
+    strtop = reach_data.loc[path, 'strtop']
+    rchlen = reach_data.loc[path, 'rchlen']
+    slopes = -np.diff(strtop)/rchlen[:-1]
+    slopes[slopes > sfrdata._maximum_slope] = sfrdata._maximum_slope
+    slopes[slopes < sfrdata._minimum_slope] = sfrdata._minimum_slope
+    outlets = reach_data.loc[path, 'outreach'] == 0
+    slopes[outlets] = sfrdata._default_slope
+    assert np.allclose(slopes.values, reach_data.loc[path, 'slope'].values[:-1])
+
     # check that inflows were written correctly
     if 'inflows' in cfg:
         m1 = sfrdata.modflow_sfr2.parent

sfrmaker/test/test_map.py

@@ -14,7 +14,7 @@ def name_path():
 
 
 @pytest.fixture(scope="module")
-def lines(test_data_path, name_path):
+def get_lines(test_data_path, name_path):
     name, path = name_path
     lns = sfrmaker.Lines.from_shapefile(
         test_data_path / f'{path}/flowlines.shp'.format(test_data_path, path),
@@ -31,15 +31,27 @@ def lines(test_data_path, name_path):
     return lns
 
 
+@pytest.fixture(scope='function')
+def lines(get_lines):
+    return copy.deepcopy(get_lines)
+
+
 @pytest.fixture(scope="module")
-def grid(test_data_path, name_path):
+def get_grid(test_data_path, name_path):
     name, path = name_path
     grd = sfrmaker.StructuredGrid.from_json('{0}/{1}/{1}/{1}_grid.json'.format(test_data_path, path))
     return grd
 
 
+@pytest.fixture(scope='function')
+def grid(get_grid):
+    return copy.deepcopy(get_grid)
+
+
 @pytest.fixture(scope="module")
-def sfr(lines, grid):
+def get_sfr(get_lines, get_grid):
+    lines = copy.deepcopy(get_lines)
+    grid = copy.deepcopy(get_grid)
     sfr = lines.to_sfr(grid=grid,
                        isfr=None,
                        cull_flowlines_to_active_area=True,
@@ -48,8 +60,14 @@ def sfr(lines, grid):
     return sfr
 
 
+@pytest.fixture(scope='function')
+def sfr(get_sfr):
+    return copy.deepcopy(get_sfr)
+
+
 @pytest.fixture(scope="module")
-def sfr_with_inflows(sfr):
+def sfr_with_inflows(get_sfr):
+    copy.deepcopy(get_sfr)
     # add some inflows
     tmp = sfr.segment_data.loc[sfr.segment_data.nseg.isin([17, 18])].sort_values(by='nseg').copy()
     dfs = []
@@ -101,10 +119,37 @@ def test_lines_to_sfr(lines, grid):
                        default_slope=0.011, minimum_slope=0.01,
                        maximum_slope=0.012,
                        )
+    # test slope computation
+    reach_data = sfr.reach_data.copy()
+    reach_data.index = reach_data['rno']
+    path = sfrmaker.routing.find_path(sfr._routing, 1)[:-1]
+    strtop = reach_data.loc[path, 'strtop']
+    rchlen = reach_data.loc[path, 'rchlen']
+    slopes = -np.diff(strtop)/rchlen[:-1]
+
     assert np.all(sfr.reach_data.loc[sfr.reach_data.outreach == 0, \
                                      'slope'] == 0.011)
     assert sfr.reach_data.slope.min() == 0.01
     assert sfr.reach_data.slope.max() == 0.012
     # check that CRS was successfully passed to SFRData object
     assert lines.crs == sfr.crs == grid.crs
+
+
+def test_reach_data_slopes(lines, grid):
+    sfr = lines.to_sfr(grid=grid,
+                       isfr=None,
+                       cull_flowlines_to_active_area=True,
+                       one_reach_per_cell=True,
+                       minimum_slope=0.,
+                       #default_slope=0.011, minimum_slope=0.01,
+                       #maximum_slope=0.012,
+                       )
+    # test slope computation
+    reach_data = sfr.reach_data.copy()
+    reach_data.index = reach_data['rno']
+    path = sfrmaker.routing.find_path(sfr.rno_routing, 1)[:-1]
+    strtop = reach_data.loc[path, 'strtop']
+    rchlen = reach_data.loc[path, 'rchlen']
+    slopes = -np.diff(strtop)/rchlen[:-1]
+    assert np.allclose(slopes.values, reach_data.loc[path, 'slope'].values[:-1])