edapy_02

In [1]:

    # edapy_02_00 clear all variables and imports various modules
    # This cell must be run first, before any of the others

    %reset -f
    import os
    from datetime import date
    from math import exp, pi, sin, sqrt, floor, ceil
    import numpy as np
    import scipy.linalg as la
    import matplotlib
    import matplotlib.pyplot as plt
    import matplotlib.cm as cm
    from matplotlib.colors import ListedColormap

    # function to make a numpy N-by-1 column vector
    # c=eda_cvec1(v) from entity v that is array-like,
    # including a number, a list of numbers, a tuple
    # of numbers, an N-by-0 np array and a N-by1 np array.
    # It promotes integers to floats, and integers
    # and floats to complex, by context. The function
    # also insures that, if v is an np array, that
    # c is a copy, as contrasted to a view, of it.
    # This version takes just one input argment.
    # whereas c=eda_cvec(v1,v2,...) concatenates
    # many argiments.
    def eda_cvec1(v):
        if isinstance(v, int):
            w = np.zeros((1,1),dtype=int);
            w[0,0] = v;
            return w;
        elif isinstance(v, float):
            w = np.zeros((1,1),dtype=float);
            w[0,0] = v;
            return w;
        elif isinstance(v, complex):
            w = np.zeros((1,1),dtype=complex);
            w[0,0] = v;
            return w;
        elif isinstance(v, np.ndarray):
            s = np.shape(v);
            if len(s) == 1:
                return np.copy(np.reshape(v,(s[0],1)));
            else:
                [r,c]=s;
                if( c==1 ):
                    return(np.copy(v));
                elif(r==1):
                    return(np.copy(v.T));
                else:
                    raise TypeError("eda_cvec: %d by %d ndarray not allowed" % (r, c));
        elif isinstance(v, list):
            r = len(v);
            t = int;
            for vi in v:
                if isinstance(vi,int):
                    pass;
                elif isinstance(vi,float):
                    t=float;
                elif isinstance(vi,complex):
                    t=complex;
                    break;
                else:
                    raise TypeError("eda_cvec: list contains unsupported type %s" % type(vi));
            w = np.zeros((r,1),dtype=t);
            w[:,0] = v;
            return w;
        elif isinstance(v, tuple):
            r = len(v);
            t = int;
            for vi in v:
                if isinstance(vi,int):
                    pass;
                elif isinstance(vi,float):
                    t=float;
                elif isinstance(vi,complex):
                    t=complex;
                    break;
                else:
                    raise TypeError("eda_cvec: tuple contains unsupported type %s" % type(vi));
            w = np.zeros((r,1),dtype=t);
            w[:,0] = v;
            return w;
        else:
            raise TypeError("eda_cvec: %s not supported" % type(v));

    # function to make a numpy N-by-1 column vector
    # c=eda_cvec(v1, v2, ...) from a list of several
    # array-like entities v1, v2, including a number
    # a list of numbers, a tuple of numbers, an N-by-0 np array
    # and a N-by-1 np array. The function
    # also insures that, if v is an np array, that
    # c is a copy, as contrasted to a view, of it
    # It promotes integers to floats, and integers
    # and floats to complex, by context.
    # This version concatenates many argments,
    # whereas c=eda_cvec1(v1) takes just one argument.
    # I recommend always using eda_cvec(v1, v2, ...)
    def eda_cvec(*argv):
        t = int;
        Nt = 0;
        for a in argv:
            v = eda_cvec1(a);
            N,M = np.shape(v);
            Nt = Nt + N;
            if( N==0 ):
                continue; # skip vector of zero length
            if (t==int) and isinstance(v[0,0],float):
                t=float;
            elif isinstance(v[0,0],complex):
                t=complex;
        w = np.zeros((Nt,1),dtype=t);
        Nt = 0;
        for a in argv:
            v = eda_cvec1(a);
            N,M = np.shape(v);
            w[Nt:Nt+N,0] = v[0:N,0];
            Nt = Nt + N;
        return w;

In [2]:

    # edapy_02_01, print numeric values of a few data
    # Black Rock Forest temperature data used as an example

    # load data from text file into matrix D
    D = np.genfromtxt('../Data/brf_temp.txt', delimiter='\t');
    t = eda_cvec( D[:,0] ); # time in col 0
    d = eda_cvec( D[:,1] ); # temperature data in col 1
    [N, M]=D.shape;  # D has rows N and columns

    print('rows ',N,'cols ',M);

    # print first few rows
    print( D[1:5,:] );

    rows  110430 cols  2
    [[  0.042 -17.85 ]
     [  0.083 -18.42 ]
     [  0.125 -18.94 ]
     [  0.167 -19.29 ]]

In [31]:

    # edapy_02_02: initial plot of the data
    # Black Rock Forest temperature data used as an example

    # load data from text file into matrix D
    D = np.genfromtxt('../Data/brf_temp.txt', delimiter='\t');
    t = eda_cvec( D[:,0] ); # time in col 0
    d = eda_cvec( D[:,1] ); # temperature data in col 1
    [N, M]=D.shape;  # determine rows N and columns M of D

    # print size and frst few rows
    print('rows ',N,'cols ',M);

    # print first few rows
    print( D[1:5,:] );

    # load data from text file
    D = np.genfromtxt('../Data/brf_temp.txt', delimiter='\t');
    [N, M]=D.shape;
    t = eda_cvec( D[:,0] );
    d = eda_cvec( D[:,1] );

    #  plot the data
    fig1 = plt.figure(1,figsize=(8,4)); # open figure 1
    ax1 = plt.subplot(1,1,1); # single plot figure
    plt.axis([0, 5000, -60, 40])  # plot bounds
    plt.plot(t,d,'k-');           # plot data
    plt.xlabel('Time (days)');    # horizontal axis label
    plt.ylabel('Teperature (C)'); # vertical axis label
    plt.title('Black Rock Forest weather station'); # title
    plt.show(); # display finished figure

In [34]:

    # edapy_02_03, enlarged plot of part of the data
    # Black Rock Forest temperature data used as an example

    # load data from text file
    D = np.genfromtxt('../Data/brf_temp.txt', delimiter='\t')
    [N, M]=D.shape;
    t = eda_cvec( D[:,0] );
    d = eda_cvec( D[:,1] );

    tc = 1560.0;        # center tine in seconds
    Dt = t[2,0]-t[1,0]; # sampling interval
    wh = 25.0;          # half width of window in days
    iwh = floor(wh/Dt); # half width of window in samples
    ii = np.where(t>tc); # find first index of t greater than tc
    i = ii[0][0];

    # close-up plot, centered at tc, 2*w2 wide
    fig1 = plt.figure(1,figsize=(8,4));
    ax3 = plt.subplot(1,1,1);
    # note: plt.axis() omitted, so plot auto-scales
    plt.plot(t[i-iwh:i+iwh,0],d[i-iwh:i+iwh,0],'k-');
    plt.plot(t[i-iwh:i+iwh,0],d[i-iwh:i+iwh,0],'k*');
    plt.xlabel('time after Jan 1, 1997 (days)');
    plt.ylabel('teperature (C)');
    plt.title('Close-up: Black Rock Forest Temp');
    plt.show();

    # whole data figure with all data in grey, selected data in black
    fig2 = plt.figure(2,figsize=(8,4));
    ax2 = plt.subplot(1,1,1);
    plt.axis([0, 5000, -60, 40])
    plt.plot(t,d,'k-',linewidth=2,color=(0.8,0.8,0.8));
    plt.plot(t[i-iwh:i+iwh,0],d[i-iwh:i+iwh,0],'k-');
    plt.xlabel('Time (days)');
    plt.ylabel('Teperature (C)');
    plt.title('Black Rock Forest weather station');
    plt.show;

In [37]:

    # edapy_02_04: enlarged plot via mouse-click
    # Black Rock Forest temperature data used as an example

    #  interactive graphics from withing the Jupyter Notebook is a bit dicey
    # if METHOD=1 fails, try METHOD=2 (or vice-versa)
    METHOD=1;

    # get current backend
    if( METHOD==1):
        bkend = matplotlib.get_backend();

    # load data from text file
    D = np.genfromtxt('../Data/brf_temp.txt', delimiter='\t')
    [N, M]=D.shape;
    t = eda_cvec( D[:,0] );
    d = eda_cvec( D[:,1] );

    #  use interactive plotting
    if( METHOD==1 ):
        matplotlib.use('Qt5Agg');
    else:
        %matplotlib qt5
    fig1 = plt.figure(1,figsize=(8,4)); # interactive figure
    ax1 = plt.subplot(1,1,1);
    plt.axis([0, 5000, -60, 40])
    plt.plot(t,d,'k-');
    plt.xlabel('Time (days)');
    plt.ylabel('Teperature (C)');
    plt.title('Black Rock Forest weather station');
    plt.show();
    Dt = t[2,0]-t[1,0]; # sampling interval
    wh = 25.0;          # half width of window in days
    iwh = floor(wh/Dt); # half width of window in samples
    v=plt.ginput(1); # get (t,d) of 1 mouse click
    plt.close();     # close interactive figure
    tc=v[0][0];      # time of click
    dc=v[0][1];      # temperature of click
    ii = np.where(t>tc); # find first index of t greater than tc
    i = ii[0][0];
    print('click at %.2f with data d[%d]=%.2f' % (tc, i, d[i,0]) );

    # switch back to standard plotting
    if( METHOD==1 ):
        matplotlib.use(bkend);
    else:
        %matplotlib inline

    fig2 = plt.figure(2,figsize=(8,4)); # close-up figure
    ax3 = plt.subplot(1,1,1);
    plt.plot(t[i-iwh:i+iwh,0],d[i-iwh:i+iwh,0],'k-');
    plt.plot(t[i-iwh:i+iwh,0],d[i-iwh:i+iwh,0],'k*');
    plt.xlabel('time after Jan 1, 1997 (days)');
    plt.ylabel('teperature (C)');
    plt.title('Enlarged Section');
    plt.show();

    # whole data figure with all data in grey, selected data in black
    fig3 = plt.figure(3,figsize=(8,4));
    ax2 = plt.subplot(1,1,1);
    plt.axis([0, 5000, -60, 40])
    plt.plot(t,d,'k-',linewidth=2,color=(0.8,0.8,0.8));
    plt.plot(t[i-iwh:i+iwh,0],d[i-iwh:i+iwh,0],'k-');
    plt.xlabel('Time (days)');
    plt.ylabel('Teperature (C)');
    plt.title('Black Rock Forest weather station');
    plt.show;

    click at 1568.10 with data d[32638]=2.79

In [42]:

    # edapy_02_05: histogram of data
    # Black Rock Forest temperature data used as an example

    # load data from text file
    D = np.genfromtxt('../Data/brf_temp.txt', delimiter='\t')
    [N, M]=D.shape;
    t = eda_cvec( D[:,0] );
    d = eda_cvec( D[:,1] );

    # histogram
    Lh = 100;            # number of bins in histogram
    dmin = np.min(d);    # minimum bin
    dmax = np.max(d);    # maximum bin
    c, e = np.histogram(d,Lh,(dmin,dmax));  # create histogram
    Nc = len(c);  # lengths of counts
    Ne = len(e);  # length of edges
    counts = eda_cvec(c);            # vector of counts
    edges = eda_cvec( e[0:Ne-1] );  # vector of edges
    centers = edges + 0.5*(edges[1,0]-edges[0,0]); # centers

    # plot histogram
    fig1 = plt.figure(1,figsize=(8,4));
    ax1 = plt.subplot(1,1,1);

    plt.axis([-60, 40, 0, np.max(c)]);
    plt.plot(centers,counts,'k-');
    plt.xlabel('temperature (C)');
    plt.ylabel('counts');
    plt.title('Histogram of Temperature Measurements');
    plt.show();

In [57]:

    # edapy_02_06: plot histogram two ways
    # Black Rock Forest temperature data used as an example

    # load data from text file
    D = np.genfromtxt('../Data/brf_temp.txt', delimiter='\t')
    [N, M]=D.shape;
    t = eda_cvec( np.copy( D[:,0] ) );
    d = eda_cvec( np.copy( D[:,1] ) );

    # histogram
    Lh = 100;         # number of bins
    dmin = np.min(d); # minimum bin
    dmax = np.max(d); # maximum bin
    c, e = np.histogram(d,Lh,(dmin,dmax)); # compute histogram
    Nc = len(c); # length of couts
    Ne = len(e); # length of edges
    counts = eda_cvec(c);          # vector of counts
    edges = eda_cvec(e[0:Ne-1]);  # vector of edges
    centers = edges + 0.5*(edges[1,0]-edges[0,0]); # centers

    # plot
    fig1 = plt.figure(1,figsize=(8,8));

    # subplot 1; traditional line plot
    ax1 = plt.subplot(1, 2, 1);
    plt.axis([0, np.max(counts), -60, 40]);
    plt.gca().invert_yaxis();  # y-axis positive down
    plt.plot(counts,centers,'k-');
    plt.xlabel('counts');
    plt.ylabel('temperature (C)');

    # subplot 1; grey-scale image
    ax2 = plt.subplot(1, 2, 2);
    mycmap = cm.gray;             # grey-scale colormap
    mycmap = mycmap.reversed();   # reverse intensities, so black=max
    plt.axis([0.0, 1.0, -60, 40]);
    myvmax = np.max(count);
    myextent=(0.0,1.0,np.max(centers),np.min(centers));
    plt.gca().invert_yaxis();  # y-axis positive down
    plt.imshow( count, cmap=mycmap, vmin=0, vmax=myvmax, extent=myextent, aspect='auto' );
    plt.colorbar(label='counts'); # draw color bar
    plt.show();

In [4]:

    # edapy_02_07, moving-window histogram of data

    # load temperature data from text fole
    D = np.genfromtxt('../Data/brf_temp.txt', delimiter='\t')
    [N, M]=D.shape;
    t = eda_cvec(D[:,0]);
    d = eda_cvec(D[:,1]);
    Dt = 1.0/24.0;

    # There are Lw histograms, each offset apart
    toffset = 30.0;   # offset of 30 days
    offset = floor(toffset/Dt); # offset in samples
    Lw = floor( N/offset ); # number of windows
    Lh = 100;  # number of bins in each histogram
    tend = toffset*Lw; # last time

    # an array to hold the histograms
    Dw = np.zeros((Lh,Lw));

    # histograms
    dmin = -60;
    dmax = 40;
    for i in range(Lw): # loop over time windows
        j=(i-1)*offset; # index of first data in window
        k=j+offset;     # index of last data in window
        c, e = np.histogram(d[j:k,0],Lh,(dmin,dmax));
        Nc = len(c);
        counts = eda_cvec(c);
        Dw[0:Lh,i:i+1]=counts;

    # plot
    fig1 = plt.figure(1,figsize=(8,4));
    ax1 = plt.subplot(1,1,1);
    plt.rcParams['xtick.bottom'] = False;
    plt.rcParams['xtick.labelbottom'] = False;
    plt.rcParams['xtick.top'] = True;
    plt.rcParams['xtick.labeltop'] = True;
    plt.xlabel('time (days)');
    plt.ylabel('temperature (C)');
    ax1.xaxis.set_label_position('top');
    plt.axis([0, tend, dmin, dmax]);
    plt.gca().invert_yaxis();

    # make a greyscale colormap that is never completely white
    bw = np.zeros((256,4));
    v = eda_cvec( 0.9*(256 - np.linspace( 0, 255, 256 ))/255 );
    w = np.ones((256,1));
    bw[0:256,0:1] = v;
    bw[0:256,1:2] = v;
    bw[0:256,2:3] = v;
    bw[0:256,3:4] = w;
    bwcmap = ListedColormap(bw);

    # alternatively, could have ued colormap
    # bwcmap = cm.gray;
    # bwcmap = bwcmap.reversed();

    plt.imshow( Dw, cmap=bwcmap, vmin=0, vmax=np.max(Dw), aspect='auto', extent=(0,tend,dmin,dmax), origin='lower' );
    plt.colorbar(label='counts');
    plt.show();

In [74]:

    # edapy_02_08: superimposed plots with different line types

    # time axis with N points, starting at tmin, with increment Dt
    Dt=1;
    N=511;
    tmin=0;
    tmax=tmin+Dt*(N-1);
    t = eda_cvec( np.linspace( tmin, tmax, N ) );

    # trig functions for sample data
    d1=np.sin(pi*t/tmax);
    d2=np.sin(2*pi*t/tmax);

    # plot the sample data
    fig1 = plt.figure(1,figsize=(8,8));
    ax1 = plt.subplot(1,1,1);
    plt.axis([tmin, tmax, -1.1, 1.1]);
    plt.plot(t,d1,'k-');
    plt.plot(t,d2,'k:');
    plt.xlabel('time (s)');
    plt.ylabel('data');
    plt.title('data consisting of sine waves');

In [78]:

    # edapy_02_09: sideways side-by-side plots

    # time axis with N points, starting at tmin, with increment Dt
    Dt=1;
    N=511;
    tmin=0;
    tmax=tmin+Dt*(N-1);
    t = eda_cvec( np.linspace( tmin, tmax, N ) );

    #  trig functions for sample data
    d1=np.sin(pi*t/tmax);
    d2=np.sin(2*pi*t/tmax);

    # plot the sample data
    fig1 = plt.figure(1,figsize=(10,8));
    # change the position of the tics
    plt.rcParams['xtick.bottom'] = False;
    plt.rcParams['xtick.labelbottom'] = False;
    plt.rcParams['xtick.top'] = True;
    plt.rcParams['xtick.labeltop'] = True;

    # subplot 1 of d1
    ax1 = plt.subplot(1,2,1);
    ax1.xaxis.set_label_position('top');
    plt.axis([-1.1, 1.1, tmin, tmax, ]);
    ax1.invert_yaxis();
    plt.plot(d1, t,'k-'); # plot data d1
    plt.plot([0,0],[tmin,tmax],'k:'); # plot zero line
    plt.xlabel('d1');
    plt.ylabel('time (s)');

    # subplot 2 of d2
    ax2 = plt.subplot(1,2,2);
    ax2.xaxis.set_label_position('top');
    plt.axis([-1.1, 1.1, tmin, tmax]);
    ax2.invert_yaxis();
    plt.plot(d2, t,'k-'); # plot data d2
    plt.plot([0,0],[tmin,tmax],'k:');  # plot zero line
    plt.xlabel('d2');
    plt.ylabel('time (s)');

In [5]:

    # edapy_02_10: colormaps
    # note you must: import matplotlib.cm as cm

    # test image
    N = 11;
    D = np.random.uniform( low=0.0, high=1.0, size=(N,N) );

    # plot with Python's gray colormap
    fig1 = plt.figure(1,figsize=(8,4));
    ax1 = plt.subplot(1,1,1);
    # Python-provided grey-scale colormap
    graycm = cm.get_cmap('gray', 256);
    plt.imshow( D, cmap=graycm, vmin=0.0, vmax=1.0 );
    plt.title('image plotted with Python''s grey colormap');
    plt.show();

    # alternaitive is to build a user-defined colormap
    n=256;
    bw = np.zeros((n,4));
    # note that this colormap is reversed; white is low
    v = eda_cvec( 0.9*(n - np.linspace( 0.0, n-1, n ))/((n-1)) );
    w = np.ones((n,1));
    bw[0:n,0:1] = v;
    bw[0:n,1:2] = v;
    bw[0:n,2:3] = v;
    bw[0:n,3:4] = w;
    usercmap = ListedColormap(bw);

    # plot
    fig2 = plt.figure(2,figsize=(8,4));
    ax1 = plt.subplot(1,1,1);
    plt.imshow( D, cmap=usercmap, vmin=0.0, vmax=1.0 );
    plt.title('image plotted with user-defined colormap');
    plt.show();

In [137]:

    # edapy_02_11, example of plotting an image
    # based on edapy_02_08: running histogram of temperature

    # load data from text fild
    D = np.genfromtxt('../Data/brf_temp.txt', delimiter='\t')
    [N, M]=D.shape;
    t = eda_cvec(D[:,0]);
    d = eda_cvec(D[:,1]);

    # moving window histogram
    offset = 1000;
    Lh = 100;
    Lw = floor( N/offset );
    Dw = np.zeros((Lh,Lw));
    tend = offset*(Lw-1); # last time

    # make image of histograms
    dmin = -60;
    dmax = 40;
    for i in range(Lw):
        j=(i-1)*offset;
        k=j+offset;
        c, e = np.histogram(d[j:k,0],Lh,(dmin,dmax));
        Nc = len(c);
        counts = eda_cvec(c);
        Dw[:,i]=counts.ravel();

    # greyscale colormap
    bw = np.zeros((256,4));
    v = eda_cvec( 0.9*(256 - np.linspace( 0, 255, 256 ))/255 );
    w = np.ones((256,1));
    bw[0:n,0:1] = v;
    bw[0:n,1:2] = v;
    bw[0:n,2:3] = v;
    bw[0:n,3:4] = w;
    bwcmap = ListedColormap(bw);

    # plot image; note image in matrix Dw
    fig1 = plt.figure(1,figsize=(8,4)); # open figure
    ax1 = plt.subplot(1,1,1);       # one plot in figure
    plt.rcParams['xtick.bottom'] = False;  # no ticks on bottom
    plt.rcParams['xtick.labelbottom'] = False; # no tick labels
                                               # on bottom
    plt.rcParams['xtick.top'] = True; # ticks on top
    plt.rcParams['xtick.labeltop'] = True; # tick labels on top
    ax1.xaxis.set_label_position('top'); # axis label on top
    vlo = 0.0;          # low value of color map
    vhi = np.max(Dw);   # high value of color map
    ex=(0,tend,dmin,dmax); # corners of image plot here
    plt.imshow(Dw, cmap=bwcmap, vmin=vlo, vmax=vhi,
      aspect='auto', extent=ex, origin='lower' ); # plot matrix Dw
    ax1.invert_yaxis();  # vertical axis increases downward
    plt.colorbar(label='counts');  # plot colrbar with label
    plt.xlabel('time (days)');  # label horizontal axis
    plt.ylabel('temperature (C)'); # label verical axis
    plt.title('moving window histogram'); # title
    plt.show(); # plot finished figure

In [149]:

    # edapy_02_12: plot and histgram of rates of change
    # using the Neuse River Hydrograph data as an example

    # load data from text fole
    D = np.genfromtxt('../Data/neuse.txt', delimiter='\t')
    [N, M]=D.shape;

    # time in column 1, discharge in colun 2
    t = eda_cvec(D[:,0]);
    d = eda_cvec(D[:,1]);

    # min, max times for plotting
    tmin=500;
    tmax = 1000;

    # compute rate
    dddt = np.divide( d[1:N,0:1]-d[0:N-1,0:1],
                            t[1:N,0:1]-t[0:N-1,0:1] );

    # plot the sample data
    fig1 = plt.figure(1,figsize=(16,8));
    plt.rcParams['xtick.bottom'] = False;
    plt.rcParams['xtick.labelbottom'] = False;
    plt.rcParams['xtick.top'] = True;
    plt.rcParams['xtick.labeltop'] = True;

    # subplot 1 of data
    ax1 = plt.subplot(1,3,1);
    ax1.xaxis.set_label_position('top');
    plt.axis([ 0, 20000, tmin, tmax,]);
    plt.plot(d,t,'k-');
    ax1.invert_yaxis();
    plt.xlabel('discharge d');
    plt.ylabel('time (days)');

    # subplot 2 of derivative
    ax2 = plt.subplot(1,3,2);
    ax2.xaxis.set_label_position('top');
    plt.axis([-5000, 5000, tmin, tmax, ]);
    ax2.invert_yaxis();
    plt.plot(dddt,t[0:N-1],'k-');
    plt.xlabel('derivative dd/dt');
    plt.ylabel('time (days)');

    # histogram
    Lh = 1000;
    dmin = np.min(dddt);
    dmax = np.max(dddt);
    c, e = np.histogram(dddt,Lh,(dmin,dmax));
    Dbin = (dmin-dmax)/(Lh-1);
    Nc = len(c);
    Ne = len(e);
    count = eda_cvec(c);
    edges = eda_cvec(e[0:Ne-1]);
    centers = edges + Dbin/2.0;

    # subplot 3 of histogram
    ax3 = plt.subplot(1, 3, 3);
    ax3.xaxis.set_label_position('top');
    plt.axis([0, np.max(c), -500, 500]);
    ax3.invert_yaxis();
    plt.plot(count,centers,'k-');
    plt.plot([0,500],[0,0],'k:');
    plt.xlabel('counts');
    plt.ylabel('dddt');

In [156]:

    # edapy_02_13: segregating data based on the sign of rate
    # using the Neuse River Hydrograph data as an example

    # load data from text fole
    D = np.genfromtxt('../Data/neuse.txt', delimiter='\t')
    [N, M]=D.shape;

    # time in column 1, discharge in colun 2
    t = D[0:N,0:1];
    d = D[0:N,1:2];

    # derivative
    dddt = np.divide( d[1:N,0:1]-d[0:N-1,0:1],
                            t[1:N,0:1]-t[0:N-1,0:1] );

    # divide (d, dddt) pairs into two groups depending
    # upon whether dddt is positive or negative
    pos = np.where( dddt>=0 );
    neg = np.where( dddt<0 );
    dpos=eda_cvec(d[pos]);
    dddtpos=eda_cvec(dddt[pos]);
    dneg=eda_cvec(d[neg]);
    dddtneg=eda_cvec(dddt[neg]);

    # plot
    fig1 = plt.figure(1,figsize=(14,12));
    plt.rcParams['xtick.bottom'] = True;
    plt.rcParams['xtick.labelbottom'] = True;
    plt.rcParams['xtick.top'] = False;
    plt.rcParams['xtick.labeltop'] = False;

    # d vs dddt scatter plot for positive rates
    ax1 = plt.subplot(2,1,1);
    plt.axis([ 0, 10000, 0, 5000]);
    plt.plot(dpos,dddtpos,'k.');
    plt.xlabel('dischagre (cfs)');
    plt.ylabel('d/dt of discharge');
    plt.title('positive rates');

    # d vs dddt scatter plot for negative rates
    ax2 = plt.subplot(2,1,2);
    plt.axis([ 0, 10000, -5000, 0]);
    plt.plot(dneg,dddtneg,'k.');
    plt.xlabel('dischagre (cfs)');
    plt.ylabel('d/dt of discharge');
    plt.title('negative rates');

    (1861,)

In [166]:

    # edapy_02_14: scatter plots of Atlantic Rock chemical data

    # load data from text fole
    D = np.genfromtxt('../Data/rocks.txt', delimiter='\t')
    [N, M]=D.shape;

    # note the use of a list to hold names
    names = [ 'SiO2', 'TiO2', 'Al203', 'FeO-total',
               'MgO', 'CaO', 'Na2O', 'K2O'];
    Ns = len(names);

    # Np: number of permulations of Ns things taken two at a time
    Np = Ns*Ns-(Ns*(Ns+1)/2);
    M = 4; # plots in a row
    Nr = ceil(Np/M); # number of rows
    i=2;

    fig1 = plt.figure(1,figsize=(4*M,4*Nr));
    sp = 0;  # counts plots
    for i in range(0,Ns-1):
        for j in range(i+1,Ns):
               sp = sp+1; # increment counter
               ax1 = plt.subplot(Nr, M, sp);
               plt.plot(D[:,i],D[:,j],'k.');
               plt.xlabel(names[i]);
               plt.ylabel(names[j]);

In [176]:

    # edapy_02_15, character strings

    # Part 1.  String variables
    my_name = "Bill";
    my_path = "C:/bill";
    print(my_name);
    print(my_path);
    print(" ");

    # Part 2.  Using lists to hold strings
    names = [ 'lithium', 'calcium', 'sodium', 'potassium' ];
    Ns = len(names);

    print("number of names:", Ns)
    print("name 0:", names[0] );
    print("name 1:", names[1] );
    print("name 2:", names[2] );
    print("name 3:", names[3] );
    print(" ");

    # Part 3. Formatting strings using the % operator

    # prints "element 2"
    i=2;
    mystring1 = "element %d" % i;
    print(mystring1);

    # prints "a=4.72"
    a=4.7213;
    mystring2 = "a=%.2f" % a;
    print(mystring2);

    # prints "row 2 column 4"
    i=2;
    j=4;
    mystring3 = "row %d column %d" % (i, j);
    print(mystring3);

    # prints "(2,4): 4.721"
    i=2;
    j=4;
    a=4.7213;
    mystring4 = "(%d,%d): %.3f" % (i, j, a);
    print(mystring4);
    print(" ");

    # Part 4. Combining lists and formats
    for i in range(Ns):
        print("name %d is %s" % (i, names[i]));

    Bill
    C:/bill

    number of names: 4
    name 0: lithium
    name 1: calcium
    name 2: sodium
    name 3: potassium

    element 2
    a=4.72
    row 2 column 4
    (2,4): 4.721

    name 0 is lithium
    name 1 is calcium
    name 2 is sodium
    name 3 is potassium

In [ ]: