From 142acb6cacb8c57753769dd6017b9287293092dc Mon Sep 17 00:00:00 2001 From: Christian Langevin Date: Tue, 22 Oct 2019 11:56:05 -0500 Subject: [PATCH 1/4] fix(ulstrd): package lists for wel, drn, ghb, etc. were not always read correctly * flopy will now read packages that have lists specified with EXTERNAL * flopy will now read packages that have lists that use the SFAC keyword * flopy now has a general list reader flopy.utils.flopy_io.ulstrd patterned after the MODFLOW ULSTRD subroutine * flopy will now use this ulstrd reader to read lists defined with parameters * With the m.write_input(), binary external files always had REPLACE tacked onto the end. This was corrected so that REPLACE is only added when output is set to True. * Closes #683 --- autotest/t003_test.py | 16 ++--- autotest/t067_ulstrd.py | 137 +++++++++++++++++++++++++++++++++++++++ flopy/modflow/mf.py | 10 ++- flopy/modflow/mfchd.py | 4 ++ flopy/modflow/mfdrn.py | 4 ++ flopy/modflow/mfghb.py | 4 ++ flopy/modflow/mfhfb.py | 8 ++- flopy/modflow/mfparbc.py | 23 ++----- flopy/modflow/mfriv.py | 4 ++ flopy/modflow/mfwel.py | 4 ++ flopy/pakbase.py | 83 +++++------------------- flopy/utils/flopy_io.py | 130 +++++++++++++++++++++++++++++++++++++ 12 files changed, 330 insertions(+), 97 deletions(-) create mode 100644 autotest/t067_ulstrd.py diff --git a/autotest/t003_test.py b/autotest/t003_test.py index 70b4ed9928..90c1c4031b 100644 --- a/autotest/t003_test.py +++ b/autotest/t003_test.py @@ -108,11 +108,11 @@ def test_load_nam_mt_nonexistant_file(): if __name__ == '__main__': - test_loadoc() - test_loadoc_nstpfail() - test_load_nam_mf_nonexistant_file() - test_load_nam_mt_nonexistant_file() - # test_loadoc_lenfail() - # test_loadfreyberg() - # test_loadoahu() - # test_loadtwrip() + #test_loadoc() + #test_loadoc_nstpfail() + #test_load_nam_mf_nonexistant_file() + #test_load_nam_mt_nonexistant_file() + #test_loadoc_lenfail() + #test_loadfreyberg() + test_loadoahu() + #test_loadtwrip() diff --git a/autotest/t067_ulstrd.py b/autotest/t067_ulstrd.py new file mode 100644 index 0000000000..7e771ea868 --- /dev/null +++ b/autotest/t067_ulstrd.py @@ -0,0 +1,137 @@ +import os +import numpy as np +import flopy + +tpth = os.path.join('temp', 't067') +if not os.path.isdir(tpth): + os.makedirs(tpth) + + +def test_ulstrd(): + + # Create an original model and then manually modify to use + # advanced list reader capabilities + ws = tpth + nlay = 1 + nrow = 10 + ncol = 10 + nper = 3 + + # create the ghbs + ghbra = flopy.modflow.ModflowGhb.get_empty(20) + l = 0 + for i in range(nrow): + ghbra[l] = (0, i, 0, 1., 100. + i) + l += 1 + ghbra[l] = (0, i, ncol - 1, 1., 200. + i) + l += 1 + ghbspd = {0: ghbra} + + # create the drains + drnra = flopy.modflow.ModflowDrn.get_empty(2) + drnra[0] = (0, 1, int(ncol / 2), .5, 55.) + drnra[1] = (0, 2, int(ncol / 2), .5, 75.) + drnspd = {0: drnra} + + # create the wells + welra = flopy.modflow.ModflowWel.get_empty(2) + welra[0] = (0, 1, 1, -5.) + welra[1] = (0, nrow - 3, ncol - 3, -10.) + welspd = {0: welra} + + m = flopy.modflow.Modflow(modelname='original', model_ws=ws, + exe_name='mf2005') + dis = flopy.modflow.ModflowDis(m, nlay=nlay, nrow=nrow, ncol=ncol, + nper=nper) + bas = flopy.modflow.ModflowBas(m) + lpf = flopy.modflow.ModflowLpf(m) + ghb = flopy.modflow.ModflowGhb(m, stress_period_data=ghbspd) + drn = flopy.modflow.ModflowDrn(m, stress_period_data=drnspd) + wel = flopy.modflow.ModflowWel(m, stress_period_data=welspd) + pcg = flopy.modflow.ModflowPcg(m) + oc = flopy.modflow.ModflowOc(m) + m.add_external('original.drn.dat', 71) + m.add_external('original.wel.bin', 72, binflag=True, output=False) + m.write_input() + + # rewrite ghb + fname = os.path.join(ws, 'original.ghb') + with open(fname, 'w') as f: + f.write('{} {}\n'.format(ghbra.shape[0], 0)) + for kper in range(nper): + f.write('{} {}\n'.format(ghbra.shape[0], 0)) + f.write('open/close original.ghb.dat\n') + + # write ghb list + sfacghb = 5 + fname = os.path.join(ws, 'original.ghb.dat') + with open(fname, 'w') as f: + f.write('sfac {}\n'.format(sfacghb)) + for k, i, j, stage, cond in ghbra: + f.write('{} {} {} {} {}\n'.format(k + 1, i + 1, j + 1, stage, cond)) + + # rewrite drn + fname = os.path.join(ws, 'original.drn') + with open(fname, 'w') as f: + f.write('{} {}\n'.format(drnra.shape[0], 0)) + for kper in range(nper): + f.write('{} {}\n'.format(drnra.shape[0], 0)) + f.write('external 71\n') + + # write drn list + sfacdrn = 1.5 + fname = os.path.join(ws, 'original.drn.dat') + with open(fname, 'w') as f: + for kper in range(nper): + f.write('sfac {}\n'.format(sfacdrn)) + for k, i, j, stage, cond in drnra: + f.write( + '{} {} {} {} {}\n'.format(k + 1, i + 1, j + 1, stage, cond)) + + # rewrite wel + fname = os.path.join(ws, 'original.wel') + with open(fname, 'w') as f: + f.write('{} {}\n'.format(drnra.shape[0], 0)) + for kper in range(nper): + f.write('{} {}\n'.format(drnra.shape[0], 0)) + f.write('external 72 (binary)\n') + + # create the wells, but use an all float dtype to write a binary file + # use one-based values + weldt = np.dtype([('k', ' 0: dt = pak_type.get_empty(1, aux_names=aux_names, structured=model.structured).dtype - pak_parms = mfparbc.load(f, nppak, dt, model.verbose) - # pak_parms = mfparbc.load(f, nppak, len(dt.names)) + pak_parms = mfparbc.load(f, nppak, dt, model, ext_unit_dict, + model.verbose) if nper is None: nrow, ncol, nlay, nper = model.get_nrow_ncol_nlay_nper() @@ -827,69 +840,8 @@ def load(f, model, pak_type, ext_unit_dict=None, **kwargs): elif itmp > 0: current = pak_type.get_empty(itmp, aux_names=aux_names, structured=model.structured) - for ibnd in range(itmp): - line = f.readline() - if "open/close" in line.lower(): - binary = False - if '(binary)' in line.lower(): - binary = True - # need to strip out existing path seps and - # replace current-system path seps - raw = line.strip().split() - fname = raw[1] - if '/' in fname: - raw = fname.split('/') - elif '\\' in fname: - raw = fname.split('\\') - else: - raw = [fname] - fname = os.path.join(*raw) - oc_filename = os.path.join(model.model_ws, fname) - msg = 'Package.load() error: open/close filename ' + \ - oc_filename + ' not found' - assert os.path.exists(oc_filename), msg - try: - if binary: - dtype2 = [] - for name in current.dtype.names: - dtype2.append((name, np.float32)) - dtype2 = np.dtype(dtype2) - d = np.fromfile(oc_filename, - dtype=dtype2, - count=itmp) - current = np.array(d, dtype=current.dtype) - else: - cd = current.dtype - current = np.loadtxt(oc_filename).transpose() - if current.ndim == 1: - current = np.atleast_2d( - current).transpose() - current = np.core.records.fromarrays(current, - dtype=cd) - current = current.view(np.recarray) - except Exception as e: - msg = 'Package.load() error loading ' + \ - 'open/close file ' + oc_filename + \ - ': ' + str(e) - raise Exception(msg) - msg = 'Package.load() error: open/close ' + \ - 'recarray from file ' + oc_filename + \ - ' shape (' + str(current.shape) + \ - ') does not match itmp: {:d}'.format(itmp) - assert current.shape[0] == itmp, msg - break - try: - t = line.strip().split() - current[ibnd] = tuple(t[:len(current.dtype.names)]) - except: - t = [] - for ivar in range(len(current.dtype.names)): - istart = ivar * 10 - istop = istart + 10 - t.append(line[istart:istop]) - current[ibnd] = tuple(t[:len(current.dtype.names)]) - - # convert indices to zero-based + current = ulstrd(f, itmp, current, model, sfac_columns, + ext_unit_dict) if model.structured: current['k'] -= 1 current['i'] -= 1 @@ -935,6 +887,7 @@ def load(f, model, pak_type, ext_unit_dict=None, **kwargs): # fill current parameter data (par_current) for ibnd, t in enumerate(data_dict): + t = tuple(t) par_current[ibnd] = tuple(t[:len(par_current.dtype.names)]) if model.structured: diff --git a/flopy/utils/flopy_io.py b/flopy/utils/flopy_io.py index 3c33c1844d..0901c7a889 100755 --- a/flopy/utils/flopy_io.py +++ b/flopy/utils/flopy_io.py @@ -1,6 +1,7 @@ """ Module for input/output utilities """ +import os import sys import numpy as np @@ -323,3 +324,132 @@ def get_url_text(url, error_msg=None): if error_msg is not None: print(error_msg) return + + +def ulstrd(f, nlist, ra, model, sfac_columns, ext_unit_dict): + """ + Read a list and allow for open/close, binary, external, sfac, etc. + + Parameters + ---------- + f : file handle + file handle for where the list is being read from + nlist : int + size of the list (number of rows) to read + ra : np.recarray + A record array of the correct size that will be filled with the list + model : model object + The model object (of type :class:`flopy.modflow.mf.Modflow`) to + which this package will be added. + sfac_columns : list + A list of strings containing the column names to scale by sfac + ext_unit_dict : dictionary, optional + If the list in the file is specified using EXTERNAL, + then in this case ext_unit_dict is required, which can be + constructed using the function + :class:`flopy.utils.mfreadnam.parsenamefile`. + + Returns + ------- + + """ + + # initialize variables + line = f.readline() + sfac = 1. + binary = False + ncol = len(ra.dtype.names) + line_list = line.strip().split() + close_the_file = False + file_handle = f + mode = 'r' + + # check for external + if line.strip().lower().startswith('external'): + inunit = int(line_list[1]) + errmsg = 'Could not find a file for unit {}'.format(inunit) + if ext_unit_dict is not None: + if inunit in ext_unit_dict: + namdata = ext_unit_dict[inunit] + file_handle = namdata.filehandle + else: + raise IOError(errmsg) + else: + raise IOError(errmsg) + if namdata.filetype == 'DATA(BINARY)': + binary = True + if not binary: + line = file_handle.readline() + + # or check for open/close + elif line.strip().lower().startswith('open/close'): + raw = line.strip().split() + fname = raw[1] + if '/' in fname: + raw = fname.split('/') + elif '\\' in fname: + raw = fname.split('\\') + else: + raw = [fname] + fname = os.path.join(*raw) + oc_filename = os.path.join(model.model_ws, fname) + msg = 'Package.load() error: open/close filename ' + \ + oc_filename + ' not found' + assert os.path.exists(oc_filename), msg + if '(binary)' in line.lower(): + binary = True + mode = 'rb' + file_handle = open(oc_filename, mode) + close_the_file = True + if not binary: + line = file_handle.readline() + + # check for scaling factor + if not binary: + line_list = line.strip().split() + if line.strip().lower().startswith('sfac'): + sfac = float(line_list[1]) + line = file_handle.readline() + + # fast binary read fromfile + if binary: + dtype2 = [] + for name in ra.dtype.names: + dtype2.append((name, np.float32)) + dtype2 = np.dtype(dtype2) + d = np.fromfile(file_handle, dtype=dtype2, count=nlist) + ra = np.array(d, dtype=ra.dtype) + ra = ra.view(np.recarray) + + # else, read ascii + else: + + for ii in range(nlist): + + # first line was already read + if ii != 0: + line = file_handle.readline() + + try: + # whitespace separated + t = line.strip().split() + ra[ii] = tuple(t[:ncol]) + except: + # fixed format + t = [] + for ivar in range(ncol): + istart = ivar * 10 + istop = istart + 10 + t.append(line[istart:istop]) + ra[ii] = tuple(t[:ncol]) + + # scale the data and check + for column_name in sfac_columns: + ra[column_name] *= sfac + if 'auxsfac' in ra.dtype.names: + ra[column_name] *= ra['auxsfac'] + + if close_the_file: + file_handle.close() + + return ra From 35253383743e7f192803b9c6b49aa178589db31f Mon Sep 17 00:00:00 2001 From: Christian Langevin Date: Tue, 22 Oct 2019 15:56:25 -0500 Subject: [PATCH 2/4] modified the flopy ulstrd to use model.free_format_input to determine how to read the line --- examples/data/parameters/Oahu_02.lpf | 3 +++ flopy/modflow/mfhfb.py | 4 +++- flopy/modflow/mfparbc.py | 3 +++ flopy/modflow/mfstr.py | 3 ++- flopy/utils/flopy_io.py | 18 ++++++++++-------- 5 files changed, 21 insertions(+), 10 deletions(-) diff --git a/examples/data/parameters/Oahu_02.lpf b/examples/data/parameters/Oahu_02.lpf index 20732b7cac..a7c6dcdc7a 100644 --- a/examples/data/parameters/Oahu_02.lpf +++ b/examples/data/parameters/Oahu_02.lpf @@ -9,3 +9,6 @@ LAY1_HK HK 1.0 1 1 HKL1 ALL LAY1_VKA VK 1.0 1 1 VKAL1 ALL +1 +1 +1 diff --git a/flopy/modflow/mfhfb.py b/flopy/modflow/mfhfb.py index c243df01ee..3c694fb865 100644 --- a/flopy/modflow/mfhfb.py +++ b/flopy/modflow/mfhfb.py @@ -298,7 +298,9 @@ def load(f, model, ext_unit_dict=None): # data set 2 and 3 if nphfb > 0: dt = ModflowHfb.get_empty(1).dtype - pak_parms = mfparbc.load(f, nphfb, dt, model.verbose) + pak_parms = mfparbc.load(f, nphfb, dt, model, + ext_unit_dict=ext_unit_dict, + verbose=model.verbose) # data set 4 bnd_output = None if nhfbnp > 0: diff --git a/flopy/modflow/mfparbc.py b/flopy/modflow/mfparbc.py index 7789a80082..e22b90bddf 100644 --- a/flopy/modflow/mfparbc.py +++ b/flopy/modflow/mfparbc.py @@ -103,6 +103,9 @@ def load(f, npar, dt, model, ext_unit_dict=None, verbose=False): instnam = 'static' ra = np.zeros(nlst, dtype=dt) + #todo: if sfac is used for parameter definition, then + # the empty list on the next line needs to be the package + # get_sfac_columns bcinst = ulstrd(f, nlst, ra, model, [], ext_unit_dict) pinst[instnam] = bcinst bc_parms[parnam] = [{'partyp': partyp, 'parval': parval, diff --git a/flopy/modflow/mfstr.py b/flopy/modflow/mfstr.py index 555352973b..c47010ce37 100644 --- a/flopy/modflow/mfstr.py +++ b/flopy/modflow/mfstr.py @@ -619,7 +619,8 @@ def load(f, model, nper=None, ext_unit_dict=None): # read parameter data if npstr > 0: dt = ModflowStr.get_empty(1, aux_names=aux_names).dtype - pak_parms = mfparbc.load(f, npstr, dt, model.verbose) + pak_parms = mfparbc.load(f, npstr, dt, model, ext_unit_dict, + model.verbose) if nper is None: nrow, ncol, nlay, nper = model.get_nrow_ncol_nlay_nper() diff --git a/flopy/utils/flopy_io.py b/flopy/utils/flopy_io.py index 0901c7a889..3b19b85735 100755 --- a/flopy/utils/flopy_io.py +++ b/flopy/utils/flopy_io.py @@ -430,17 +430,19 @@ def ulstrd(f, nlist, ra, model, sfac_columns, ext_unit_dict): if ii != 0: line = file_handle.readline() - try: + if model.free_format_input: # whitespace separated t = line.strip().split() - ra[ii] = tuple(t[:ncol]) - except: + if len(t) < ncol: + t = t + (ncol - len(t)) * [0.0] + else: + t = t[:ncol] + t = tuple(t) + ra[ii] = t + else: # fixed format - t = [] - for ivar in range(ncol): - istart = ivar * 10 - istop = istart + 10 - t.append(line[istart:istop]) + t = read_fixed_var(line, ncol=ncol) + t = tuple(t) ra[ii] = tuple(t[:ncol]) # scale the data and check From f1a33dddc667d288fbf5baae17bf3bc9bf042288 Mon Sep 17 00:00:00 2001 From: Christian Langevin Date: Tue, 22 Oct 2019 16:15:35 -0500 Subject: [PATCH 3/4] Minor clean up --- autotest/t003_test.py | 14 +++++++------- flopy/modflow/mf.py | 1 - flopy/utils/flopy_io.py | 2 +- 3 files changed, 8 insertions(+), 9 deletions(-) diff --git a/autotest/t003_test.py b/autotest/t003_test.py index 90c1c4031b..e547ad2147 100644 --- a/autotest/t003_test.py +++ b/autotest/t003_test.py @@ -108,11 +108,11 @@ def test_load_nam_mt_nonexistant_file(): if __name__ == '__main__': - #test_loadoc() - #test_loadoc_nstpfail() - #test_load_nam_mf_nonexistant_file() - #test_load_nam_mt_nonexistant_file() - #test_loadoc_lenfail() - #test_loadfreyberg() + test_loadoc() + test_loadoc_nstpfail() + test_load_nam_mf_nonexistant_file() + test_load_nam_mt_nonexistant_file() + test_loadoc_lenfail() + test_loadfreyberg() test_loadoahu() - #test_loadtwrip() + test_loadtwrip() diff --git a/flopy/modflow/mf.py b/flopy/modflow/mf.py index d540a6ffba..7450b7bb3a 100644 --- a/flopy/modflow/mf.py +++ b/flopy/modflow/mf.py @@ -433,7 +433,6 @@ def write_name_file(self): self.external_binflag, self.external_output): if u == 0: continue - # fr = os.path.relpath(f, self.model_ws) replace_text = '' if o: replace_text = 'REPLACE' diff --git a/flopy/utils/flopy_io.py b/flopy/utils/flopy_io.py index 3b19b85735..531abc7ddf 100755 --- a/flopy/utils/flopy_io.py +++ b/flopy/utils/flopy_io.py @@ -443,7 +443,7 @@ def ulstrd(f, nlist, ra, model, sfac_columns, ext_unit_dict): # fixed format t = read_fixed_var(line, ncol=ncol) t = tuple(t) - ra[ii] = tuple(t[:ncol]) + ra[ii] = t # scale the data and check for column_name in sfac_columns: From a685c3c48b98c683f970cfdfa55307bd0ccb3b74 Mon Sep 17 00:00:00 2001 From: Christian Langevin Date: Tue, 22 Oct 2019 16:19:36 -0500 Subject: [PATCH 4/4] t067 needed to have _test in the file name --- autotest/{t067_ulstrd.py => t067_test_ulstrd.py} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename autotest/{t067_ulstrd.py => t067_test_ulstrd.py} (100%) diff --git a/autotest/t067_ulstrd.py b/autotest/t067_test_ulstrd.py similarity index 100% rename from autotest/t067_ulstrd.py rename to autotest/t067_test_ulstrd.py