Demonstrates functionality of Flopy SFR module using the example documented by Prudic and others (2004):
In [1]:
import sys
import platform
import os
import numpy as np
import glob
import shutil
import matplotlib as mpl
import matplotlib.pyplot as plt
import flopy
import flopy.utils.binaryfile as bf
In [2]:
#Set name of MODFLOW exe
# assumes executable is in users path statement
exe_name = 'mf2005'
if platform.system() == 'Windows':
exe_name = 'mf2005.exe'
% matplotlib inline
mpl.rcParams['figure.figsize'] = (11, 8.5)
In [3]:
path = 'data'
gpth = os.path.join('..', 'data', 'mf2005_test', 'test1ss.*')
for f in glob.glob(gpth):
shutil.copy(f, path)
In [4]:
m = flopy.modflow.Modflow.load('test1ss.nam', version='mf2005', exe_name=exe_name,
model_ws=path, load_only=['ghb', 'evt', 'rch', 'dis', 'bas6', 'oc', 'sip', 'lpf'])
Reach data (Item 2 in the SFR input instructions), are input and stored in a numpy record array
http://docs.scipy.org/doc/numpy/reference/generated/numpy.recarray.html
This allows for reach data to be indexed by their variable names, as described in the SFR input instructions.
For more information on Item 2, see the Online Guide to MODFLOW:
http://water.usgs.gov/nrp/gwsoftware/modflow2000/MFDOC/index.html?sfr.htm
In [5]:
rpth = os.path.join('..', 'data', 'sfr_examples', 'test1ss_reach_data.csv')
reach_data = np.genfromtxt(rpth, delimiter=',', names=True)
reach_data
Out[5]:
In [6]:
spth = os.path.join('..', 'data', 'sfr_examples', 'test1ss_segment_data.csv')
ss_segment_data = np.genfromtxt(spth, delimiter=',', names=True)
segment_data = {0: ss_segment_data}
segment_data[0][0:1]['width1']
Out[6]:
dataset 6e is stored in a nested dictionary keyed by stress period and segment,
with a list of the following lists defined for each segment with icalc == 4
FLOWTAB(1) FLOWTAB(2) ... FLOWTAB(NSTRPTS)
DPTHTAB(1) DPTHTAB(2) ... DPTHTAB(NSTRPTS)
WDTHTAB(1) WDTHTAB(2) ... WDTHTAB(NSTRPTS)
In [7]:
channel_flow_data = {0: {1: [[0.5, 1.0, 2.0, 4.0, 7.0, 10.0, 20.0, 30.0, 50.0, 75.0, 100.0],
[0.25, 0.4, 0.55, 0.7, 0.8, 0.9, 1.1, 1.25, 1.4, 1.7, 2.6],
[3.0, 3.5, 4.2, 5.3, 7.0, 8.5, 12.0, 14.0, 17.0, 20.0, 22.0]]}}
dataset 6d is stored in a nested dictionary keyed by stress period and segment,
with a list of the following lists defined for each segment with icalc == 4
FLOWTAB(1) FLOWTAB(2) ... FLOWTAB(NSTRPTS)
DPTHTAB(1) DPTHTAB(2) ... DPTHTAB(NSTRPTS)
WDTHTAB(1) WDTHTAB(2) ... WDTHTAB(NSTRPTS)
In [8]:
channel_geometry_data = {0: {7: [[0.0, 10.0, 80.0, 100.0, 150.0, 170.0, 240.0, 250.0],
[20.0, 13.0, 10.0, 2.0, 0.0, 10.0, 13.0, 20.0]],
8: [[0.0, 10.0, 80.0, 100.0, 150.0, 170.0, 240.0, 250.0],
[25.0, 17.0, 13.0, 4.0, 0.0, 10.0, 16.0, 20.0]]}}
In [9]:
nstrm = len(reach_data) # number of reaches
nss = len(segment_data[0]) # number of segments
nsfrpar = 0 # number of parameters (not supported)
nparseg = 0
const = 1.486 # constant for manning's equation, units of cfs
dleak = 0.0001 # closure tolerance for stream stage computation
istcb1 = 53 # flag for writing SFR output to cell-by-cell budget (on unit 53)
istcb2 = 81 # flag for writing SFR output to text file
dataset_5 = {0: [nss, 0, 0]} # dataset 5 (see online guide)
In [10]:
sfr = flopy.modflow.ModflowSfr2(m, nstrm=nstrm, nss=nss, const=const, dleak=dleak, istcb1=istcb1, istcb2=istcb2,
reach_data=reach_data,
segment_data=segment_data,
channel_geometry_data=channel_geometry_data,
channel_flow_data=channel_flow_data,
dataset_5=dataset_5)
In [11]:
sfr.reach_data[0:1]
Out[11]:
In [12]:
sfr.plot(key='iseg');
In [13]:
chk = sfr.check()
In [14]:
m.external_fnames = [os.path.split(f)[1] for f in m.external_fnames]
m.external_fnames
Out[14]:
In [15]:
m.write_input()
In [16]:
m.run_model()
Out[16]:
In [17]:
sfr_outfile = os.path.join('..', 'data', 'sfr_examples', 'test1ss.flw')
names = ["layer", "row", "column", "segment", "reach", "Qin",
"Qaquifer", "Qout", "Qovr", "Qprecip", "Qet", "stage", "depth", "width", "Cond", "gradient"]
In [18]:
sfrresults = np.genfromtxt(sfr_outfile, skip_header=8, names=names, dtype=None)
sfrresults[0:1]
Out[18]:
In [19]:
import pandas as pd
df = pd.read_csv(sfr_outfile, delim_whitespace=True, skiprows=8, names=names, header=None)
df
Out[19]:
In [20]:
inds = df.segment == 3
ax = df.ix[inds, ['Qin', 'Qaquifer', 'Qout']].plot(x=df.reach[inds])
ax.set_ylabel('Flow, in cubic feet per second')
ax.set_xlabel('SFR reach')
Out[20]:
In [21]:
streambed_top = m.sfr.segment_data[0][m.sfr.segment_data[0].nseg == 3][['elevup', 'elevdn']][0]
streambed_top
Out[21]:
In [22]:
df['model_top'] = m.dis.top.array[df.row.values - 1, df.column.values -1]
fig, ax = plt.subplots()
plt.plot([1, 6], list(streambed_top), label='streambed top')
ax = df.ix[inds, ['stage', 'model_top']].plot(ax=ax, x=df.reach[inds])
ax.set_ylabel('Elevation, in feet')
plt.legend()
Out[22]:
In [23]:
bpth = os.path.join('data', 'test1ss.cbc')
cbbobj = bf.CellBudgetFile(bpth)
cbbobj.list_records()
In [24]:
sfrleak = cbbobj.get_data(text=' STREAM LEAKAGE')[0]
sfrleak[sfrleak == 0] = np.nan # remove zero values
In [25]:
im = plt.imshow(sfrleak[0], interpolation='none', cmap='coolwarm', vmin = -3, vmax=3)
cb = plt.colorbar(im, label='SFR Leakage, in cubic feet per second');
In [26]:
sfrQ = sfrleak[0].copy()
sfrQ[sfrQ == 0] = np.nan
sfrQ[df.row.values-1, df.column.values-1] = df[['Qin', 'Qout']].mean(axis=1).values
im = plt.imshow(sfrQ, interpolation='none')
plt.colorbar(im, label='Streamflow, in cubic feet per second');
In [ ]: