diff --git a/sfrmaker/elevations.py b/sfrmaker/elevations.py index feae6b34..7f355ef1 100644 --- a/sfrmaker/elevations.py +++ b/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): """ diff --git a/sfrmaker/sfrdata.py b/sfrmaker/sfrdata.py index a4003601..4008906c 100644 --- a/sfrmaker/sfrdata.py +++ b/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,11 +1021,13 @@ 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. @@ -1030,28 +1035,33 @@ def get_slopes(self, default_slope=0.001, minimum_slope=0.0001, ---------- 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: diff --git a/sfrmaker/test/test_from_yaml.py b/sfrmaker/test/test_from_yaml.py index bd6af424..397d6e69 100644 --- a/sfrmaker/test/test_from_yaml.py +++ b/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 diff --git a/sfrmaker/test/test_map.py b/sfrmaker/test/test_map.py index 319a068f..580ca2ce 100644 --- a/sfrmaker/test/test_map.py +++ b/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]) \ No newline at end of file