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.
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
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:
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()
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:
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()
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]:
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)
There are endless ways of producing animations; we will explore two of them.
For the first case, let make python output a set of images:
In [11]:
!mkdir animation_frames
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,
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:
where -r 10 indicates the frames per second and 1000k the bitrate (quality) of the video.
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.
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))
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.
1) 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.
2) 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.
3) 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) $$4) 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$.
5) 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).
6) 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]: