.. _03notes-01visualization:

================================
Scientific Visualization
================================

The main focus here is scientific visualization. 

-----

.. index:: Bokeh

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:`Download Notebook <./images/InteractivePlottingExamples.ipynb.zip>`. 

Install Bokeh (via Anaconda)
------------------------------

.. code-block:: bash

    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
'''''''''''''''''

.. code-block:: python 

    # 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
''''''''''''''''''

.. code-block:: python

    # 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)

.. :download:`Click to see the result <./images/xyplot.html>`

.. raw:: html

    <iframe src="../_static/xyplot.html" height="700px" width="100%"></iframe>

Customize the Figure
'''''''''''''''''''''''

.. code-block:: python

    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)

.. :download:`Click to see the result <./images/customized.html>`

.. raw:: html

    <iframe src="../_static/customized.html" height="700px" width="130%"></iframe>

Filled Contour
''''''''''''''''''

.. code-block:: python

    # 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)

:download:`Click to see the result <./images/contour.html>`

.. .. raw:: html
        <iframe src="../_static/contour.html" height="850px" width="150%"></iframe>

Slider Bar
''''''''''''''''''

Initial plot

.. code-block:: python

    # 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

.. code-block:: python

    # 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)

:download:`Click to see the result <./images/sliderBar.html>`

.. .. raw:: html
        <iframe src="../_static/sliderBar.html" height="850px" width="150%"></iframe>

Linked Plots
''''''''''''''''''

Initial plot

.. code-block:: python

    # 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

.. code-block:: python

    # 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)

:download:`Click to see the result <./images/linkedPlots.html>`

..  .. raw:: html
    
    <iframe src="../_static/linkedPlots.html" height="700px" width="200%"></iframe>

-----

.. index:: matplotlib

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 <https://psu.box.com/s/vfoddxh72ihpmb30e6vj6kry5z3ofegg>`__ .

X-Y Plot
------------------------------

Import Packages and Set Plotting Parameters
''''''''''''''''''''''''''''''''''''''''''''''''''''''

.. code-block:: python

    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
''''''''''''''''''''''''''''''''''''

.. code-block:: python

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

Plot X-Y
''''''''''''''''''''''''''''''''''''

.. code-block:: python

    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.

.. image:: ./images/liquefaction.png

Date Plot
------------------------------

Load Data for Date Plot
''''''''''''''''''''''''''''''''''''

.. code-block:: python

    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
''''''''''''''''''''''''''''''''''''

.. code-block:: python

    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.

.. image:: ./images/largeEQDate.png

-----

.. index:: GMT

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.

.. code-block:: bash

    #!/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.

.. image:: ./images/globalEQ.png