diff --git a/autotest/t041_test.py b/autotest/t041_test.py index 14a7366f2d..8a7b594739 100644 --- a/autotest/t041_test.py +++ b/autotest/t041_test.py @@ -217,7 +217,8 @@ def test_filenames(): pkglst = ['dis', 'bas6', 'pcg', 'lpf'] m = flopy.modflow.Modflow.load(modelname + '.nam', model_ws=pth, check=False, load_only=pkglst, - verbose=False, exe_name=exe_name) + verbose=False, exe_name=exe_name, + forgive=False) obs = flopy.modflow.HeadObservation(m, layer=0, row=5, column=5, time_series_data=[[1., 54.4], @@ -229,48 +230,6 @@ def test_filenames(): obs_data=[obs], options=['NOPRINT'], filenames=filenames) - - # add DRN package - spd = {0: [[0, 5, 5, .5, 8e6], - [0, 8, 8, .7, 8e6]]} - drn = flopy.modflow.ModflowDrn(m, 53, stress_period_data=spd) - - # flow observation - - # Lists of length nqfb - nqobfb = [1, 1] - nqclfb = [1, 1] - - # Lists of length nqtfb - obsnam = ['drob_1', 'drob_2'] - irefsp = [0, 0] - toffset = [0, 0] - flwobs = [0., 0.] - - # Lists of length (nqfb, nqclfb) - layer = [[0], [0]] - row = [[5], [8]] - column = [[5], [8]] - factor = [[1.], [1.]] - - drob = flopy.modflow.ModflowFlwob(m, - nqfb=len(nqclfb), - nqcfb=np.sum(nqclfb), - nqtfb=np.sum(nqobfb), - nqobfb=nqobfb, - nqclfb=nqclfb, - obsnam=obsnam, - irefsp=irefsp, - toffset=toffset, - flwobs=flwobs, - layer=layer, - row=row, - column=column, - factor=factor, - flowtype='drn', - options=['NOPRINT'], - filenames=['flwobs_simple.drob', - 'flwobs_simple.obd']) # Write the model input files m.write_input() @@ -315,8 +274,95 @@ def test_multilayerhob_prfail(): return +def test_flwob_load(): + """ + test041 create, write, and load ModflowFlwob package. + """ + # load the modflow model + opth = os.path.join(cpth, 'tc1-true', 'orig') + m = flopy.modflow.Modflow.load('tc1-true.nam', verbose=True, + model_ws=opth, exe_name=exe_name) + + npth = os.path.join(cpth, 'tc1-true', 'flwob') + m.change_model_ws(new_pth=npth, reset_external=True) + + # write the lgr model in to the new path + m.write_input() + + # add DRN package + spd = {0: [[0, 5, 5, .5, 8e6], + [0, 8, 8, .7, 8e6]]} + drn = flopy.modflow.ModflowDrn(m, 53, stress_period_data=spd) + + # flow observation + + # Lists of length nqfb + nqobfb = [1, 1] + nqclfb = [1, 1] + + # Lists of length nqtfb + obsnam = ['drob_1', 'drob_2'] + irefsp = [0, 0] + toffset = [0, 0] + flwobs = [-5.678, -6.874] + + # Lists of length (nqfb, nqclfb) + layer = [[0], [0]] + row = [[5], [8]] + column = [[5], [8]] + factor = [[1.], [1.]] + + drob = flopy.modflow.ModflowFlwob(m, + nqfb=len(nqclfb), + nqcfb=np.sum(nqclfb), + nqtfb=np.sum(nqobfb), + nqobfb=nqobfb, + nqclfb=nqclfb, + obsnam=obsnam, + irefsp=irefsp, + toffset=toffset, + flwobs=flwobs, + layer=layer, + row=row, + column=column, + factor=factor, + flowtype='drn', + options=['NOPRINT']) + # Write the model input files + m.write_input(check=False) + + # Load the DROB package and compare it to the original + pkglst = ['drob'] + m = flopy.modflow.Modflow.load('tc1-true.nam', model_ws=npth, + check=False, load_only=pkglst, + verbose=False, exe_name=exe_name, + forgive=False) + + # check variables were read properly + s = 'nqfb loaded from {} read incorrectly'.format(m.drob.fn_path) + assert(drob.nqfb == m.drob.nqfb), s + s = 'nqcfb loaded from {} read incorrectly'.format(m.drob.fn_path) + assert (drob.nqcfb == m.drob.nqcfb), s + s = 'nqtfb loaded from {} read incorrectly'.format(m.drob.fn_path) + assert (drob.nqtfb == m.drob.nqtfb), s + s = 'obsnam loaded from {} read incorrectly'.format(m.drob.fn_path) + assert (list([n for n in drob.obsnam]) == + list([n for n in m.drob.obsnam])), s + s = 'flwobs loaded from {} read incorrectly'.format(m.drob.fn_path) + assert np.array_equal(drob.flwobs, m.drob.flwobs), s + s = 'layer loaded from {} read incorrectly'.format(m.drob.fn_path) + assert np.array_equal(drob.layer, m.drob.layer), s + s = 'row loaded from {} read incorrectly'.format(m.drob.fn_path) + assert np.array_equal(drob.row, m.drob.row), s + s = 'column loaded from {} read incorrectly'.format(m.drob.fn_path) + assert np.array_equal(drob.column, m.drob.column), s + + return + + if __name__ == '__main__': test_hob_simple() test_obs_create_and_write() test_obs_load_and_write() test_filenames() + test_flwob_load() diff --git a/flopy/modflow/mf.py b/flopy/modflow/mf.py index 4b58e23fa6..a930539565 100644 --- a/flopy/modflow/mf.py +++ b/flopy/modflow/mf.py @@ -203,6 +203,10 @@ def __init__(self, modelname='modflowtest', namefile_ext='nam', "swt": flopy.modflow.ModflowSwt, "hyd": flopy.modflow.ModflowHyd, "hob": flopy.modflow.ModflowHob, + "chob": flopy.modflow.ModflowFlwob, + "gbob": flopy.modflow.ModflowFlwob, + "drob": flopy.modflow.ModflowFlwob, + "rvob": flopy.modflow.ModflowFlwob, "vdf": flopy.seawat.SeawatVdf, "vsc": flopy.seawat.SeawatVsc } diff --git a/flopy/modflow/mfflwob.py b/flopy/modflow/mfflwob.py index d47247d5ad..beeb08f28d 100755 --- a/flopy/modflow/mfflwob.py +++ b/flopy/modflow/mfflwob.py @@ -1,5 +1,8 @@ +import os +import sys import numpy as np from ..pakbase import Package +from ..utils import parsenamefile class ModflowFlwob(Package): @@ -135,41 +138,67 @@ def __init__(self, model, nqfb=0, nqcfb=0, nqtfb=0, iufbobsv=0, if extension is None: extension = ['chob', 'obc', 'gbob', 'obg', 'drob', 'obd', 'rvob', 'obr'] - if unitnumber is None: - unitnumber = [40, 140, 41, 141, 42, 142, 43, 143] + pakunits = {'chob': 40, + 'gbob': 41, + 'drob': 42, + 'rvob': 43} + outunits = {'chob': 140, + 'gbob': 141, + 'drob': 142, + 'rvob': 143} + # if unitnumber is None: + # unitnumber = [40, 140, 41, 141, 42, 142, 43, 143] if flowtype.upper().strip() == 'CHD': name = ['CHOB', 'DATA'] extension = extension[0:2] - unitnumber = unitnumber[0:2] - iufbobsv = unitnumber[1] + # unitnumber = unitnumber[0:2] + # iufbobsv = unitnumber[1] + self._ftype = 'CHOB' self.url = 'chob.htm' self.heading = '# CHOB for MODFLOW, generated by Flopy.' elif flowtype.upper().strip() == 'GHB': name = ['GBOB', 'DATA'] extension = extension[2:4] - unitnumber = unitnumber[2:4] - iufbobsv = unitnumber[1] + # unitnumber = unitnumber[2:4] + # iufbobsv = unitnumber[1] + self._ftype = 'GBOB' self.url = 'gbob.htm' self.heading = '# GBOB for MODFLOW, generated by Flopy.' elif flowtype.upper().strip() == 'DRN': name = ['DROB', 'DATA'] extension = extension[4:6] - unitnumber = unitnumber[4:6] - iufbobsv = unitnumber[1] + # unitnumber = unitnumber[4:6] + # iufbobsv = unitnumber[1] + self._ftype = 'DROB' self.url = 'drob.htm' self.heading = '# DROB for MODFLOW, generated by Flopy.' elif flowtype.upper().strip() == 'RIV': name = ['RVOB', 'DATA'] extension = extension[6:8] - unitnumber = unitnumber[6:8] - iufbobsv = unitnumber[1] + # unitnumber = unitnumber[6:8] + # iufbobsv = unitnumber[1] + self._ftype = 'RVOB' self.url = 'rvob.htm' self.heading = '# RVOB for MODFLOW, generated by Flopy.' else: msg = 'ModflowFlwob: flowtype must be CHD, GHB, DRN, or RIV' raise KeyError(msg) + if unitnumber is None: + unitnumber = [pakunits[name[0].lower()], + outunits[name[0].lower()]] + elif isinstance(unitnumber, int): + unitnumber = [unitnumber] + if len(unitnumber) == 1: + if unitnumber[0] in outunits.keys(): + unitnumber = [pakunits[name[0].lower()], + unitnumber[0]] + else: + unitnumber = [unitnumber[0], + outunits[name[0].lower()]] + iufbobsv = unitnumber[1] + # set filenames if filenames is None: filenames = [None, None] @@ -201,10 +230,14 @@ def __init__(self, model, nqfb=0, nqcfb=0, nqtfb=0, iufbobsv=0, self.factor = factor # -create empty arrays of the correct size - self.layer = np.zeros((self.nqfb, max(self.nqclfb)), dtype='int32') - self.row = np.zeros((self.nqfb, max(self.nqclfb)), dtype='int32') - self.column = np.zeros((self.nqfb, max(self.nqclfb)), dtype='int32') - self.factor = np.zeros((self.nqfb, max(self.nqclfb)), dtype='float32') + self.layer = np.zeros((self.nqfb, max(np.abs(self.nqclfb))), + dtype='int32') + self.row = np.zeros((self.nqfb, max(np.abs(self.nqclfb))), + dtype='int32') + self.column = np.zeros((self.nqfb, max(np.abs(self.nqclfb))), + dtype='int32') + self.factor = np.zeros((self.nqfb, max(np.abs(self.nqclfb))), + dtype='float32') self.nqobfb = np.zeros((self.nqfb), dtype='int32') self.nqclfb = np.zeros((self.nqfb), dtype='int32') self.irefsp = np.zeros((self.nqtfb), dtype='int32') @@ -238,6 +271,9 @@ def __init__(self, model, nqfb=0, nqcfb=0, nqtfb=0, iufbobsv=0, # add checks for input compliance (obsnam length, etc.) self.parent.add_package(self) + def ftype(self): + return self._ftype + def write_file(self): """ Write the package file @@ -275,10 +311,10 @@ def write_file(self): # Loop through observation times for the groups for j in range(self.nqobfb[i]): # write section 4 - line = '{}{:10d}{:10.4g} {:10.4g}\n'.format(self.obsnam[c], - self.irefsp[c] + 1, - self.toffset[c], - self.flwobs[c]) + line = '{:12}'.format(self.obsnam[c]) + line += '{:8d}'.format(self.irefsp[c] + 1) + line += '{:16.10g}'.format(self.toffset[c]) + line += ' {:10.4g}\n'.format(self.flwobs[c]) f_fbob.write(line) c += 1 # index variable @@ -314,3 +350,230 @@ def write_file(self): # swm: END hack for writing standard file return + + @staticmethod + def load(f, model, ext_unit_dict=None, check=True): + """ + Load an existing package. + + Parameters + ---------- + f : filename or file handle + File to load. + model : model object + The model object (of type :class:`flopy.modflow.mf.Modflow`) to + which this package will be added. + ext_unit_dict : dictionary, optional + If the arrays in the file are specified using EXTERNAL, + or older style array control records, then `f` should be a file + handle. In this case ext_unit_dict is required, which can be + constructed using the function + :class:`flopy.utils.mfreadnam.parsenamefile`. + check : boolean + Check package data for common errors. (default True) + + Returns + ------- + flwob : ModflowFlwob package object + ModflowFlwob package object. + + Examples + -------- + + >>> import flopy + >>> m = flopy.modflow.Modflow() + >>> hobs = flopy.modflow.ModflowFlwob.load('test.drob', m) + + """ + + if model.verbose: + sys.stdout.write('loading flwob package file...\n') + + if not hasattr(f, 'read'): + filename = f + f = open(filename, 'r') + + # dataset 0 -- header + while True: + line = f.readline() + if line[0] != '#': + break + + # read dataset 1 -- NQFB NQCFB NQTFB IUFBOBSV Options + t = line.strip().split() + nqfb = int(t[0]) + nqcfb = int(t[1]) + nqtfb = int(t[2]) + iufbobsv = int(t[3]) + options = [] + if len(t) > 4: + options = t[4:] + + # read dataset 2 -- TOMULTFB + line = f.readline() + t = line.strip().split() + tomultfb = float(t[0]) + + nqobfb = np.zeros(nqfb, dtype=np.int32) + nqclfb = np.zeros(nqfb, dtype=np.int32) + obsnam = [] + irefsp = [] + toffset = [] + flwobs = [] + + layer = [] + row = [] + column = [] + factor = [] + + # read datasets 3, 4, and 5 for each of nqfb groups + # of cells + nobs = 0 + while True: + + # read dataset 3 -- NQOBFB NQCLFB + line = f.readline() + t = line.strip().split() + nqobfb[nobs] = int(t[0]) + nqclfb[nobs] = int(t[1]) + + # read dataset 4 -- OBSNAM IREFSP TOFFSET FLWOBS + ntimes = 0 + while True: + line = f.readline() + t = line.strip().split() + obsnam.append(t[0]) + irefsp.append(int(t[1])) + toffset.append(float(t[2])) + flwobs.append(float(t[3])) + ntimes += 1 + if ntimes == nqobfb[nobs]: + break + + # read dataset 5 -- Layer Row Column Factor + k = np.zeros(abs(nqclfb[nobs]), np.int32) + i = np.zeros(abs(nqclfb[nobs]), np.int32) + j = np.zeros(abs(nqclfb[nobs]), np.int32) + fac = np.zeros(abs(nqclfb[nobs]), np.float32) + + ncells = 0 + while True: + line = f.readline() + t = line.strip().split() + k[ncells] = int(t[0]) + i[ncells] = int(t[1]) + j[ncells] = int(t[2]) + fac[ncells] = float(t[3]) + + ncells += 1 + if ncells == abs(nqclfb[nobs]): + layer.append(k) + row.append(i) + column.append(j) + factor.append(fac) + break + + nobs += 1 + if nobs == nqfb: + break + + irefsp = np.array(irefsp) - 1 + layer = np.array(layer) - 1 + row = np.array(row) - 1 + column = np.array(column) - 1 + factor = np.array(factor) + + # close the file + f.close() + + # get ext_unit_dict if none passed + if ext_unit_dict is None: + namefile = os.path.join(model.model_ws, model.namefile) + ext_unit_dict = parsenamefile(namefile, model.mfnam_packages) + + flowtype, ftype = _get_ftype_from_filename(f.name, ext_unit_dict) + + # set package unit number + unitnumber = None + filenames = [None, None] + if ext_unit_dict is not None: + unitnumber, filenames[0] = \ + model.get_ext_dict_attr(ext_unit_dict, + filetype=ftype.upper()) + if iufbobsv > 0: + iu, filenames[1] = \ + model.get_ext_dict_attr(ext_unit_dict, unit=iufbobsv) + model.add_pop_key_list(iufbobsv) + + # create ModflowFlwob object instance + flwob = ModflowFlwob(model, iufbobsv=iufbobsv, tomultfb=tomultfb, + nqfb=nqfb, nqcfb=nqcfb, + nqtfb=nqtfb, nqobfb=nqobfb, nqclfb=nqclfb, + obsnam=obsnam, irefsp=irefsp, toffset=toffset, + flwobs=flwobs, layer=layer, row=row, + column=column, factor=factor, options=options, + flowtype=flowtype, unitnumber=unitnumber, + filenames=filenames) + + return flwob + + +def _get_ftype_from_filename(fn, ext_unit_dict=None): + """ + Returns the boundary flowtype and filetype for a given ModflowFlwob + package filename. + + Parameters + ---------- + fn : str + The filename to be parsed. + ext_unit_dict : dictionary, optional + If the arrays in the file are specified using EXTERNAL, + or older style array control records, then `f` should be a file + handle. In this case ext_unit_dict is required, which can be + constructed using the function + :class:`flopy.utils.mfreadnam.parsenamefile`. + + Returns + ------- + flowtype : str + Corresponds to the type of the head-dependent boundary package for + which observations are desired (e.g. "CHD", "GHB", "DRN", or "RIV"). + ftype : str + Corresponds to the observation file type (e.g. "CHOB", "GBOB", + "DROB", or "RVOB"). + """ + + ftype = None + + # determine filetype from filename using ext_unit_dict + if ext_unit_dict is not None: + for key, value in ext_unit_dict.items(): + if value.filename == fn: + ftype = value.filetype + break + + # else, try to infer filetype from filename extension + else: + ext = fn.split('.')[-1].lower() + if 'ch' in ext.lower(): + ftype = 'CHOB' + elif 'gb' in ext.lower(): + ftype = 'GBOB' + elif 'dr' in ext.lower(): + ftype = 'DROB' + elif 'rv' in ext.lower(): + ftype = 'RVOB' + + msg = 'ModflowFlwob: filetype cannot be inferred ' \ + 'from file name {}'.format(fn) + if ftype is None: + raise AssertionError(msg) + + flowtype_dict = {'CHOB': 'CHD', + 'GOBO': 'GHB', + 'DROB': 'DRN', + 'RVOB': 'RIV'} + flowtype = flowtype_dict[ftype] + + return flowtype, ftype