Example 4

Let's use IPython Notebook to download model output from WMT and examine the results.

Set up with pylab magic, plus other global imports:



In [ ]:

    
%pylab inline
import os

Switch to the directory examples/4-WMT to use as our working directory:



In [ ]:

    
os.chdir(os.path.join('..', 'examples', '4-WMT'))
os.getcwd()

The model output is stored as a tar.gz file on the CSDMS website. Here's the name and location of the results of the experiment Multidim Parameter Study MP-2:



In [ ]:

    
run_id = '0b6296d2-ccdd-4717-b347-96be60bfe8e7'
download_file = run_id + '.tar.gz'

You can download and unpack the model output with shell commands. Or, here's a pure Python solution:



In [ ]:

    
def wmt_download_and_unpack(download_file):
    import requests
    import tarfile
    download_url = 'http://csdms.colorado.edu/pub/users/wmt/' + download_file
    r = requests.get(download_url)
    with open(download_file, 'w') as fp:
        fp.write(r.content)
    tar = tarfile.open(download_file)
    tar.extractall()
    tar.close()

Here's the call (it takes a few seconds):



In [ ]:

    
wmt_download_and_unpack(download_file)

Change to the directory containing the unpacked output and get a listing:



In [ ]:

    
os.chdir(run_id)
%ls

This is the standard WMT packaging for a model run. Change to the parameter study directory:



In [ ]:

    
os.chdir('multidim_parameter_study')

Read the Dakota tabular data file:



In [ ]:

    
data = numpy.loadtxt('dakota.dat', skiprows=1, unpack=True, usecols=[0,2,3,4])
data.shape

Reshape the variables from the data file to set up a surface plot:



In [ ]:

    
m = len(set(data[1,]))
n = len(set(data[2,]))
T2 = data[1,].reshape(n,m)
P2 = data[2,].reshape(n,m)
Qs2 = data[3,].reshape(n,m)

Make a surface plot with Axes3D.plot_surface:



In [ ]:

    
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = Axes3D(fig)
ax.plot_surface(T2, P2, Qs2, rstride=2, cstride=2)
ax.view_init(elev=30, azim=-75)
ax.set_xlabel('$T \,(^\circ C)$')
ax.set_ylabel('$P \,(m)$')
ax.set_zlabel('$Q_s \,(kg \, s^{-1})$')