notebook.community

Edit and run 





In [1]:



    


import os
import numpy as np
import scipy.ndimage
import matplotlib.pyplot as plt
import flopy

%matplotlib inline



    











In [2]:



    


fn = os.path.join('data', 'uspb', 'results', 'USPB_capture_fraction_04_01.dat')
cf = np.loadtxt(fn)
print(cf.shape)



    


















    






(110, 80)

















In [3]:



    


cf2 = scipy.ndimage.zoom(cf, 4, order=0)
print(cf2.shape)
#cf2[cf2==0.] = np.nan



    


















    






(440, 320)

















In [4]:



    


c = plt.imshow(cf2, cmap='jet')
plt.colorbar(c)



    


















    
Out[4]:








<matplotlib.colorbar.Colorbar instance at 0x108e2d248>










    
























In [5]:



    


ws = os.path.join('data', 'uspb', 'flopy')
ml = flopy.modflow.Modflow.load('DG.nam', model_ws=ws, verbose=False)



    











In [6]:



    


nlay, nrow, ncol = ml.nlay, ml.dis.nrow, ml.dis.ncol
xmax, ymax = ncol * 250., nrow * 250.



    











In [7]:



    


plt.rcParams.update({'font.size': 6})
fig = plt.figure(figsize=(3.25,4.47))
ax1 = plt.gca()
ax1.set_aspect('equal')
mm1 = flopy.plot.ModelMap(model=ml, layer=4)
plt.xlim(0, xmax)
plt.ylim(0, ymax)
mm1.plot_inactive(color_noflow='0.75')
c = plt.imshow(cf2, cmap='jet', extent=[0, ncol*250., 0, nrow*250.])
cb = plt.colorbar(c, shrink=0.5)
cb.ax.set_ylabel('Layer 4 capture fraction')
mm1.plot_bc(ftype='STR', plotAll=True)
plt.plot([-10000], [-10000], marker='s', ms=10, lw=0.0, mec='0.2', mfc='white', 
         label='Maximum active model extent')
plt.plot([-10000,0], [-10000,0], color='purple', lw=0.75, label='STR reaches (all layers)')
leg = plt.legend(loc='upper left', numpoints=1, prop={'size':6})
leg.draw_frame(False)
plt.xticks([0, 20000, 40000, 60000, 80000])
plt.tight_layout()
plt.savefig(os.path.join('..', 'figs', 'capture_fraction_010y.png'), dpi=300);



    


















    
























In [8]:



    


hedObj = flopy.utils.HeadFile(os.path.join('data', 'uspb', 'flopy', 'DG.hds'), precision='double')
h = hedObj.get_data(kstpkper=(0,0))
cbcObj = flopy.utils.CellBudgetFile(os.path.join('data', 'uspb', 'flopy', 'DG.cbc'), precision='double')
#cbcObj.list_records()
frf = cbcObj.get_data(kstpkper=(0,0), text='FLOW RIGHT FACE')[0]
fff = cbcObj.get_data(kstpkper=(0,0), text='FLOW FRONT FACE')[0]



    











In [9]:



    


cnt = np.arange(1200, 1700, 100)
f, (ax1, ax2) = plt.subplots(1, 2, figsize=(6.75, 4.47))
ax1.set_xlim(0, xmax)
ax1.set_ylim(0, ymax)
ax2.set_xlim(0, xmax)
ax2.set_ylim(0, ymax)
ax1.set_aspect('equal')
ax2.set_aspect('equal')

mm1 = flopy.plot.ModelMap(model=ml, ax=ax1, layer=3)
h1 = mm1.plot_array(h, masked_values=[-888, -999], vmin=1100, vmax=1700)
mm1.plot_inactive(color_noflow='0.75')
mm1.plot_bc(ftype='STR')
q1 = mm1.plot_discharge(frf, fff, istep=5, jstep=5, normalize=True, 
                        color='0.4', scale=70,
                        headwidth=3, headlength=3, headaxislength=3)
c1 = mm1.contour_array(h, masked_values=[-888, -999], colors='black', levels=cnt, 
                       linewidths=0.5)
ax1.clabel(c1, fmt='%.0f', inline_spacing=0.5)

mm2 = flopy.plot.ModelMap(model=ml, ax=ax2, layer=4)
h2 = mm2.plot_array(h, masked_values=[-888, -999], vmin=1100, vmax=1700)
mm2.plot_inactive(color_noflow='0.75')
mm2.plot_bc(ftype='STR')
q2 = mm2.plot_discharge(frf, fff, istep=5, jstep=5, normalize=True, 
                        color='0.4', scale=70,
                        headwidth=3, headlength=3, headaxislength=3)
c2 = mm2.contour_array(h, masked_values=[-888, -999], colors='black', levels=cnt, 
                       linewidths=0.5)
ax2.clabel(c2, fmt='%.0f', inline_spacing=0.5)

ax3 = f.add_axes([0.08, 0.125, 0.01, 0.15])
cb = plt.colorbar(h2, cax=ax3)
cb.ax.set_ylabel('Simulated head, m')

ax1.plot([-10000,0], [-10000,0], color='purple', lw=0.75, label='STR reaches')
ax1.plot([-10000], [-10000], marker='s', ms=10, lw=0.0, mec='black', mfc='None', 
         label='inactive areas')
leg = ax1.legend(loc='upper left', numpoints=1, prop={'size':6})
leg.draw_frame(False)

ax1.text(0.0, 1.01, 'Model layer 4', ha='left', va='bottom',
         transform=ax1.transAxes)
ax2.text(0.98, 0.02, '100 m contour interval', ha='right', va='bottom',
         transform=ax2.transAxes)
ax2.text(0.0, 1.01, 'Model layer 5', ha='left', va='bottom',
         transform=ax2.transAxes)

f.tight_layout()
plt.savefig(os.path.join('..', 'figs', 'uspb_heads.png'), dpi=300);



    


















    






/Users/jdhughes/anaconda/lib/python2.7/site-packages/matplotlib/figure.py:1653: UserWarning: This figure includes Axes that are not compatible with tight_layout, so its results might be incorrect.
  warnings.warn("This figure includes Axes that are not "
/Users/jdhughes/anaconda/lib/python2.7/site-packages/matplotlib/text.py:52: UnicodeWarning: Unicode equal comparison failed to convert both arguments to Unicode - interpreting them as being unequal
  if rotation in ('horizontal', None):
/Users/jdhughes/anaconda/lib/python2.7/site-packages/matplotlib/text.py:54: UnicodeWarning: Unicode equal comparison failed to convert both arguments to Unicode - interpreting them as being unequal
  elif rotation == 'vertical':










    
























In [10]:



    


fn = os.path.join('data', 'uspb', 'results', 'USPB_capture_fraction_04_10.dat')
cf = np.loadtxt(fn)
cf2 = scipy.ndimage.zoom(cf, 4, order=0)



    











In [11]:



    


fig = plt.figure(figsize=(3.25,4.47))
ax1 = plt.gca()
ax1.set_aspect('equal')
mm1 = flopy.plot.ModelMap(model=ml, layer=4)
plt.xlim(0, xmax)
plt.ylim(0, ymax)
mm1.plot_inactive(color_noflow='0.75')
c = plt.imshow(cf2, cmap='jet', extent=[0, ncol*250., 0, nrow*250.])
cb = plt.colorbar(c, shrink=0.5)
cb.ax.set_ylabel('Layer 4 capture fraction')
mm1.plot_bc(ftype='STR', plotAll=True)
plt.plot([-10000,0], [-10000,0], color='purple', lw=0.75, label='STR reaches (all layers)')
plt.plot([-10000], [-10000], marker='s', ms=10, lw=0.0, mec='black', mfc='None', 
         label='Layer 5 inactive area')
leg = plt.legend(loc='upper left', numpoints=1, prop={'size':6})
leg.draw_frame(False)
plt.xticks([0, 20000, 40000, 60000, 80000])
plt.tight_layout()
plt.savefig(os.path.join('..', 'figs', 'capture_fraction_100y.png'), dpi=300);



    


















    
























In [12]:



    


fn = os.path.join('data', 'uspb', 'results', 'USPB_capture_fraction_02_01.dat')
cf_2 = np.loadtxt(fn)
print(cf.shape)



    


















    






(110, 80)

















In [13]:



    


cp_2 = scipy.ndimage.zoom(cf_2, 2, order=0)
print(cf2.shape)



    


















    






(440, 320)

















In [14]:



    


fig = plt.figure(figsize=(3.25,4.47))
ax1 = plt.gca()
ax1.set_aspect('equal')
mm1 = flopy.plot.ModelMap(model=ml, layer=4)
plt.xlim(0, xmax)
plt.ylim(0, ymax)
mm1.plot_inactive(color_noflow='0.75')
c = plt.imshow(cp_2, cmap='jet', extent=[0, ncol*250., 0, nrow*250.])
cb = plt.colorbar(c, shrink=0.5)
cb.ax.set_ylabel('Layer 4 capture fraction')
mm1.plot_bc(ftype='STR', plotAll=True)
plt.plot([-10000,0], [-10000,0], color='purple', lw=0.75, label='STR reaches (all layers)')
plt.plot([-10000], [-10000], marker='s', ms=10, lw=0.0, mec='black', mfc='None', 
         label='Layer 5 inactive area')
leg = plt.legend(loc='upper left', numpoints=1, prop={'size':6})
leg.draw_frame(False)
plt.xticks([0, 20000, 40000, 60000, 80000])
plt.tight_layout()
plt.savefig(os.path.join('..', 'figs', 'capture_fraction_2x_010y.png'), dpi=300);



    


















    
























In [ ]:

Content source: brclark-usgs/flopy

Similar notebooks:

notebook.community | gallery | about