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 commandline python
import numpy as np
import os
import sys
import matplotlib.pyplot as plt
import time
Now that we have some packages to work with, let's see what files we have. These sample files, provided by (more correctly stolen from) David, and are part of a mode-2 simulation.
- rho.15: binary file with the output 15 of the density field
- [x,y,z]grid: binary file with the grids
- matpy: a package with some useful functions (such as darkjet, just for you Marek)
- spinspy: a package with the SPINS readers
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)
In [3]:
help(spy)
Now that we have all of our tools, let's start doing! To begin, let's load our density field.
In [8]:
cd /scratch/orca/scratch/kglamb/ddeepwel/asymkh/
In [9]:
ls rho.*
In [34]:
rho = spy.reader('rho', 30, [0,-1], [0,-1], [0,-1])
What the heck did that just do? Let's break it down and look at the inputs.
'rho': This tells the function which variable to load (in this case, rho)15: This is the index, so we will load rho.15. x,y,z variables done take an index.[0,-1]: The x-range to load. This is shorthand and is equivalent to ':' from MATLAB.[0,-1]: The y-range[0,-1]: The z-rangeHere we have loaded the entire thing. Alternatively, you could have given a single value or a list of indices. That is, giving
0would take the first slice,[0,100]takes the first 100 elements, and[0,2,4,6]takes the first, third, fifth, and, seventh elements.To load an x-z slice we could have done:
rho = spy.reader('rho',15,[0,-1],0,[0,-1])
In [28]:
rho.shape
Out[28]:
Note that this is ordered as (Nx,Ny,Nz). If we had wanted MATLAB style ordering (Ny,Nx,Nz), we could have used the optional argument
ordering = 'matlab'(useful for some plotting). The default isordering = 'natural'(useful because it feels right).If we had high resolution or wanted to load multiply variables into memory simultaneously, we might start running into memory issues. For example, let's look at the size of rho.
In [29]:
(rho.nbytes)/(1024.0**2)
Out[29]:
So rho is taking up about 800 mb of memory. We could easily imagine a higer resolution grid taking up much more, quickly eating up all available memory. Instead, it may be useful to only load a single 1D or 2D slice. This is particularly true with the grids (since we only need a 1D vector to describe each dimension) and when we want to make 2D plots (such as pcolor or contour). Let's look at that now.
To begin, let's load in our grid. We know that they only need to be 1D objects, so they should be pretty small. Remember, for loading
x,y, orzthere is no index value (the 15 from the rho loading call). Notice that the reader automatically squeezes out any singleton dimensions, so our grid vectors are actually vectors!
In [30]:
x = spy.reader('x', [0,-1],0,0)
y = spy.reader('y', 0,[0,-1],0)
z = spy.reader('z', 0,0,[0,-1])
In [31]:
print('The shape of x is {0} and it takes up {1}kb.'.format(x.shape, x.nbytes/(1024.0)))
print('The shape of y is {0} and it takes up {1}kb.'.format(y.shape, y.nbytes/(1024.0)))
print('The shape of z is {0} and it takes up {1}kb.'.format(z.shape, z.nbytes/(1024.0)))
Perfect, now we have our grid vectors and they are fairly small (a few kb). The readers package also provides a function to let you load in all of the grid vectors in one fell swoop.
x,y,z = spy.grid()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 [32]:
data = spy.get_shape()
print(data.Nx,data.Ny,data.Nz)
print('---')
data.display()
We now have our grid vectors and rho 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 plotting functions. For example, here we have a pcolor and contour plot of an x-z slice of rho.
In [35]:
rho = spy.reader('rho', 30, [0,-1], [0,-1], [0,-1])
plt.figure(1)
t0 = time.clock()
plt.contour(x,z,rho[:,192/2,:].transpose());
t1 = time.clock()
print('Time to plot: {0:1.2e}'.format(t1-t0))
plt.title('Mode-2 Waviness');
plt.xlabel('x');
plt.ylabel('z');
In [37]:
rho = spy.reader('rho', 30, [0,-1], [0,-1], [0,-1])
plt.figure(2)
t0 = time.clock()
plt.pcolor(x,z,rho[:,1,:].transpose());
t1 = time.clock()
print('Time to plot: {0:1.2e}'.format(t1-t0))
plt.title('Mode-2 Waviness');
plt.xlabel('x');
plt.ylabel('z');
In [39]:
rho = spy.reader('rho', 30, [0,-1], [0,-1], [0,-1])
plt.figure(3)
t0 = time.clock()
plt.pcolormesh(x,z,rho[:,192/2,:].transpose(), cmap='darkjet')
t1 = time.clock()
print('Time to plot: {0:1.2e}'.format(t1-t0))
plt.title('Mode-2 Waviness')
plt.xlabel('x')
plt.ylabel('z')
plt.axis('tight')
plt.show()
#plt.gca().set_xlim((data.xlim[0],data.xlim[1]))
#plt.gca().set_ylim((data.zlim[0],data.zlim[1]))
We've just seen three examples of 2D plots:
contour,pcolor, andpcolormesh(contourfalso exists). In general, if you'll be working with high resolution rectangular grids, choosepcolormeshoverpcoloras it is cheaper on memory and time. The time taken for each plot was printed when you made the plots and pcolormesh was two orders of magnitude faster!Also, both
pcolorandpcolormeshhave the unfortunate habit of adding extra white space, so you will often need to specify the axis limits.Most MATLAB colour maps are also available in matplotlib, you just need to set the optional argument
cmap='colourmap_name'.
In [21]:
In [62]:
In [62]:
In [62]:
In [62]: