This notebook demonstrates the simple array and results plotting capabilities of flopy. It demonstrates these capabilities by loading and running an existing model, and then showing how the .plot() method can be used to make simple plots of the model data and model results.
In [1]:
%matplotlib inline
from IPython.display import Image
import sys
import os
import platform
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import flopy
print(sys.version)
print('numpy version: {}'.format(np.__version__))
print('matplotlib version: {}'.format(mpl.__version__))
print('flopy version: {}'.format(flopy.__version__))
In [2]:
#Set name of MODFLOW exe
# assumes executable is in users path statement
version = 'mf2005'
exe_name = 'mf2005'
if platform.system() == 'Windows':
exe_name = 'mf2005.exe'
mfexe = exe_name
#Set the paths
loadpth = os.path.join('..', 'data', 'secp')
modelpth = os.path.join('data')
#make sure modelpth directory exists
if not os.path.exists(modelpth):
os.makedirs(modelpth)
files = ['secp.hds']
A model called the "Southeast Coastal Plain Model" is located in the loadpth folder. In the following code block, we load that model, then change into a new workspace (modelpth) where we recreate and run the model. For this to work properly, the MODFLOW-2005 executable (mf2005) must be in the path. We verify that it worked correctly by checking for the presence of secp.hds.
In [3]:
ml = flopy.modflow.Modflow.load('secp.nam', model_ws=loadpth,
exe_name=exe_name, version=version)
ml.change_model_ws(new_pth=modelpth)
ml.write_input()
success, buff = ml.run_model(silent=True)
if not success:
print ('Something bad happened.')
# confirm that the model files have been created
for f in files:
if os.path.isfile(os.path.join(modelpth, f)):
msg = 'Output file located: {}'.format(f)
print (msg)
else:
errmsg = 'Error. Output file cannot be found: {}'.format(f)
print (errmsg)
In [4]:
ml.dis.top.plot()
Out[4]:
As you can see, the .plot() methods returns a matplotlib.pyplot axis object, which can be used to add additional data to the figure. Below we will add black contours to the axis returned in the first line.
In [5]:
ax = ml.dis.top.plot()
ml.dis.top.plot(axes=ax, contour=True, pcolor=False)
Out[5]:
You will notice that we passed several keywords in the second line. There are a number of keywords that can be passed to the .plot() method to control plotting. Available keywords are:
axes - if you already have plot axes you can pass them to the methodpcolor - turns pcolor on if pcolor=True or off if pcolor=False, default is pcolor=Truecolorbar - turns on colorbar if colorbar=True or off if colorbar=False, default is colorbar=False and is only used if pcolor=Trueinactive - turns on a black inactive cell overlay if inactive=True or turns off the inactive cell overlay if inactive=False, default is inactive=Truecontour - turns on contours if contour=True or off if contour=False, default is contour=Falseclabel - turns on contour labels if clabel=True or off if clabel=False, default is clabel=False and is only used if contour=Truegrid - turns on model grid if grid=True or off if grid=False, default is grid=Falsemasked_values - list with unique values to be excluded from the plot (for example, HNOFLO)mflay - for three-dimensional data (for example layer bottoms or simulated heads) mflay can be used to plot data for a single layer - note mflay is zero-basedkper - for transient two-dimensional data (for example recharge package data) kper can be used to plot data for a single stress period - note kper is zero-basedfilename_base - a base file name that will be used to automatically generate file names for two-dimensional, three-dimensional, and transient two-dimensional data, default is filename_base=Nonefile_extension - valid matplotlib file extension, default is png and is only used if filename_base is specifiedmatplotlib.pyplot keywords are also accepted The previous code block is recreated in a single line using keywords in the code block below.
In [6]:
ml.dis.top.plot(contour=True)
Out[6]:
We can save the same image to a file.
In [7]:
fname = os.path.join(modelpth, 'secp')
ml.dis.top.plot(contour=True, filename_base=fname)
The image file that was just created is shown below
In [8]:
fname = os.path.join(modelpth, 'secp_model_top.png')
Image(filename=fname)
Out[8]:
Three-dimensional data can be plotted using the .plot() method. User's do not actually need to know that the data are two- or three-dimensional. The .plot() method is attached to the two- and three-dimensional data objects so it knows how to process the model data. Examples of three-dimensional data are horizontal hydraulic conductivity (hk), layer bottoms (botm), specific yield (sy), etc.
Here we plot the horizontal hydraulic conductivity for each layer. We are also masking areas where the horizontal hydraulic conductivity is zero and adding a color bar.
In [9]:
ml.lpf.hk.plot(masked_values=[0.], colorbar=True)
Out[9]:
In addition to the plots of horizontal hydraulic conductivity you can see that the .plot() method returned a list containing 16 axis objects (one for each layer).
If the mflay= keyword is provided to the plot.() method then data for an individual layer is plotted. Remember mflay is zero-based.
Here we plot the horizontal hydraulic conductivity for layer 11 (mflay=10).
In [10]:
ml.lpf.hk.plot(mflay=10, masked_values=[0.], colorbar=True)
Out[10]:
Transient two-dimensional data can be plotted using the .plot() method. User's do not actually need to know that the data are two- or three-dimensional. The .plot() method is attached to the two- and three-dimensional, and transient two-dimensional data objects so it knows how to process the model data. Examples of transient two-dimensional data are recharge rates (rch.rech) and evapotranspiration rates (evt.evtr).
Here we plot recharge rates for all six stress periods in the model. We are also masking areas where the recharge rate is zero and adding a color bar.
In [11]:
ml.rch.rech.plot(kper='all', masked_values=[0.], colorbar=True)
Out[11]:
In addition to the plots of recharge rates you can see that the .plot() method returned a list containing 6 axis objects (one for each stress period).
If the kper= keyword is provided to the plot.() method then data for an individual stress period is plotted. Remember kper is zero-based.
Here we plot the recharge rate for stress period 6 (kper=5).
In [12]:
ml.rch.rech.plot(kper=5, masked_values=[0.], colorbar=True)
Out[12]:
We can also save the image to a file by provided the filename_base keyword with an appropriate base file name.
In [13]:
fr = os.path.join(modelpth, 'secp_rch')
ml.rch.rech.plot(kper=5, masked_values=[0.], colorbar=True,
filename_base=fr)
Out[13]:
If the kper keyword is not provided images are saved for each stress period in the model.
The image file that was just created of recharge rates for stress period 6 is is shown below.
In [14]:
fname = os.path.join(modelpth, 'secp_rch_00006.png')
Image(filename=fname)
Out[14]:
Simulated model results can be plotted using the .plot() method.
First we create an instance of the HeadFile class with the simulated head file (secp.hds) and extract the simulation times available in the binary head file using the .get_times() method. Here we plot last simulated heads in the binary heads file (totim=times[-1]). We are also masking cells having the HDRY (-9999.0) value and adding a color bar.
In [15]:
print(modelpth, files)
fname = os.path.join(modelpth, files[0])
hdobj = flopy.utils.HeadFile(fname, model=ml)
times = hdobj.get_times()
head = hdobj.plot(totim=times[-1],
masked_values=[-9999.0],
colorbar=True)
In [16]:
head = hdobj.plot(totim=times[1], mflay=1,
masked_values=[-9999.0],
colorbar=True,
contour=True, colors='black')
We can also save the plots of the head results for every layer (or a single layer) to a file by provided the filename_base keyword with an appropriate base file name.
In [17]:
fh = os.path.join(modelpth, 'secp_head')
head = hdobj.plot(totim=times[-1],
masked_values=[-9999.0],
colorbar=True,
contour=True, colors='black',
filename_base=fh)
The image file that was just created of the simulated head for layer 2 for stress period 6 is is shown below.
In [18]:
fname = os.path.join(modelpth, 'secp_head_Layer2.png')
Image(filename=fname)
Out[18]:
matplotlib.pyplot keywords to .plot() methodsWe can also pass matplotlib.pyplot keywords to .plot() methods attached to the model input data arrays. For example you can pass a matplotlib colormap (cmap=) keyword to the .plot() method to plot contours of simulated heads over a color flood of hk. We can also use the norm=LogNorm() keyword to use a log color scale when plotting hydraulic conductivity.
Available matplotlib colormaps can be found at http://matplotlib.org/examples/color/colormaps_reference.html
In [19]:
from matplotlib.colors import LogNorm
ax = ml.lpf.hk.plot(mflay=2, cmap='GnBu', norm=LogNorm(), colorbar=True)
t = hdobj.plot(axes=ax, mflay=2, masked_values=[-9999.0],
pcolor=False, contour=True,
colors='black')
The input data for a model or an individual package can also be plotted using the .plot() method. The .plot() methods attached to a model or an individual package are meant to provide a method to quickly evaluate model or package input. As a result, there is limited ability to customize the plots. Example of using the .plot() method with a model or and individual packages is demonstrated below.
All input data for a package can be plotted using the .plot() method. Below all of the data for the lpf package is plotted.
In [20]:
ml.dis.plot()
Out[20]:
In [21]:
ml.lpf.plot(mflay=11)
Out[21]:
In [22]:
ap = ml.plot()
In [23]:
ap = ml.plot(mflay=11)
This notebook demonstrates some of the simple plotting functionality available with flopy. Although not described here, the plotting functionality tries to be general by passing keyword arguments passed to the plot() and plot_data() methods down into the matplot.pyplot routines that do the actual plotting. For those looking to customize these plots, it may be necessary to search for the available keywords. The plot() method return the matplot.pyplot axis objects that are created (or passed). These axes objects can be used to plot additional data (except when plots are saved as image files).
Hope this gets you started!
In [ ]: