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 [12]:

    
%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 [3]:

    
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 help 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 [4]:

    
help(spy)
spy.set_prefix('Data/2d')









    



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_shape
    grid
    isdim
    make_movie
    reader
    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_shape()
        ## Determine simulation parameters
        ## Purpose:
        ##     If spins.conf exists
        ##         parse spins.conf
        ##     Else
        ##         Raise error
        ##
        ## Usage:
        ##     data = spinspy.get_shape()
        ##
        ##     data.display() prints a summary
        ##                    of known values
        ## ------
    
    grid(type='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(type='full')
        ## ------
    
    reader(var, *args, **kwargs)

DATA
    __all__ = ['get_shape', 'grid', 'reader', 'spinspy_classes', 'get_diag...
    __author__ = 'Ben Storer <bastorer@uwaterloo.ca>'
    __date__ = '29th of April, 2015'

DATE
    29th of April, 2015

AUTHOR
    Ben Storer <bastorer@uwaterloo.ca>


Help on function set_prefix in module spinspy.set_prefix:

set_prefix(prefix)

Reading Full Files

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



In [5]:

    
ls Data/2d









    



q.1         q.12        q.4         q.7         spins.conf
q.10        q.2         q.5         q.8         xgrid
q.11        q.3         q.6         q.9         ygrid

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



In [38]:

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

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

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

0 : This is the index, so we will load h.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.



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 [34]:

    
#x,y = spy.grid()
x = spy.reader('x',0,[0,-1])
y = spy.reader('y',[0,-1],0)
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 [10]:

    
data = spy.get_shape()
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: 1.000e+05
    Bounds of Domain: n/a
    Type: None
  y-Dimension:
    Number of Points: 128
    Length of Domain: 1.000e+05
    Bounds of Domain: n/a
    Type: None
Other parameters:
  beta: 8.42506982132
  dt: 45.0
  rd: 10000000.0
  tmax: 86400.0
  tplot: 6750.0
Settings:

Plotting

We now have our grid vectors and height field 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 [39]:

    
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.23e-01



In [40]:

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









    



Time to plot: 8.42e-01



In [41]:

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









    



Time to plot: 2.18e-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]: