fix(SFRData.create_modflow_sfr2): if no model is specified but the SF…

…R package is MODFLOW 2005-style and transient (for example, with inflows specified by stress period), add stress periods and model grid information to flopy model instance (that is created so that flopy can be used to write the SFR package).

examples/tylerforks/inflows.csv

@@ -0,0 +1,3 @@
+end_datetime,inflow_m3d,SP,line_id_COMID
+10/31/1990,15000.0000,0,1815035
+11/30/1990,16000.0000,1,1815035

examples/tylerforks/tf_sfrmaker_config2.yml

@@ -10,6 +10,12 @@ flowlines:
   PlusFlowlineVAA: NHDPlus/NHDPlusAttributes/PlusFlowlineVAA.dbf
   PlusFlow: NHDPlus/NHDPlusAttributes/PlusFlow.dbf
   elevslope: NHDPlus/NHDPlusAttributes/elevslope.dbf
+inflows:  # see sfrmaker.data.add_to_perioddata for arguments
+  filename: inflows.csv
+  #segment_column: Segment
+  line_id_column: line_id_COMID
+  period_column: SP  # column with model stress periods
+  data_column: inflow_m3d  # column with flow values
 dem:
   filename: dem_26715.tif
   elevation_units: meters

sfrmaker/sfrdata.py

@@ -643,6 +643,11 @@ def create_modflow_sfr2(self, model=None, const=None,
             model_ws = '.' if model is None else model.model_ws
             m = fm.Modflow(modelname=self.package_name, model_ws=model_ws,
                            structured=self.structured)
+            nper = self.segment_data['per'].max() + 1
+            if nper > 1:
+                dis = fm.ModflowDis(m, nrow=self.grid.nrow, ncol=self.grid.ncol, 
+                                    nlay=self.grid.nlay, nper=nper)
+
             if model is not None:
                 nper = 1
                 if 'tdis' in model.simulation.package_key_dict.keys():

sfrmaker/test/test_from_yaml.py

@@ -9,6 +9,7 @@
 import pytest
 import gisutils
 import sfrmaker
+import flopy
 from sfrmaker.fileio import read_mf6_block
 
 
@@ -179,3 +180,24 @@ def test_tylerforks_from_config(config_file, dem):
     rd = sfrdata.reach_data.dropna(subset=['asum'], axis=0)
     assert np.allclose(rd.strtop.values.min(), 1125, rtol=0.01)
     assert sfrdata.grid.crs.srs == 'EPSG:26715'
+
+    # check that inflows were written correctly
+    if 'inflows' in cfg:
+        m1 = sfrdata.modflow_sfr2.parent
+        m2 = flopy.modflow.Modflow(model_ws='.')
+        nper = sfrdata.segment_data['per'].max() + 1
+        dis = flopy.modflow.ModflowDis(m2, 
+            nrow=sfrdata.grid.nrow, ncol=sfrdata.grid.ncol, 
+            nlay=sfrdata.grid.nlay, nper=nper)
+        sfr = flopy.modflow.ModflowSfr2.load('tf.sfr', m2)
+        for i in range(nper):
+            assert np.allclose(sfr.segment_data[i]['flow'], 
+                               sfrdata.modflow_sfr2.segment_data[i]['flow'])
+        original_inflow_data = pd.read_csv(cfg['inflows']['filename'])
+        for i in range(nper):
+            per_col = cfg['inflows']['period_column']
+            q_col = cfg['inflows']['data_column']
+            q_orig = original_inflow_data.loc[original_inflow_data[per_col] == i, q_col].values[0]
+            idx = np.where(sfr.segment_data[i]['flow'] > 0)[0][0]
+            q_input = sfr.segment_data[i]['flow'][idx]
+            assert np.allclose(q_input, q_orig)