This notebook demonstrates how to use the ZoneBudget class to extract budget information from the cell by cell budget file using an array of zones.
First set the path and import the required packages. The flopy path doesn't have to be set if you install flopy from a binary installer. If you want to run this notebook, you have to set the path to your own flopy path.
In [1]:
%matplotlib inline
import os
import sys
import platform
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import pandas as pd
import flopy
print(sys.version)
print('numpy version: {}'.format(np.__version__))
print('matplotlib version: {}'.format(mpl.__version__))
print('pandas version: {}'.format(pd.__version__))
print('flopy version: {}'.format(flopy.__version__))
In [2]:
# Set path to example datafiles
loadpth = os.path.join('..', 'data', 'zonbud_examples')
cbc_f = os.path.join(loadpth, 'freyberg_mlt', 'freyberg.gitcbc')
In [3]:
from flopy.utils import read_zbarray
zone_file = os.path.join(loadpth, 'zonef_mlt')
zon = read_zbarray(zone_file)
nlay, nrow, ncol = zon.shape
fig = plt.figure(figsize=(10, 4))
for lay in range(nlay):
ax = fig.add_subplot(1, nlay, lay+1)
im = ax.pcolormesh(zon[lay, :, :])
cbar = plt.colorbar(im)
plt.gca().set_aspect('equal')
plt.show()
np.unique(zon)
Out[3]:
In [4]:
# Create a ZoneBudget object and get the budget record array
zb = flopy.utils.ZoneBudget(cbc_f, zon, kstpkper=(0, 1096))
zb.get_budget()
Out[4]:
In [5]:
# Get a list of the unique budget record names
zb.get_record_names()
Out[5]:
In [6]:
# Look at a subset of fluxes
names = ['RECHARGE_IN', 'ZONE_1_IN', 'ZONE_3_IN']
zb.get_budget(names=names)
Out[6]:
In [7]:
# Look at fluxes in from zone 2
names = ['RECHARGE_IN', 'ZONE_1_IN', 'ZONE_3_IN']
zones = ['ZONE_2']
zb.get_budget(names=names, zones=zones)
Out[7]:
In [8]:
# Look at all of the mass-balance records
names = ['TOTAL_IN', 'TOTAL_OUT', 'IN-OUT', 'PERCENT_DISCREPANCY']
zb.get_budget(names=names)
Out[8]:
In [9]:
cmd = flopy.utils.ZoneBudget(cbc_f, zon, kstpkper=(0, 0))
cfd = cmd / 35.3147
inyr = (cfd / (250 * 250)) * 365 * 12
cmdbud = cmd.get_budget()
cfdbud = cfd.get_budget()
inyrbud = inyr.get_budget()
names = ['RECHARGE_IN']
rowidx = np.in1d(cmdbud['name'], names)
colidx = 'ZONE_1'
print('{:,.1f} cubic meters/day'.format(cmdbud[rowidx][colidx][0]))
print('{:,.1f} cubic feet/day'.format(cfdbud[rowidx][colidx][0]))
print('{:,.1f} inches/year'.format(inyrbud[rowidx][colidx][0]))
In [10]:
cmd is cfd
Out[10]:
In [11]:
aliases = {1: 'SURF', 2:'CONF', 3: 'UFA'}
zb = flopy.utils.ZoneBudget(cbc_f, zon, totim=[1097.], aliases=aliases)
zb.get_budget()
Out[11]:
In [12]:
zon = np.ones((nlay, nrow, ncol), np.int)
zon[1, :, :] = 2
zon[2, :, :] = 3
aliases = {1: 'SURF', 2:'CONF', 3: 'UFA'}
zb = flopy.utils.ZoneBudget(cbc_f, zon, kstpkper=None, totim=None, aliases=aliases)
df = zb.get_dataframes()
print(df.head())
print(df.tail())
Slice the multi-index dataframe to retrieve a subset of the budget.
NOTE: We can now pass "names" directly to the get_dataframes() method to return a subset of reocrds. By omitting the "_IN" or "_OUT" suffix we get both.
In [13]:
dateidx1 = 1092.
dateidx2 = 1097.
names = ['RECHARGE_IN', 'WELLS_OUT', 'CONSTANT_HEAD']
zones = ['SURF', 'CONF']
df = zb.get_dataframes(names=names)
df.loc[(slice(dateidx1, dateidx2), slice(None)), :][zones]
Out[13]:
Look at pumpage (WELLS_OUT) as a percentage of recharge (RECHARGE_IN)
In [14]:
dateidx1 = 1092.
dateidx2 = 1097.
zones = ['SURF']
# Pull out the individual records of interest
rech = df.loc[(slice(dateidx1, dateidx2), ['RECHARGE_IN']), :][zones]
pump = df.loc[(slice(dateidx1, dateidx2), ['WELLS_OUT']), :][zones]
# Remove the "record" field from the index so we can
# take the difference of the two DataFrames
rech = rech.reset_index()
rech = rech.set_index(['totim'])
rech = rech[zones]
pump = pump.reset_index()
pump = pump.set_index(['totim'])
pump = pump[zones] * -1
# Compute pumping as a percentage of recharge
pump_as_pct = (pump / rech) * 100.
pump_as_pct
Out[14]:
Pass start_datetime and timeunit keyword arguments to return a dataframe with a datetime multi-index
In [15]:
dateidx1 = pd.Timestamp('1972-12-01')
dateidx2 = pd.Timestamp('1972-12-06')
names = ['RECHARGE_IN', 'WELLS_OUT', 'CONSTANT_HEAD']
zones = ['SURF', 'CONF']
df = zb.get_dataframes(start_datetime='1970-01-01', timeunit='D', names=names)
df.loc[(slice(dateidx1, dateidx2), slice(None)), :][zones]
Out[15]:
Pass index_key to indicate which fields to use in the multi-index (defualt is "totim"; valid keys are "totim" and "kstpkper")
In [16]:
df = zb.get_dataframes(index_key='kstpkper')
df.head()
Out[16]:
In [17]:
zb = flopy.utils.ZoneBudget(cbc_f, zon, kstpkper=[(0, 0), (0, 1096)])
zb.to_csv(os.path.join(loadpth, 'zonbud.csv'))
# Read the file in to see the contents
fname = os.path.join(loadpth, 'zonbud.csv')
try:
import pandas as pd
print(pd.read_csv(fname).to_string(index=False))
except:
with open(fname, 'r') as f:
for line in f.readlines():
print('\t'.join(line.split(',')))
In [18]:
zon = np.ones((nlay, nrow, ncol), np.int)
zon[1, :, :] = 2
zon[2, :, :] = 3
aliases = {1: 'SURF', 2:'CONF', 3: 'UFA'}
zb = flopy.utils.ZoneBudget(cbc_f, zon, kstpkper=None, totim=None, aliases=aliases)
cfd = zb.get_budget(names=['STORAGE', 'WELLS'], zones=['SURF', 'UFA'], net=True)
cfd
Out[18]:
In [19]:
df = zb.get_dataframes(names=['STORAGE', 'WELLS'], zones=['SURF', 'UFA'], net=True)
df.head(6)
Out[19]:
In [20]:
def tick_label_formatter_comma_sep(x, pos):
return '{:,.0f}'.format(x)
def volumetric_budget_bar_plot(values_in, values_out, labels, **kwargs):
if 'ax' in kwargs:
ax = kwargs.pop('ax')
else:
ax = plt.gca()
x_pos = np.arange(len(values_in))
rects_in = ax.bar(x_pos, values_in, align='center', alpha=0.5)
x_pos = np.arange(len(values_out))
rects_out = ax.bar(x_pos, values_out, align='center', alpha=0.5)
plt.xticks(list(x_pos), labels)
ax.set_xticklabels(ax.xaxis.get_majorticklabels(), rotation=90)
ax.get_yaxis().set_major_formatter(mpl.ticker.FuncFormatter(tick_label_formatter_comma_sep))
ymin, ymax = ax.get_ylim()
if ymax != 0:
if abs(ymin) / ymax < .33:
ymin = -(ymax * .5)
else:
ymin *= 1.35
else:
ymin *= 1.35
plt.ylim([ymin, ymax * 1.25])
for i, rect in enumerate(rects_in):
label = '{:,.0f}'.format(values_in[i])
height = values_in[i]
x = rect.get_x() + rect.get_width() / 2
y = height + (.02 * ymax)
vertical_alignment = 'bottom'
horizontal_alignment = 'center'
ax.text(x, y, label, ha=horizontal_alignment, va=vertical_alignment, rotation=90)
for i, rect in enumerate(rects_out):
label = '{:,.0f}'.format(values_out[i])
height = values_out[i]
x = rect.get_x() + rect.get_width() / 2
y = height + (.02 * ymin)
vertical_alignment = 'top'
horizontal_alignment = 'center'
ax.text(x, y, label, ha=horizontal_alignment, va=vertical_alignment, rotation=90)
# horizontal line indicating zero
ax.plot([rects_in[0].get_x() - rects_in[0].get_width() / 2,
rects_in[-1].get_x() + rects_in[-1].get_width()], [0, 0], "k")
return rects_in, rects_out
In [21]:
fig = plt.figure(figsize=(16, 5))
times = [2., 500., 1000., 1095.]
for idx, time in enumerate(times):
ax = fig.add_subplot(1, len(times), idx + 1)
zb = flopy.utils.ZoneBudget(cbc_f, zon, kstpkper=None, totim=time, aliases=aliases)
recname = 'STORAGE'
values_in = zb.get_dataframes(names='{}_IN'.format(recname)).T.squeeze()
values_out = zb.get_dataframes(names='{}_OUT'.format(recname)).T.squeeze() * -1
labels = values_in.index.tolist()
rects_in, rects_out = volumetric_budget_bar_plot(values_in, values_out, labels, ax=ax)
plt.ylabel('Volumetric rate, in Mgal/d')
plt.title('totim = {}'.format(time))
plt.tight_layout()
plt.show()