Notebook

VISUALIZATION TECHNIQUES IN PYTHON¶

Scientific computing with python seminar¶

by Mauricio J. Del Razo S. 2014

Although visualization might seem as a secondary issue in terms of scientific computing, the reality is that visualization can be fundamental for understanding a problem. When giant data sets are analyzed, appropriate visualization might uncover results thet would remain unknown otherwise. Additionally, visualization is how you present your results and how people understand them. Good visualization might improve your chances of your work being published understood and cited.

The visualization covered in this lesson might be still introductory, though it will introduce some very useful methods and tools that can improve the presentation of your work, and it will point out a good directions to go if you want to go deeper into the rabbit hole.

1. MATPLOTLIB BASICS¶

Here we will explore some of the basic visualization techniques using the matplotlib library. As an extended reference, check Jake's notebook on Advanced Matplolib.

In order to import matplolib library within the ipython notebook; it is only required to add the line: "%pylab inline" in the begginign in the code:

In [1]:

%pylab inline

Populating the interactive namespace from numpy and matplotlib

Or equivalently in a python script one could type at the beggining of the document:

In [2]:

import numpy as np
import matplotlib.pyplot as plt
from pylab import *

We can now produce a simple figure

In [3]:

# Create values to plot
x = linspace(-5,5,50)
y = np.sin(x)
# Plot figure
plt.plot(x,sin(x))
plt.show() # Not neccesary on ipython notebooks, but neccesary on python scripts

Or an even nicer figure by adding titles and some nicer formatting, which is fundamental for papers

In [4]:

# Create nicer plot, legends are self explanatory
plt.plot(x,sin(x), linewidth = 4)
plt.axis([-5,5,-1.3,1.3])
plt.xlabel('time', fontsize = 23)
plt.ylabel('amplitude', fontsize = 23)
plt.title('NOT REALLY MY FIRST PLOT', fontsize = 25)
xticks = linspace(-5,5,5)
yticks = linspace(-1,1,5)
plt.xticks(xticks, fontsize = 15);
plt.yticks(yticks, fontsize = 15);

However, now we know object oriented programming, so lets use it!. It will allows to do more complicated and better organized scripts

In [5]:

# Create a figure object
fig = plt.figure(figsize=(10,7))

# We can add several subplots to the figure with add_subplot(rows,columns,plotnumber)
ax1 = fig.add_subplot(2,1,1) #
ax1.set_title('My top subplot', fontsize = 20, color = 'b', fontweight = 'bold')
ax1.plot(x,np.sin(x),linewidth = 3, color = 'b')
ax1.plot(x,np.cos(x),linewidth = 3, color = 'r', linestyle = '--')
ax1.axis([-5,5,-1.3,1.3])
ax1.set_xticks([])
ax1.yaxis.set_tick_params(labelsize=25)

# Middle subplot
ax2 = fig.add_subplot(2,1,2)
ax2.set_title('My bottom subplot', fontsize = 22, color = 'r', fontweight = 'bold')
ax2.plot(x,np.cos(x),linewidth = 3, color = 'r')
ax2.fill_between(x,-2, np.cos(x), facecolor='blue', alpha=0.2)
ax2.axis([-5,5,-1.3,1.3])
ax2.set_xticks(xticks);
ax2.set_yticks(yticks);
ax2.xaxis.set_tick_params(labelsize=25)
ax2.yaxis.set_tick_params(labelsize=25)

Now lets do a more realistic example of what would be a good way to present your results in a talk or in a paper:

Remember to always create a script (or an ipython notebook) to create your figures, so if you need to change anything after

a review you can reproduce the exact same figure easily.

If you are planning to create different plots for different parameters is a good idea to create a function

In [6]:

def plot_figure1(time=0,xcut=6,ycut=-5):
    # Lets create some data to visualize
    nx, ny = (150,150)
    xlim = 10
    ylim = 15
    x = np.linspace(-xlim,xlim,nx)
    y = np.linspace(-ylim,ylim,ny)
    # Create a mesh grid of the two arrays x and y
    xm, ym = meshgrid(x, y)
    # Compute a value on that mesh
    z = sin(np.sqrt(xm**2 + 3*ym**2) - time)/(xm**2 + ym**2)**(0.25)
    
    # Transform x_cut and y_cut
    x_cut = int(nx*xcut/(2*xlim) + nx/2.0)
    y_cut = int(ny*ycut/(2*ylim) + ny/2.0)

    # Lets create a plotting grid with gridspec
    # Create a 3 x 3 plotting grid object with horizontal and vertical spacing wspace and hspace
    gs = plt.GridSpec(2, 2, wspace=0.2, hspace=0.2) 
    # Create a figure
    fig = plt.figure(figsize=(16, 6))

    # SUBFIGURE 1
    # Create subfigure 1 (takes over two rows (0 to 1) and column 0 of the grid)
    ax1 = fig.add_subplot(gs[:2, 0])
    # Plot using pcolor and arrange colorbar
    colormap = 'Blues'#'YlGnBu' or 'afmhot' or 'RdBu' or many others
    s1 = ax1.pcolor(xm,ym,z, cmap = colormap, vmin=-0.2, vmax=0.7) # vmin and vmax are min and max values of colormap
    cbar = fig.colorbar(s1) # Show colorbar
    cbar_ticks = np.linspace(-0.2,0.7,5) 
    cbar.set_ticks(cbar_ticks) # adjust colorbar ticks
    # Add 10 black contours (dashed means negative)
    ax1.contour(xm,ym,z, 10, colors = 'k')
    # Add xcut and ycut slice line
    ax1.plot([-xlim, xlim], [xcut, xcut], color = 'r', linewidth = 2)
    ax1.plot([ycut, ycut], [-ylim, ylim], color = 'g', linewidth = 2)
    # Arrange axis and labels
    ax1.axis([-xlim,xlim,-ylim,ylim])
    ax1.set_xlabel('x', fontsize = 20)
    ax1.set_ylabel('y', fontsize = 20)
    ticks = np.linspace(-xlim,xlim,5)
    ax1.set_xticks(ticks);
    ax1.set_yticks(ticks);
    ax1.xaxis.set_tick_params(labelsize=20)
    ax1.yaxis.set_tick_params(labelsize=20)
    ax1.set_title('Sine wave amplitude', fontsize = 20)

    # SUBFIGURE 2
    # Create subfigure 2 (takes over row 0 and column 1 of the grid)
    ax2 = fig.add_subplot(gs[0, 1])
    ax2.plot(xm[0,:],z[x_cut,:],linewidth = 3, color = 'r')
    ax2.fill_between(xm[0,:],-2, z[x_cut,:], facecolor='red', alpha=0.1)
    ax2.set_title('Horizontal slice', fontsize = 20)
    ax2.axis([-xlim,xlim,-1.3,1.3])
    ax1.set_xlabel('x', fontsize = 20)
    ax2.xaxis.set_tick_params(labelsize=15)
    ax2.yaxis.set_tick_params(labelsize=15)

    # SUBFIGURE 3
    # Create subfigure 3 (takes over row 1 and column 1 of the grid)    
    ax3 = fig.add_subplot(gs[1, 1])
    ax3.plot(ym[:,0],z[:,y_cut],linewidth = 3, color = 'g')
    ax3.fill_between(ym[:,0],-2, z[:,y_cut], facecolor='green', alpha=0.1)
    ax3.set_title('Vertical slice', fontsize = 20)
    ax3.axis([-ylim,ylim,-1.3,1.3])
    ax3.set_xlabel('y', fontsize = 20)
    ax3.xaxis.set_tick_params(labelsize=15)
    ax3.yaxis.set_tick_params(labelsize=15)
    
from IPython.html.widgets import interact
interact(plot_figure1, time=(0,5,0.1), xcut = (-10,10,1), ycut = (-10,10,1));

You can always use this template as a reference to produce nice pictures for your works. Also, remeber you don't need the ipython notebook to produce these plots. You can use the function we just created directly in a python script (just remeber importing the libraries at the beggining of the script as I did in the beggining of this document). For instance, if you add these commands at the end of your python script, you'll obtain the corresponding figure

In [7]:

plot_figure1(0,-7,5)
plt.show()

2. STEPPING OUT INTO THE THIRD DIMENSION¶

There are several very good tools for three dimensional visualization. Some of these packagaes are big independent software packages and some of them you can interface with python. The bigger packages might be an overkill for simple plots. Some examples are:

MayaVi can easily be used with python, and it is quite powerful easy and fairly easy to use. VisIt and ParaView are more hardcore, so you can try them out if you want to do something more heavy duty. In this lesson, we will employ mplot3D, which is native to python. Lets start with a couple of examples from the mplot3D documentation website:

Plotting a parametric curve¶

In [8]:

from mpl_toolkits.mplot3d import Axes3D
# Create figure
fig = plt.figure()
ax = fig.gca(projection='3d')
# Create parametric curve
theta = np.linspace(-4 * np.pi, 4 * np.pi, 100)
z = np.linspace(-2, 2, 100)
r = z**2 + 1
x = r * np.sin(theta)
y = r * np.cos(theta)
ax.plot(x, y, z, label='parametric curve', linewidth = 3)
ax.legend()

plt.show()

Surface plotting with contour projections¶

In [9]:

from mpl_toolkits.mplot3d import axes3d
from matplotlib import cm

# Create figure and get data
fig = plt.figure()
ax = fig.gca(projection='3d')
X, Y, Z = axes3d.get_test_data(0.05)

# Plot surface and contours
ax.plot_surface(X, Y, Z, rstride=8, cstride=8, alpha=0.3)
cset = ax.contourf(X, Y, Z, zdir='z', offset=-100, cmap=cm.coolwarm)
cset = ax.contour(X, Y, Z, zdir='x', offset=-40, cmap=cm.coolwarm)
cset = ax.contour(X, Y, Z, zdir='y', offset=40, cmap=cm.coolwarm)

ax.set_xlabel('X')
ax.set_xlim(-40, 40)
ax.set_ylabel('Y')
ax.set_ylim(-40, 40)
ax.set_zlabel('Z')
ax.set_zlim(-100, 100)

Out[9]:

(-100, 100)

A 3D function of time and viewpoint¶

Now we will create a function that plots a surface as a function of time, and also as a function of the camera viewpoint. We will use this function to produce our animations, so the camera can rotate while times passes by.

In [10]:

# Import 3D libraries
from mpl_toolkits.mplot3d import Axes3D

# It's always good to create a function for your plots
def plot3dsine(theta, phi, time):
    
    # Create figure
    fig = plt.figure(figsize=(12, 6))
    ax = fig.gca(projection='3d')

    # Create data to plot (same as before)
    nx, ny = (150,150)
    xlim = 10
    ylim = 10
    x = np.linspace(-xlim,xlim,nx)
    y = np.linspace(-ylim,ylim,ny)
    # Create a mesh grid of the two arrays x and y
    xm, ym = meshgrid(x, y)
    # Compute a value on that mesh
    z = sin(np.sqrt(xm**2 + 3*ym**2) - time)/(xm**2 + ym**2)**(0.25)

    # Produce 3d plot, set axis and colorbar
    specs = {'cmap' : 'seismic', 'vmin' : -0.6, 'vmax' : 1.0}
    surf = ax.plot_surface(xm, ym, z, rstride = 1, cstride = 2, linewidth=0, **specs)
    cset = ax.contourf(xm, ym, z, zdir='z', offset=3, alpha = 0.2, **specs)
    cset = ax.contour(xm, ym, z, zdir='y', offset=-15, **specs)
    cset = ax.contour(xm, ym, z, zdir='x', offset=-15, **specs)

    ax.set_zlim(-1.01, 1.01)
    ax.set_axis_off() # Remove background axis and grid

    ax.view_init(elev=phi, azim=theta)
    fig.colorbar(surf, shrink=0.5, aspect=5)
    

plot3dsine(45,290,0)

3. HOW TO PRODUCE AN ANIMATION¶

There are endless ways of producing animations; we will explore two of them.

Output a set of images from ipython and use another software to write them into a video format
Use Jake's package JSAnimation , to produce html files where you can reproduce the movie.

Produce video file from a set of images¶

For the first case, let make python output a set of images:

In [11]:

!mkdir animation_frames

mkdir: animation_frames: File exists

In [12]:

# Output a set of images using plot3dsine()
nframes = 10 # Use 100 to produce bigger animation
print('Frame number:')
for i in xrange(0,nframes,1):
    # 3D plot parameters as function of iterator i
    theta = 45 - 0.3*i
    phi = 290 + 0.4*i
    time = i/25.0
    # Call plotting function and save figure into file
    plot3dsine(theta,phi,time);
    plt.savefig('animation_frames/frame{0:07d}.png'.format(i+1));
    plt.close()
    print i,

Frame number:
0 1 2 3 4 5 6 7 8 9

If you're using linux, now we can open a terminal and use a third party program like: avconv, ffmpeg or mencoder (maybe even on a Mac). On an Mac you can also even use iMovie, on windows ... Good luck! Lets do an example in Ubuntu, on a terminal do the following:

Install avconv: sudo apt-get install avconv
avconv -r 10 -i frame%07d.png -b:v 1000k amovie.mp4

where -r 10 indicates the frames per second and 1000k the bitrate (quality) of the video.

Produce an animation html file with JSAnimation¶

Now we will use a tool created by Jake Vanderplas called JSAnimation, whch stands for JavaScript animation. It allow us to create an animation within an html file by employing javascript. The advantage of this tool is that we do not need to dive into the JavaScript code, we can employ python uniquely to create the animations. In order to use JSAnimation, we first need to download it and install it:

Go to the JSAnimation GitHub webapge and clone the repository in your favorite local folder by typing: git clone https://github.com/jakevdp/JSAnimation.git
Go into the JSAnimation folder and install it by typing: python setup.py install

You might need administrator rights, so you might need to type sudo before: sudo python setup.py install

Here is an example, where we will use the images we already used to produce the video:

In [13]:

# Import JSAnimation 
from JSAnimation import IPython_display, HTMLWriter
from matplotlib import animation
import Image
import glob

# Choose location of files to use
filenames=sorted(glob.glob('animation_frames/frame*.png'))

#Create figure
fig = plt.figure(figsize=(10,5))
ax = fig.add_axes([0, 0, 1, 1])
ax.axis('off')
im = plt.imshow(Image.open(filenames[0]))

# Create initializatio function
def init():
    im.set_data(Image.open(filenames[0]))
    return im,

# Create animation function
def animate(i):
    image=Image.open(filenames[i])
    im.set_data(image)
    return im,

anim = animation.FuncAnimation(fig, animate, init_func=init,
                               frames=nframes, interval=20, blit=True)

#set embed_frames=True to embed base64-encoded frames directly in the HTML
anim.save('animation_frames/animation.html', writer=HTMLWriter(embed_frames=True))

# Show animation in ipython notebook
anim

Out[13]:

Note you can decide if you embed the files into your html file by setting embed_frames to True or False. If False is chosen a new folder with all the image files will be created.

Display an animation in the ipython notebook¶

We can also create an animation within the ipython notebook without even creating the image files in the first place

In [14]:

# JSAnimation import available at https://github.com/jakevdp/JSAnimation
from JSAnimation import IPython_display
from matplotlib import animation

# create a simple animation
fig = plt.figure()
ax = plt.axes(xlim=(0, 10), ylim=(-2, 2))
line, = ax.plot([], [], lw=3)

x = np.linspace(0, 10, 1000)

def init():
    line.set_data([], [])
    return line,

def animate(i):
    line.set_data(x, np.cos(i * 0.02 * np.pi) * np.sin(x - i * 0.02 * np.pi))
    return line,

animation.FuncAnimation(fig, animate, init_func=init,
                        frames=100, interval=20, blit=True)

Out[14]:

If we want to save it, we only need to do:

In [15]:

# Save the animation into an html file
anim = animation.FuncAnimation(fig, animate, init_func=init,
                        frames=100, interval=20, blit=True)
anim.save('animation_frames/animation_cos.html', writer=HTMLWriter(embed_frames=True))

4. ASSIGNMENT: PLOTTING AND ANIMATING THE DRUMHEAD VIBRATIONAL MODES¶

Whenever you hit a drum, the drumhead vibrates following the two dimensional wave equation in polar coordinates,

$$ \frac{\partial u(r,\theta,t)}{\partial t} = c^2 \nabla^2 u(r,\theta,t) \\[5mm] u(R_0,\theta,t) = 0 $$

with $R_0 = 1$ the radius of the drumhead and $u(r,\theta,t)$ the elevation of the drumhead. The solution is given in terms of a superposition of its vibrational modes:

$$ u_{n,k}(r,\theta,t) = \cos(t) \cos(n\theta) \mathrm{J}_n(\lambda_{n,k} r) $$

where $\lambda_{n,k}$ is the $k^{th}$ positive root of $\mathrm{J}_n$ and we assumed $c=1$. Don't worry, if you don't understand the equation, in this assignment we will only plot the solutions step by step.

Create a function called drumhead_height that calculates $u_{n,k}(r,\theta,t)$. This function should take as argument $n, k,$ the radius $r$, the angle $\theta$ and the time $t$. In order to calculate the Bessel functions import:

from scipy import *

from scipy.special import jn, jn_zeros

The function jn(n,x) gives the Bessel function $J_n(x)$, and the function jn_zeros(n,k) gives you a list with the $k$ zeros of $J_n$ function.

Create a fuction called plot_Bessel that takes time $t$ as an argument. Inside this function create a figure with the GridSpec() function as in the example in the notebook; the grid should be 3x3 units. You are going to make four subplots. The first subplot is going to occupy the three rows and the two first columns. The next subplots will each occupy one row of the third column of your grid.
At the beggining of the plot_Bessel function, create two arrays with linspace. One for the radius $r$ from $0$ to $1$ and another one for the angle $\theta$ from $0$ to $2\pi$. Then use the function meshgrid, to create the polar coordintes grid ($r,\theta$). Use the output just obtained and transform it to cartesian grid using,

$$ x=r\cos(\theta)\\ y=r\sin(\theta) $$

For each of the subplots, calculate the value of the elevation of the drumhead using the function you created in 1). Remeber the drumhead_height function takes as input the grid in polar coordinates. For the first subplot use the values $n=2$ and $k=3$ to calculate your solution. For the other three subplots, you can choose your favorite values of $n$ and $k$.
In each of the subplots, use the function plot_surface from the 3D plotting library mplot3d to plot the surface (just as in the class examples). Note you will require to input the cartesian grid in the plotting function (not the polar grid).
Create 5 of these figures for different times by calling the plot_Bessel for 5 different times. Then install and use the JSAnimation package to produce an html animation that loops through the 5 different frames, make sure to embed the images into the html file and upload it as part of your homework. You can run it for 100 or more time frames in your computer, but please only submit one with 5 frames. We don't want giant files in github.

In [15]: