SPINSpy Tutorial (2D):

This document will (hopefully) show you everything that you need to know to be able to process your SPINS outputs using python (why python? Because it's not MATLAB). This tutorial will assume that you have a basic understanding of python syntax.

To begin, we import some useful packages. Base python is fairly bare-bones, but you'll find that there's a package for almost anything that you want.



In [1]:

    
%matplotlib inline 
# Tells the system to plot in-line, only necessary for iPython notebooks,
# not regular command-line python
import numpy as np
import os
import sys
import matplotlib.pyplot as plt
import time



In [2]:

    
# Now that we have our packages, we need data. The file 'make_2d_data.py' will
# generate a sample data set. Let's run that now (may take a minute)

execfile('make_2d_data.py')









    



Creating 2D data...
Done!

If we want to use our shiny python scripts, we'll need to import them too.



In [2]:

    
import spinspy as spy
import matpy as mp

If we want a quick man-page style summary, we can call help(spy). We can also call help(spy.<func>) for more information on the function <func>. Since the data that we're looking at is in a different directory, let's specify that now.

Note: This means that you don't need to be in the same directory as your data.



In [5]:

    
help(spy.set_path)
print('=======')
print('=======')
print('=======')
help(spy)
spy.set_path('Data/2d')









    



Help on function set_path in module spinspy.set_path:

set_path(path)
    ## Specify the baseline path for the SPINSpy
    ## functions. In particular, data will be read
    ## from path/u.0, etc.

=======
=======
=======
Help on package spinspy:

NAME
    spinspy

FILE
    /u/bastorer/Documents/PyLab/SPINSpy/spinspy/__init__.py

DESCRIPTION
    # SPINSPY
    #   This module contains functions that are
    #   designed to handle SPINS-type outputs.
    #   The provided functions are listed below,
    #   along with basic usage information.

PACKAGE CONTENTS
    get_diagnostics
    get_grid
    get_gridparams
    get_paramgrid
    get_params
    get_shape
    grid
    isdim
    make_movie
    nearestindex
    reader
    set_path
    set_prefix
    spinspy_classes

FUNCTIONS
    get_diagnostics(fp='diagnostics.txt')
        ## ------------------
        ## Purpose: Parse diagnostics.txt and return a
        ##          class that contains each diagnostic
        ##          variable.
        ## Inputs:
        ##    fp (optional): Specifies the filename of
        ##                   of the diagnostics file.
        ## Usage: 
        ##    diag = spinspy.get_diagnostics()
    
    get_grid(style='vector')
        ## Read in the grid as either vectors (1)
        ## or full matrices (2).
        ## (1) x,y[,z] = spinspy.grid()
        ## (2) x,y[,z] = spinspy.grid(style='full')
        ## ------
    
    get_gridparams(style='vector')
        ## Determine simulation parameters and make grid
        ## Purpose:
        ##     - Read spins.conf if it exists
        ##     - Read the grid
        ##     - Check and add other grid parameters (such as grid spacing) into parameters object
        ##
        ## Usage:
        ##     x,y[,z],params = spinspy.get_paramgrid()
        ## ------
    
    get_params()
        ## Determine simulation parameters
        ## Purpose:
        ##     If spins.conf exists
        ##         parse spins.conf
        ##     Else
        ##         Raise error
        ##
        ## Usage:
        ##     data = spinspy.get_params()
        ##
        ##     data.display() prints a summary
        ##                    of known values
        ## ------
    
    nearestindex(list, value)
        ## Find the index of list that is closest to value
        ## assumes list is monotonic and is a vector (though it will work otherwise)
        ## ------
    
    reader(var, *args, **kwargs)
        ## Load SPINS outputs into numpy arrays.
        ## Usage: field = spinspy.reader(var, *args, **kwargs)
        ## Inputs:
        ##    var: string indicating the field to be read, 'rho', 'x', etc.
        ##    args: Required and ordered inputs.
        ##        index: integer indicating output for non-grid fields
        ##            i.e. 15 to load field.15
        ##        slicing: how to slice, one for each dimension
        ##    kwargs: Optional keyword arguments
        ##        type: 'ndarray' (default) or 'memmap'
        ##              determines the output type
        ##        ordering: 'natural' (default) or 'matlab'
        ##              determines dimension ordering (x,y,z) vs (y,x,z)
        ##
        ## Examples:
        ## 1) To load rho.15 for a full 3D grid
        ##      rho = spinspy.reader('rho',15,[0,-1],[0,-1],[0,-1])
        ## 2) To load rho.15 on a 3D grid but at first z
        ##      rho = spinspy.reader('rho',15,[0,-1],[0,-1],0)
        ## 3) To load x (grid) as a vector
        ##        x = spinspy.reader('x',[0,-1],0,0)
        ## 4) To load rho.15 on a 2D grid, but in matlab ordering
        ##      rho = spinspy.reader('rho',15,[0,-1],[0,-1],ordering='matlab')
        ## ------
    
    set_path(path)
        ## Specify the baseline path for the SPINSpy
        ## functions. In particular, data will be read
        ## from path/u.0, etc.

DATA
    __all__ = ['spinspy_classes', 'get_params', 'get_grid', 'get_gridparam...
    __author__ = ['Ben Storer <bastorer@uwaterloo.ca>', 'David Deepwell, <...
    __date__ = '29th of July, 2015'

DATE
    29th of July, 2015

AUTHOR
    ['Ben Storer <bastorer@uwaterloo.ca>', 'David Deepwell, <ddeepwel@uwaterloo.ca>']

Reading Full Files

Now that we have all of our tools, let's start doing! Let's see what we have first.



In [6]:

    
ls Data/2d









    



q.0         q.14        q.2         q.25        q.30        q.7
q.1         q.15        q.20        q.26        q.31        q.8
q.10        q.16        q.21        q.27        q.32        q.9
q.11        q.17        q.22        q.28        q.4         spins.conf
q.12        q.18        q.23        q.29        q.5         xgrid
q.13        q.19        q.24        q.3         q.6         ygrid

We have q (potential vorticity, this is a QG example). Let's read and plot our initial q field.



In [7]:

    
q = spy.reader('q', 0, [0,-1], [0,-1])

What did that just do? Let's break it down and look at the inputs.

'q' : This tells the function which variable to load (in this case, q)

0 : This is the index, so we will load q.0.

[0,-1] : The x-range to load. This is shorthand and is equivalent to ':' from MATLAB.

[0,-1] : The y-range

Here we have loaded the entire thing. Alternatively, you could have given a single value or a list of indices. That is, giving 0 would take the first slice, [0,100] takes the first 100 elements, and [0,2,4,6] takes the first, third, fifth, and, seventh elements.

Lazy loading had also been implemented so that q=spy.reader('q',0) would have also worked.



In [8]:

    
q.shape









    Out[8]:





(128, 128)

Note that this is ordered as (Nx,Ny). If we had wanted MATLAB style ordering (Ny,Nx), we could have used the optional argument ordering = 'matlab' (useful for some plotting). The default is ordering = 'natural' (useful because it feels right).

Loading Grids and Simulation Information

To begin, let's load in our grid. Since we aren't dealing with mapped grids, we only need vectors (the default behaviour).



In [11]:

    
x,y = spy.get_grid()
print('The shape of x is {0} and the shape of y is {1}'.format(x.shape,y.shape))









    



The shape of x is (128,) and the shape of y is (128,)

Perfect, now we have our grid vectors. We can also load in grid information, such as the domain size and limits. This is stored in a class that has a method called display. These are illustrated below.



In [12]:

    
data = spy.get_params()
print(data.Nx,data.Ny,data.Nz)
print('---')
data.display()









    



(128, 128, 1)
---
2-dimensional simulation:
  x-Dimension:
    Number of Points: 128
    Length of Domain: 1e+05
    Bounds of Domain: n/a
    Type: None
  y-Dimension:
    Number of Points: 128
    Length of Domain: 1e+05
    Bounds of Domain: n/a
    Type: None
Other parameters:
  beta: 8.42507
  dt: 45
  rd: 1e+07
  tmax: 86400
  tplot: 2700
Settings:

Plotting

We now have our grid vectors and potential vorticity loaded into memory, so it's time to try some plotting! A standard plotting library in python is matplotlib, which replicates many of MATLAB's 2D plotting functions. For example, here we have a pcolor and contour plot.



In [16]:

    
q = spy.reader('q', 10)



In [17]:

    
plt.figure(1)
t0 = time.clock()
plt.contour(x,y,q.T) # Transpose for plotting order
t1 = time.clock()
print('Time to plot: {0:1.2e}'.format(t1-t0))
plt.title('q');
plt.xlabel('x');
plt.ylabel('y');









    



Time to plot: 1.13e-01



In [20]:

    
plt.figure(2)
t0 = time.clock()
plt.pcolor(x,y,q.transpose()); # Regrettably, the default colormap is still jet...
t1 = time.clock()
print('Time to plot: {0:1.2e}'.format(t1-t0))
plt.colorbar()
plt.title('Height Field');
plt.xlabel('x');
plt.ylabel('y');









    



Time to plot: 1.09e+00



In [21]:

    
plt.figure(3)
t0 = time.clock()
plt.pcolormesh(x,y,q.T, cmap='darkjet') # So we made darkjet!
t1 = time.clock()
print('Time to plot: {0:1.2e}'.format(t1-t0))
plt.colorbar()
plt.title('Height Field')
plt.xlabel('x')
plt.ylabel('y')
plt.axis('tight')
plt.show()









    



Time to plot: 2.06e-01

We've just seen three examples of 2D plots: contour, pcolor, and pcolormesh (contourf also exists). In general, if you'll be working with high resolution rectangular grids, choose pcolormesh over pcolor as it is cheaper on memory and time. The time taken for each plot was printed when you made the plots and pcolormesh was significantly faster! (Note: pcolormesh will not work for mapped grids)

Also, both pcolor and pcolormesh have the unfortunate habit of adding extra white space, so you will often need to specify the axis limits or call plt.axis('tight')

Most MATLAB colour maps are also available in matplotlib, you just need to set the optional argument cmap='colourmap_name'. The package matpy provides darkjet.



In [62]: