Scientific Visualization

The main focus here is scientific visualization.


Interactive Plotting Using Bokeh

This section is prepared for the Penn State Geoscience Computer Users Group and presented on 02/24/2016. Some examples are based on the Bokeh documentation (http://bokeh.pydata.org/).

You can download the Python Notebook that I used during the meeting here Download Notebook.

Install Bokeh (via Anaconda)

conda install bokeh

Why use interactive plotting?

  • View Data or Model at Different Scales
  • Visualize Changes in Data or Model
  • Show Links Between Data Sets and Model
  • Present High Demensional Results

Some Examples

Load Libraries

# Basic libaries
from bokeh.plotting import Figure, show, output_file
# Use to combine multiple figures
from bokeh.io import HBox, VBox
# For interactions
from bokeh.models import CustomJS, ColumnDataSource, Slider
# For data generation
import numpy as np

X-Y Plot

# some data points
x = [1, 2, 3, 4, 5]
y = [6, 7, 2, 4, 5]
# create a blank figure
p = Figure()
# plot a line
p.line(x,y, line_width=2)
# save to a html file
output_file('xyplot.html')
show(p)

Customize the Figure

x = [1, 2, 3, 4, 5]
y = [6, 7, 2, 4, 5]
# select tools
mytools = ['pan,wheel_zoom,box_zoom,reset,resize,save,crosshair']
# set up the frame
p = Figure(plot_width=800, plot_height=600, \
           tools=mytools, title='My Bokeh Plot')
# plot data
p.line(x,y, line_width=2, legend='Line')
p.circle(x,y,size=15,color='black',legend='Dots')
# add labels for axes
p.xaxis.axis_label = 'X'
p.yaxis.axis_label = 'Y'
# add legend
p.legend
# change fontsize for legend labels
p.legend.label_text_font_size = '10p'
# save to a html file
output_file('customized.html',title='Line Plot')
show(p)

Filled Contour

# generate some data points
N = 1000
x = np.linspace(0, 10, N)
y = np.linspace(0, 10, N)
xx, yy = np.meshgrid(x, y)
d = np.sin(xx)*np.cos(yy)
# set up the frame
p = Figure(x_range=(0, 10), y_range=(0, 10),\
           plot_width=800, plot_height=800)
# plot the image
p.image(image=[d], x=0, y=0, dw=10, dh=10, \
        palette='Spectral11')
# save to a html file
output_file('contour.html',title='Contour Plot')
show(p)

Click to see the result

Slider Bar

Initial plot

# create some data points
x = [x*0.005 for x in range(0, 200)]
y = x
# save the data in a bokeh format
dataSource = ColumnDataSource(data=dict(x=x, y=y))
# set up the frame
plot = Figure(plot_width=800, plot_height=600)
# plot a line, use data from 'source'
plot.line('x', 'y', source=dataSource, line_width=3, line_alpha=0.6)

Add Interactions

# define interactions
callback = CustomJS(args=dict(source=dataSource), code="""
        var data = source.get('data');
        var f = cb_obj.get('value')
        x = data['x']
        y = data['y']
        for (i = 0; i < x.length; i++) {
            y[i] = Math.pow(x[i], f)
        }
        source.trigger('change');
    """)
# set up the slider bar
slider = Slider(start=0.1, end=4, value=1, step=.1, title="power", \
                callback=callback)
layout = HBox(VBox(slider, plot))
# save to a html file
output_file('sliderBar.html')
show(layout)

Click to see the result

Linked Plots

Initial plot

# generate some data
x = np.random.uniform(0,1,500)
y = np.random.uniform(0,1,500)
# save into bokeh format
dataSet = ColumnDataSource(data=dict(x=x, y=y))
# plot left figure
leftFigure = Figure(plot_width=600, plot_height=600,
    tools="lasso_select", title="Select Here")
leftFigure.circle('x', 'y', source=dataSet, size=10, alpha=0.6)
# plot initial figure
selectedData = ColumnDataSource(data=dict(x=[], y=[]))
rightFigure = Figure(plot_width=600, plot_height=600, \
    x_range=(0, 1), y_range=(0, 1), tools="", title="Watch Here")
rightFigure.circle('x', 'y', source=selectedData, size=10, alpha=0.6

Add Interactions

# define the interaction
dataSet.callback = CustomJS(args=dict(sD=selectedData), code="""
        var inds = cb_obj.get('selected')['1d'].indices;
        var allData = cb_obj.get('data');
        var selectedData = sD.get('data');
        selectedData['x'] = []
        selectedData['y'] = []
        for (i = 0; i < inds.length; i++) {
            selectedData['x'].push(allData['x'][inds[i]])
            selectedData['y'].push(allData['y'][inds[i]])
        }
        sD.trigger('change');
    """)
# control the layout
layout = HBox(leftFigure, rightFigure, width=1600)
# save figures to a html file
output_file('linkedPlots.html')
show(layout)

Click to see the result


Plotting Data Using a Python Package - matplotlib

This section is prepared for the Penn State Geoscience Computer Users Group and presented on 11/04/2015. You may get the code and data (and examples from two other students) from this shared folder Link .

X-Y Plot

Import Packages and Set Plotting Parameters

import numpy as np
#%matplotlib inline #for python notebook
import matplotlib.pyplot as plt
from obspy import UTCDateTime
import datetime
plt.rcParams['xtick.labelsize'] = 12
plt.rcParams['ytick.labelsize'] = 12
plt.rcParams['font.sans-serif'] = 'Helvetica'
plt.rcParams['xtick.major.size'] = 8
plt.rcParams['xtick.minor.size'] = 4
plt.rcParams['ytick.major.size'] = 8
plt.rcParams['ytick.minor.size'] = 4

Load Data for X-Y Plot

data = np.loadtxt('liquefactionData.txt',skiprows=1)
x = data[:,0]
y = data[:,1]

Plot X-Y

fig = plt.figure(0,figsize=(6,4.5))
ax = plt.subplot(111)
ax.plot(x,y,'o',color='gray',markeredgecolor='black',alpha=1.0)
#ax.set_xlim(0,470)
ax.set_ylabel('Distance (km)',fontsize=12)
ax.set_xlabel('Magnitude',fontsize=12)
ax.minorticks_on()
ax.set_yscale('log')
plt.savefig('liquefaction.pdf')

Here is the result.

../_images/liquefaction.png

Date Plot

Load Data for Date Plot

eventTime, eventLat, eventLon = [], [], []
eventDepth, eventMag, eventMagType = [], [], []
eventPlace = []
fid = open('usgsLargeEQ-1900-2015.txt','r')
fid.readline()
for aline in fid:
    words = aline.split(',')
    eventTime.append(UTCDateTime(words[0]).datetime)
    eventLat.append(float(words[1]))
    eventLon.append(float(words[2]))
    eventDepth.append(float(words[3]))
    eventMag.append(float(words[4]))
    eventMagType.append(words[5])
    eventPlace.append(words[13])
fid.close()
# prepare data file for the map
fid = open('usgsLargeEQ-1900-2015-gmt.txt','w')
fid.write('      {0:12s} {1:12s} {2:12s} {3:12s} \n'.format('Lon','Lat',\
                                        'Depth','Mag'))
for i in range(len(eventLat)):
    fid.write('{0:12.4f} {1:12.4f} {2:12.4f} {3:12.4f} \n'.format(\
                eventLon[i],eventLat[i],eventDepth[i],eventMag[i]))
fid.close()

Plot Figure

from matplotlib.dates import YearLocator, MonthLocator, DateFormatter
from matplotlib.ticker import MultipleLocator
fig, ax = plt.subplots(figsize=(9,4.0))
for i in range(len(eventTime)):
    time = eventTime[i]
    mag = eventMag[i]
    size = mag**4/600
    if mag > 8:
        ax.plot_date([time,time],[5,mag],'gray',alpha=1.0,lw=1)
    ax.plot_date(time,mag,'o',color='gray',alpha=1.0,ms=size,\
                 markeredgecolor='black')
ax.set_ylim(6,10)
ax.set_xlim(datetime.date(1897,1,1),datetime.date(2020,1,1))
years = YearLocator(15)
ax.xaxis.set_major_locator(years)
years = YearLocator(1)
ax.minorticks_on()
ax.xaxis.set_minor_locator(years)
#
minorLocator = MultipleLocator(0.1)
ax.yaxis.set_minor_locator(minorLocator)
ax.set_ylabel('Magnitude',fontsize=12)
ax.set_xlabel('Year',fontsize=12)
plt.savefig('largeEQDate.pdf')

Here is the result.

../_images/largeEQDate.png

Plotting Earthquake Data on a Map Using GMT

This section is prepared for the Penn State Geoscience Computer Users Group and presented on 11/04/2015.

#!/bin/sh
gmt set FONT_ANNOT_PRIMARY 17p,Helvetica
gmt set FONT_LABEL 17p,Helvetica
gmt set PS_MEDIA 16ix8i
gmt set PROJ_LENGTH_UNIT p
out=`basename $0 .sh`.ps
#
PROJ=-JN13i
BOUNDS=-Rg0/360/-90/90
CPT=us.cpt
#
gmt psbasemap $PROJ $BOUNDS -K -B90/15 -P -Xc-0.5i -Yc > $out
# plot elevation using ETOPO1, requires topography data
#gmt grdimage ETOPO1_Bed_g_gmt4_01.grd -IETOPO1.grd -C$CPT $PROJ $BOUNDS -K -O -t50 >> $out
# plot plate boundaries, requires plate boundary data
gmt psxy $PROJ $BOUNDS PB2002_boundaries.gmt -W2p,gray -Sf2.25/12p -K -O >> $out
# plot coast lines and political boundaries
gmt pscoast $PROJ $BOUNDS -W1p -Dl -N1/2p -A50000 -K -O >> $out
#
gmt makecpt -Cseis -T-100/700/10 -D > temp.cpt
# plot earthquakes
gmt psxy $PROJ $BOUNDS pltTemp.txt -Sc -Ctemp.cpt -K -O -W0.1p -t20 >> $out
#
# plot colorbar
gmt psscale -D13.7i/3.5i/4i/0.1i -Ctemp.cpt -E -Ba100g20:"Earthquake Depth":/:km: -G0/4500 -K -O >> $out
ps2pdf $out
outpdf=`basename $0 .sh`.pdf

Here is the result.

../_images/globalEQ.png