fix(mover.py::get_sfr_package_connections): refactor to use GWFGWF ex…

…changes and a threshold distance to determine SFR Package connections between a parent and inset/child LGR model. This fixes previous issue where SFR streams meandering in and out of the inset model were not getting connected to the parent model (due to SFRmaker connecting the segments within either model).

mfsetup/fileio.py

@@ -265,7 +265,8 @@ def set_cfg_paths_to_absolute(cfg, config_file_location):
                     if input_block == 'lgr':
                         for model_name, config in pckg[input_block].items():
                             if 'filename' in config:
-                                file_keys = _parse_file_path_keys_from_source_data(config)
+                                file_keys = _parse_file_path_keys_from_source_data(
+                                    {model_name: config})
                     else:
                         file_keys = _parse_file_path_keys_from_source_data(pckg[input_block])
                     for key in file_keys:

mfsetup/mf6model.py

@@ -939,6 +939,25 @@ def setup_simulation_mover(self):
                     if inset.get_package('sfr'):
                         inset_perioddata = get_mover_sfr_package_input(self, inset)
                         perioddata_dfs.append(inset_perioddata)
+                        # for each SFR reach with a connection
+                        # to a reach in another model
+                        # set the SFR Package downstream connection to 0
+                        for i, r in inset_perioddata.iterrows():
+                            rd = self.simulation.get_model(r['mname1']).sfrdata.reach_data
+                            rd.loc[rd['rno'] == r['id1']+1, 'outreach'] = 0
+                            # fix flopy connectiondata as well
+                            sfr_package = self.simulation.get_model(r['mname1']).sfr
+                            cd = sfr_package.connectiondata.array.tolist()
+                            # there should be no downstream reaches
+                            # (indicated by negative numbers)
+                            cd[r['id1']] = tuple(v for v in cd[r['id1']] if v > 0)
+                            sfr_package.connectiondata = cd
+                        # re-write the shapefile exports with corrected routing
+                        inset.sfrdata.write_shapefiles(f'{inset._shapefiles_path}/{inset_name}')
+
+                self.sfrdata.write_shapefiles(f'{self._shapefiles_path}/{self.name}')
+
+
         if len(perioddata_dfs) > 0:
             perioddata = pd.concat(perioddata_dfs)
             if len(perioddata) > 0:

mfsetup/mfmodel.py

@@ -1836,8 +1836,8 @@ def setup_from_cfg(cls, cfg, verbose=False):
             for k, v in m.inset.items():
                 if v._is_lgr:
                     v.setup_packages()
-            m.setup_simulation_mover()
             m.setup_lgr_exchanges()
+            m.setup_simulation_mover()
 
         print('finished setting up model in {:.2f}s'.format(time.time() - t0))
         print('\n{}'.format(m))

mfsetup/mover.py

@@ -4,6 +4,8 @@
 import numpy as np
 import pandas as pd
 from scipy.spatial.distance import cdist
+from sfrmaker.routing import find_path
+from shapely.geometry import Point
 
 
 def get_connections(from_features, to_features, distance_threshold=250):
@@ -36,19 +38,22 @@ def get_connections(from_features, to_features, distance_threshold=250):
     return connections
 
 
-def get_sfr_package_connections(reach_data1, reach_data2, distance_threshold=1000):
-    """Connect SFR reaches between two packages (for example, in a parent and inset model).
-    Connections are made when a headwater reach in one package is within distance_threshold
-    of an outlet in the other package.
+def get_sfr_package_connections(gwfgwf_exchangedata,
+                                reach_data1, reach_data2, distance_threshold=1000):
+    """Connect SFR reaches between SFR packages in a pair of groundwater flow models linked
+    by the GWFGWF (simulation-level mover) Package.
 
     Parameters
     ----------
+    gwfgwf_exchangedata : flopy recarray or pandas DataFrame
+        Exchange data from the GWFGWF package
+        (listing cell connections between two groundwater flow models).
     reach_data1 : DataFrame, similar to sfrmaker.SFRData.reach_data
-        Reach information for first package to connect.
+        SFR reach information for model 1 in gwfgwf_exchangedata.
         Must contain reach numbers and 'geometry' column of shapely geometries
         for each reach (can be LineStrings or Polygons)
     reach_data2 : DataFrame, similar to sfrmaker.SFRData.reach_data
-        Reach information for second package to connect.
+        SFR reach information for model 2 in gwfgwf_exchangedata.
         Must contain reach numbers and 'geometry' column of shapely geometries
         for each reach (can be LineStrings or Polygons)
     distance_threshold : float
@@ -60,34 +65,104 @@ def get_sfr_package_connections(reach_data1, reach_data2, distance_threshold=100
     connections1 : dictionary of connections from package 1 to package 2
     connections2 : dictionary of connections from package 2 to package 1
     """
-    outlets1 = reach_data1.loc[reach_data1.outreach == 0]
-    headwaters1 = set(reach_data1.rno).difference(reach_data1.outreach)
-    headwaters1 = reach_data1.loc[reach_data1.rno.isin(headwaters1)]
-    outlets2 = reach_data2.loc[reach_data2.outreach == 0]
-    headwaters2 = set(reach_data2.rno).difference(reach_data2.outreach)
-    headwaters2 = reach_data2.loc[reach_data2.rno.isin(headwaters2)]
-
-    # get the connections between outlets and headwaters that are less than distance_threshold
-    connections1_idx = get_connections(outlets1.geometry, headwaters2.geometry,
-                                       distance_threshold=distance_threshold)
-    connections2_idx = get_connections(outlets2.geometry, headwaters1.geometry,
-                                       distance_threshold=distance_threshold)
-    # map those back to actual reach numbers
-    connections1 = {outlets1.rno.values[k]: headwaters2.rno.values[v]
-                    for k, v in connections1_idx.items()}
-    connections2 = {outlets2.rno.values[k]: headwaters1.rno.values[v]
-                    for k, v in connections2_idx.items()}
-
-    return connections1, connections2
+    gwfgwf_exchangedata = pd.DataFrame(gwfgwf_exchangedata)
+    all_reach_data1 = reach_data1.copy()
+    all_reach_data2 = reach_data2.copy()
+
+    # add cellids to connect reaches with GWFGWF exchanges
+    # ignore layers in case there are any mismatches
+    # due to pinchouts at the different model resolutions
+    # or different layer assignments by SFRmaker
+    # and because we only expect one SFR reach at each i, j location
+    gwfgwf_exchangedata['cellidm1_2d'] =\
+        [(i, j) for k, i, j in gwfgwf_exchangedata['cellidm1']]
+    gwfgwf_exchangedata['cellidm2_2d'] =\
+        [(i, j) for k, i, j in gwfgwf_exchangedata['cellidm2']]
+    for df in all_reach_data1, all_reach_data2:
+        if 'cellid' not in df.columns and\
+            {'k', 'i', 'j'}.intersection(df.columns):
+                df['cellid_2d'] =\
+                    list(df[['i', 'j']].itertuples(index=False, name=None))
+        else:
+            df['cellid_2d'] = [(i, j) for k, i, j in df['cellid']]
+    reach_data1 = all_reach_data1.loc[all_reach_data1['cellid_2d'].isin(
+        gwfgwf_exchangedata['cellidm1_2d'])]
+    reach_data2 = all_reach_data2.loc[all_reach_data2['cellid_2d'].isin(
+        gwfgwf_exchangedata['cellidm2_2d'])]
+
+    # starting locations of each reach along the parent/inset interface
+    reach_data1['reach_start'] = [Point(g.coords[0]) for g in reach_data1['geometry']]
+    reach_data2['reach_start'] = [Point(g.coords[0]) for g in reach_data2['geometry']]
+
+    # neighboring cells (in a structured grid)
+    def neighbors(i, j):
+        return {(i-1, j-1), (i-1, j), (i-1, j+1),
+                (i, j+1), (i+1, j+1), (i+1, j),
+                (i+1, j-1), (i, j-1)}
+
+    # make the connections
+    parent_to_inset = dict()
+    inset_to_parent = dict()
+    for _, r in reach_data1.iterrows():
+        # check for connections to the other model
+        # if the next reach is in another cell
+        # along the parent/inset model interface,
+        # or if the curret reach is an outlet
+        if r['outreach'] in reach_data1['rno'].values or r['outreach'] == 0:
+            # if the outreach is in a neighboring cell,
+            # assume it is correct
+            # (that the next reach is not in the other model)
+            out_reach = reach_data1.loc[reach_data1['rno'] == r['outreach']]
+            if r['outreach'] != 0 and\
+                out_reach['cellid_2d'].iloc[0] in neighbors(*r['cellid_2d']):
+                continue
+            # if the distance to the next reach is greater than 1 cell
+            # consider this each to be an outlet
+            reach_end = Point(r['geometry'].coords[-1])
+            distances = reach_end.distance(reach_data2['reach_start'])
+            if np.min(distances) < distance_threshold:
+                next_reach = reach_data2.iloc[np.argmin(distances)]
+                parent_to_inset[r['rno']] = next_reach['rno']
+        # next reach is somewhere else in the parent model
+        else:
+            continue
+    # repeat for inset to parent connections
+    for i, r in reach_data2.iterrows():
+        # check for connections to the other model
+        # if the next reach is in another cell
+        # along the parent/inset model interface,
+        # or if the curret reach is an outlet
+        if r['outreach'] in reach_data2['rno'].values or r['outreach'] == 0:
+            # if the outreach is in a neighboring cell,
+            # assume it is correct
+            # (that the next reach is not in the other model)
+            out_reach = reach_data2.loc[reach_data2['rno'] == r['outreach']]
+            if r['outreach'] != 0 and\
+                out_reach['cellid_2d'].iloc[0] in neighbors(*r['cellid_2d']):
+                continue
+            # if the distance to the next reach is greater than 1 cell
+            # consider this each to be an outlet
+            reach_end = Point(r['geometry'].coords[-1])
+            distances = reach_end.distance(reach_data1['reach_start'])
+            if np.min(distances) < distance_threshold:
+                next_reach = reach_data1.iloc[np.argmin(distances)]
+                inset_to_parent[r['rno']] = next_reach['rno']
+        # next reach is somewhere else in this model
+        else:
+            continue
+    return parent_to_inset, inset_to_parent
 
 
 def get_mover_sfr_package_input(parent, inset, convert_to_zero_based=True):
 
     grid_spacing = parent.dis.delc.array[0]
     connections = []
-    to_inset, to_parent = get_sfr_package_connections(parent.sfrdata.reach_data,
+    # use 2x grid spacing for distance threshold
+    # because reaches could be small fragments in opposite corners of two adjacent cells
+    to_inset, to_parent = get_sfr_package_connections(parent.simulation.gwfgwf.exchangedata.array,
+                                                      parent.sfrdata.reach_data,
                                                       inset.sfrdata.reach_data,
-                                                      distance_threshold=grid_spacing)
+                                                      distance_threshold=2*grid_spacing)
     # convert to zero-based if reach_data aren't
     # corrections are quantities to subtract off
     parent_rno_correction = parent.sfrdata.reach_data.rno.min()

mfsetup/tests/test_lgr.py

@@ -230,7 +230,7 @@ def test_lgr_grid_setup(lgr_grid_spacing, layer_refinement_input,
         ax.set_ylim(b, t)
         ax.set_xlim(l, r)
         ax.set_aspect(1)
-        plt.savefig(project_root_path / '/docs/source/_static/pleasant_lgr.png',
+        plt.savefig(project_root_path / 'docs/source/_static/pleasant_lgr.png',
                     bbox_inches='tight')
 
 
@@ -300,8 +300,9 @@ def test_mover_get_sfr_package_connections(pleasant_lgr_setup_from_yaml):
 
     parent_reach_data = m.sfrdata.reach_data
     inset_reach_data = inset_model.sfrdata.reach_data
-    to_inset, to_parent = get_sfr_package_connections(parent_reach_data, inset_reach_data,
-                                                      distance_threshold=200)
+    to_inset, to_parent = get_sfr_package_connections(
+        m.simulation.gwfgwf.exchangedata.array,
+        parent_reach_data, inset_reach_data, distance_threshold=200)
     assert len(to_inset) == 0
     # verify that the last reaches in the two segments are keys
     last_reaches = inset_model.sfrdata.reach_data.groupby('iseg').last().rno
@@ -326,7 +327,12 @@ def test_meandering_sfr_connections(shellmound_cfg, project_root_path, tmpdir):
     between parent and child models."""
     cfg = deepcopy(shellmound_cfg)
     default_cfg = load(project_root_path / 'mfsetup/mf6_defaults.yml')
+    del cfg['dis']['source_data']['idomain']
+    cfg['dis']['griddata']['idomain'] = 1
     lgr_inset_cfg = deepcopy(default_cfg)
+    del cfg['sfr']['sfrmaker_options']['to_riv']
+    inset_sfr_config = deepcopy(cfg['sfr'])
+    del inset_sfr_config['sfrmaker_options']['add_outlets']
     specified_lgr_inset_cfg = {
         'simulation': {
         'sim_name': 'shellmound',
@@ -337,7 +343,7 @@ def test_meandering_sfr_connections(shellmound_cfg, project_root_path, tmpdir):
         'model': {
         'simulation': 'shellmound',
         'modelname': 'shellmound_lgr',
-        'packages': ['dis'],
+        'packages': ['dis', 'sfr'],
         'list_filename_fmt': '{}.list'
         },
         'setup_grid': {
@@ -354,11 +360,18 @@ def test_meandering_sfr_connections(shellmound_cfg, project_root_path, tmpdir):
                 'nrow': 4,
                 'ncol': 16
             },
-            # 'griddata':
-        }
+            'griddata': {
+                'delr': 500,
+                'delc': 500,
+                'idomain': 1
+            },
+            'source_data': deepcopy(cfg['dis']['source_data'])
+
+        },
+        'sfr': inset_sfr_config
+
     }
     update(lgr_inset_cfg, specified_lgr_inset_cfg)
-
     cfg['simulation']['sim_ws'] = tmpdir / 'shellmound_lgr'
     cfg['setup_grid']['lgr'] = {
         'shellmound_lgr': {
@@ -374,8 +387,38 @@ def test_meandering_sfr_connections(shellmound_cfg, project_root_path, tmpdir):
         if v._is_lgr:
             v.setup_dis()
             v.setup_sfr()
-    m.setup_simulation_mover()
-    j=2
+    gwfgwf = m.setup_lgr_exchanges()
+    mvr = m.setup_simulation_mover()
+
+    # just test the connections for period 0
+    exchangedata = pd.DataFrame(m.simulation.mvr.perioddata.array[0])
+    # the reach numbers for these are liable to change
+    # check SFR package shapefile output in test output folder
+    # to verify that these are correct
+    expected_connections = {
+    ('shellmound', 173, 'shellmound_lgr', 7),
+    ('shellmound', 235, 'shellmound_lgr', 0),
+    ('shellmound', 293, 'shellmound_lgr', 5),
+    ('shellmound_lgr', 4, 'shellmound', 258),
+    ('shellmound_lgr', 14, 'shellmound', 186),
+    ('shellmound_lgr', 17, 'shellmound', 170)
+    }
+    assert set(exchangedata[['mname1', 'id1', 'mname2', 'id2']].
+               itertuples(index=False, name=None)) == expected_connections
+    cd1 = pd.DataFrame(m.simulation.get_model('shellmound').sfr.connectiondata.array)
+    cd1.index = cd1['ifno']
+    # None of the reaches should have downstream connections
+    # within their SFR Package
+    # (no negative numbers indicates an outlet condition)
+    assert not any(cd1.loc[[173, 235, 293]].iloc[:, 1] < 0)
+    rd1 = m.sfrdata.reach_data
+    assert all(rd1.loc[rd1['rno'].isin(np.array([173, 235, 293])+1), 'outreach'] == 0)
+    inset = m.simulation.get_model('shellmound_lgr')
+    cd2 = pd.DataFrame(inset.sfr.connectiondata.array)
+    cd2.index = cd2['ifno']
+    assert not any(cd2.loc[[4, 14, 17]].iloc[:, 1] < 0)
+    rd2 = inset.sfrdata.reach_data
+    assert all(rd2.loc[rd2['rno'].isin(np.array([4, 14, 17])+1), 'outreach'] == 0)
 
 
 def test_lgr_bottom_elevations(pleasant_vertical_lgr_setup_from_yaml, mf6_exe):