Skip to content

Commit

Permalink
fix(binaryfile_utils, CellbudgetFile): update for IMETH == 6 (#823)
Browse files Browse the repository at this point in the history 
* Update(plot_bc): updated plot_bc for MAW, UZF, SFR, and multiple bc packages

* fixed bug with mf2k loading #802
* added shapefile support for output_helper #807
* updated test cases

* Fix(ModflowSfr2.export): Update ModflowSfr2.export() to use modelgrid instead of SpatialReference

* Fix(MFOutputRequester): full3D=False when imeth == 6

* Update(CellBudgetFile): added methods to determine CellBudgetFile precision automatically

* Fix(CellBudgetFile.get_ts()): update get_ts() for IMETH==6.

* Codacy update
  • Loading branch information
@jlarsen-usgs
jlarsen-usgs committed Apr 10, 2020
1 parent 0336b5c commit f514a22
Show file tree
Hide file tree
Showing 3 changed files with 182 additions and 67 deletions.
8 changes: 4 additions & 4 deletions autotest/t504_test.py
View file
Expand Up @@ -69,7 +69,7 @@ def test001a_tharmonic():

# get expected results
budget_file = os.path.join(os.getcwd(), expected_cbc_file_a)
budget_obj = bf.CellBudgetFile(budget_file, precision='double')
budget_obj = bf.CellBudgetFile(budget_file, precision='auto')
budget_obj.list_records()
budget_frf_valid = np.array(budget_obj.get_data(text=' FLOW JA FACE', full3D=True))

Expand Down Expand Up @@ -112,7 +112,7 @@ def test001a_tharmonic():

# get expected results
budget_file = os.path.join(os.getcwd(), expected_cbc_file_b)
budget_obj = bf.CellBudgetFile(budget_file, precision='double')
budget_obj = bf.CellBudgetFile(budget_file, precision='auto')
budget_frf_valid = np.array(budget_obj.get_data(text=' FLOW JA FACE', full3D=True))

# compare output to expected results
Expand Down Expand Up @@ -167,7 +167,7 @@ def test003_gwfs_disv():

# get expected results
budget_file = os.path.join(os.getcwd(), expected_cbc_file_a)
budget_obj = bf.CellBudgetFile(budget_file, precision='double')
budget_obj = bf.CellBudgetFile(budget_file, precision='auto')
budget_fjf_valid = np.array(budget_obj.get_data(text=' FLOW JA FACE', full3D=True))

head_file = os.path.join(os.getcwd(), expected_head_file_a)
Expand Down Expand Up @@ -357,7 +357,7 @@ def test006_gwf3():

# get expected results
budget_file = os.path.join(os.getcwd(), expected_cbc_file_b)
budget_obj = bf.CellBudgetFile(budget_file, precision='double')
budget_obj = bf.CellBudgetFile(budget_file, precision='auto')
budget_fjf_valid = np.array(budget_obj.get_data(text=' FLOW JA FACE', full3D=True))
jaentries = budget_fjf_valid.shape[-1]
budget_fjf_valid.shape = (-1, jaentries)
Expand Down
62 changes: 38 additions & 24 deletions flopy/mf6/utils/binaryfile_utils.py
View file
Expand Up @@ -48,25 +48,28 @@ def __getitem__(self, index):

class MFOutputRequester:
"""
MFOutputRequest class is a helper function to enable the user to query
binary data from the SimulationDict() object on the fly without
actually storing it in the SimulationDict() object.
MFOutputRequest class is a helper function to enable the user to query
binary data from the SimulationDict() object on the fly without
actually storing it in the SimulationDict() object.
Parameters:
----------
mfdict: local instance of the SimulationDict() object
path: pointer to the MFSimulationPath object
key: user requested data key
Methods:
-------
MFOutputRequester.querybinarydata
returns: Xarray object
Parameters:
----------
mfdict: OrderedDict
local instance of the SimulationDict() object
path:
pointer to the MFSimulationPath object
key: tuple
user requested data key
Methods:
-------
MFOutputRequester.querybinarydata
returns: Xarray object
Examples:
--------
>>> data = MFOutputRequester(mfdict, path, key)
>>> data.querybinarydata
Examples:
--------
>>> data = MFOutputRequester(mfdict, path, key)
>>> data.querybinarydata
"""

def __init__(self, mfdict, path, key):
Expand Down Expand Up @@ -103,8 +106,11 @@ def _querybinarydata(self, key):
bindata = self._get_binary_file_object(path, bintype, key)

if bintype == 'CBC':
return np.array(bindata.get_data(text=key[-1], full3D=True))

try:
return np.array(bindata.get_data(text=key[-1], full3D=True))
except ValueError:
# imeth == 6
return np.array(bindata.get_data(text=key[-1], full3D=False))
else:
return np.array(bindata.get_alldata())

Expand All @@ -119,17 +125,22 @@ def _querybinarydata_vertices(self, mfdict, key):
if bintype == 'CBC':
if key[-1] == 'FLOW-JA-FACE':
data = np.array(bindata.get_data(text=key[-1]))
# todo: uncomment line to remove unnecessary dimensions from
# uncomment line to remove extra dimensions from data
# data data.shape = (len(times), -1)
return data

else:
data = np.array(bindata.get_data(text=key[-1], full3D=True))

try:
data = np.array(bindata.get_data(text=key[-1],
full3D=True))
except ValueError:
# imeth == 6
data = np.array(bindata.get_data(text=key[-1],
full3D=False))
else:
data = np.array(bindata.get_alldata())

# todo: uncomment line to remove extra dimensions from data
# uncomment line to remove extra dimensions from data
# data = _reshape_binary_data(data, 'V')
return data

Expand All @@ -141,7 +152,10 @@ def _querybinarydata_unstructured(self, key):
bindata = self._get_binary_file_object(path, bintype, key)

if bintype == 'CBC':
data = np.array(bindata.get_data(text=key[-1], full3D=True))
try:
data = np.array(bindata.get_data(text=key[-1], full3D=True))
except ValueError:
data = np.array(bindata.get_data(text=key[-1], full3D=False))
else:
data = bindata.get_alldata()

Expand Down
179 changes: 140 additions & 39 deletions flopy/utils/binaryfile.py
View file
Expand Up @@ -529,6 +529,10 @@ def __init__(self, filename, text='concentration', precision='auto',
return


class BudgetIndexError(Exception):
pass


class CellBudgetFile(object):
"""
CellBudgetFile Class.
Expand Down Expand Up @@ -564,7 +568,7 @@ class CellBudgetFile(object):
"""

def __init__(self, filename, precision='single', verbose=False, **kwargs):
def __init__(self, filename, precision='auto', verbose=False, **kwargs):
self.filename = filename
self.precision = precision
self.verbose = verbose
Expand All @@ -589,50 +593,55 @@ def __init__(self, filename, precision='single', verbose=False, **kwargs):
self.imethlist = []
self.paknamlist = []
self.nrecords = 0
h1dt = [('kstp', 'i4'), ('kper', 'i4'), ('text', 'a16'),
('ncol', 'i4'), ('nrow', 'i4'), ('nlay', 'i4')]

if precision == 'single':
self.realtype = np.float32
ffmt = 'f4'
elif precision == 'double':
self.realtype = np.float64
ffmt = 'f8'
else:
raise Exception('Unknown precision specified: ' + precision)
h2dt0 = [('imeth', 'i4'), ('delt', ffmt), ('pertim', ffmt),
('totim', ffmt)]
h2dt = [('imeth', 'i4'), ('delt', ffmt), ('pertim', ffmt),
('totim', ffmt), ('modelnam', 'a16'), ('paknam', 'a16'),
('modelnam2', 'a16'), ('paknam2', 'a16')]

self.dis = None
self.sr = None
self.modelgrid = None
if 'model' in kwargs.keys():
self.model = kwargs.pop('model')
self.sr = self.model.sr
self.modelgrid = self.model.modelgrid
self.dis = self.model.dis
if 'dis' in kwargs.keys():
self.dis = kwargs.pop('dis')
self.sr = self.dis.parent.sr
self.modelgrid = self.dis.parent.modelgrid
if 'sr' in kwargs.keys():
self.sr = kwargs.pop('sr')
from ..utils import SpatialReferenceUnstructured
from ..discretization import StructuredGrid, UnstructuredGrid
sr = kwargs.pop('sr')
if isinstance(sr, SpatialReferenceUnstructured):
self.modelgrid = UnstructuredGrid(vertices=sr.verts,
iverts=sr.iverts,
xcenters=sr.xc,
ycenters=sr.yc,
ncpl=sr.ncpl)
else:
self.modelgrid = StructuredGrid(delc=sr.delc, delr=sr.delr,
xoff=sr.xll, yoff=sr.yll,
angrot=sr.rotation)
if 'modelgrid' in kwargs.keys():
self.modelgrid = kwargs.pop('modelgrid')
if len(kwargs.keys()) > 0:
args = ','.join(kwargs.keys())
raise Exception('LayerFile error: unrecognized kwargs: ' + args)

self.header1_dtype = np.dtype(h1dt)
self.header2_dtype0 = np.dtype(h2dt0)
self.header2_dtype = np.dtype(h2dt)
hdt = h1dt + h2dt
self.header_dtype = np.dtype(hdt)
if precision == 'auto':
success = self._set_precision('single')
if not success:
success = self._set_precision('double')
if not success:
s = "Budget precision could not be auto determined"
raise BudgetIndexError(s)
elif precision == 'single':
success = self._set_precision(precision)
elif precision == 'double':
success = self._set_precision(precision)
else:
raise Exception('Unknown precision specified: ' + precision)

# read through the file and build the pointer index
self._build_index()
if not success:
s = "Budget file could not be read using {} " \
"precision".format(precision)
raise Exception(s)

# allocate the value array
# self.value = np.empty((self.nlay, self.nrow, self.ncol),
# dtype=self.realtype)
return

def __enter__(self):
Expand All @@ -641,6 +650,60 @@ def __enter__(self):
def __exit__(self, *exc):
self.close()

def __reset(self):
"""
Reset indexing lists when determining precision
"""
self.file.seek(0, 0)
self.times = []
self.kstpkper = []
self.recordarray = []
self.iposheader = []
self.iposarray = []
self.textlist = []
self.imethlist = []
self.paknamlist = []
self.nrecords = 0

def _set_precision(self, precision='single'):
"""
Method to set the budget precsion from a CBC file. Enables
Auto precision code to work
Parameters
----------
precision : str
budget file precision (accepts 'single' or 'double')
"""
success = True
h1dt = [('kstp', 'i4'), ('kper', 'i4'), ('text', 'a16'),
('ncol', 'i4'), ('nrow', 'i4'), ('nlay', 'i4')]
if precision == 'single':
self.realtype = np.float32
ffmt = 'f4'
else:
self.realtype = np.float64
ffmt = 'f8'

h2dt0 = [('imeth', 'i4'), ('delt', ffmt), ('pertim', ffmt),
('totim', ffmt)]
h2dt = [('imeth', 'i4'), ('delt', ffmt), ('pertim', ffmt),
('totim', ffmt), ('modelnam', 'a16'), ('paknam', 'a16'),
('modelnam2', 'a16'), ('paknam2', 'a16')]
self.header1_dtype = np.dtype(h1dt)
self.header2_dtype0 = np.dtype(h2dt0)
self.header2_dtype = np.dtype(h2dt)
hdt = h1dt + h2dt
self.header_dtype = np.dtype(hdt)

try:
self._build_index()
except BudgetIndexError:
success = False
self.__reset()

return success

def _totim_from_kstpkper(self, kstpkper):
if self.dis is None:
return -1.0
Expand All @@ -667,6 +730,10 @@ def _build_index(self):
Build the ordered dictionary, which maps the header information
to the position in the binary file.
"""
asciiset = ' '
for i in range(33, 127):
asciiset += chr(i)

header = self._get_header()
self.nrow = header["nrow"]
self.ncol = header["ncol"]
Expand Down Expand Up @@ -696,6 +763,16 @@ def _build_index(self):
if kstpkper not in self.kstpkper:
self.kstpkper.append(kstpkper)
if header['text'] not in self.textlist:
# check the precision of the file using text records
try:
text = header['text']
if isinstance(text, bytes):
text = text.decode()
for t in text:
if t.upper() not in asciiset:
raise Exception()
except:
raise BudgetIndexError("Improper precision")
self.textlist.append(header['text'])
self.imethlist.append(header['imeth'])
if header['paknam'] not in self.paknamlist:
Expand Down Expand Up @@ -1214,14 +1291,38 @@ def get_ts(self, idx, text=None, times=None):
result[idx, 0] = t

for itim, k in enumerate(kk):
v = self.get_data(kstpkper=k, text=text, full3D=True)
# skip missing data - required for storage
if len(v) > 0:
v = v[0]
istat = 1
for k, i, j in kijlist:
result[itim, istat] = v[k, i, j].copy()
istat += 1
try:
v = self.get_data(kstpkper=k, text=text, full3D=True)
# skip missing data - required for storage
if len(v) > 0:
v = v[0]
istat = 1
for k, i, j in kijlist:
result[itim, istat] = v[k, i, j].copy()
istat += 1
except ValueError:
v = self.get_data(kstpkper=k, text=text)
# skip missing data - required for storage
if len(v)> 0:
if self.modelgrid is None:
s = "A modelgrid instance must be provided during " \
"instantiation to get IMETH=6 timeseries data"
raise AssertionError(s)

if self.modelgrid.grid_type == 'structured':
ndx = [lrc[0] * (self.modelgrid.nrow *
self.modelgrid.ncol) +
lrc[1] * self.modelgrid.ncol +
(lrc[2] + 1) for lrc in kijlist]
else:
ndx = [lrc[0] * self.modelgrid.ncpl +
(lrc[-1] + 1) for lrc in kijlist]

for vv in v:
field = vv.dtype.names[2]
dix = np.where(np.isin(vv['node'], ndx))[0]
if len(dix) > 0:
result[itim, 1:] = vv[field][dix]

return result

Expand Down

0 comments on commit f514a22

Please sign in to comment.