{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# edapy_02_00 clear all variables and imports various modules\n",
    "# This cell must be run first, before any of the others\n",
    "\n",
    "%reset -f\n",
    "import os\n",
    "from datetime import date\n",
    "from math import exp, pi, sin, sqrt, floor, ceil\n",
    "import numpy as np\n",
    "import scipy.linalg as la\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.cm as cm\n",
    "from matplotlib.colors import ListedColormap\n",
    "\n",
    "# function to make a numpy N-by-1 column vector\n",
    "# c=eda_cvec1(v) from entity v that is array-like,\n",
    "# including a number, a list of numbers, a tuple\n",
    "# of numbers, an N-by-0 np array and a N-by1 np array.\n",
    "# It promotes integers to floats, and integers\n",
    "# and floats to complex, by context. The function\n",
    "# also insures that, if v is an np array, that\n",
    "# c is a copy, as contrasted to a view, of it.\n",
    "# This version takes just one input argment.\n",
    "# whereas c=eda_cvec(v1,v2,...) concatenates\n",
    "# many argiments.\n",
    "def eda_cvec1(v):\n",
    "    if isinstance(v, int):\n",
    "        w = np.zeros((1,1),dtype=int);\n",
    "        w[0,0] = v;\n",
    "        return w;\n",
    "    elif isinstance(v, float):\n",
    "        w = np.zeros((1,1),dtype=float);\n",
    "        w[0,0] = v;\n",
    "        return w;\n",
    "    elif isinstance(v, complex):\n",
    "        w = np.zeros((1,1),dtype=complex);\n",
    "        w[0,0] = v;\n",
    "        return w;\n",
    "    elif isinstance(v, np.ndarray):\n",
    "        s = np.shape(v);\n",
    "        if len(s) == 1:\n",
    "            return np.copy(np.reshape(v,(s[0],1)));\n",
    "        else:\n",
    "            [r,c]=s;\n",
    "            if( c==1 ):\n",
    "                return(np.copy(v));\n",
    "            elif(r==1):\n",
    "                return(np.copy(v.T));\n",
    "            else:\n",
    "                raise TypeError(\"eda_cvec: %d by %d ndarray not allowed\" % (r, c));\n",
    "    elif isinstance(v, list):\n",
    "        r = len(v);\n",
    "        t = int;\n",
    "        for vi in v:\n",
    "            if isinstance(vi,int):\n",
    "                pass;\n",
    "            elif isinstance(vi,float):\n",
    "                t=float;\n",
    "            elif isinstance(vi,complex):\n",
    "                t=complex;\n",
    "                break;\n",
    "            else:\n",
    "                raise TypeError(\"eda_cvec: list contains unsupported type %s\" % type(vi));\n",
    "        w = np.zeros((r,1),dtype=t);\n",
    "        w[:,0] = v;\n",
    "        return w;\n",
    "    elif isinstance(v, tuple):\n",
    "        r = len(v);\n",
    "        t = int;\n",
    "        for vi in v:\n",
    "            if isinstance(vi,int):\n",
    "                pass;\n",
    "            elif isinstance(vi,float):\n",
    "                t=float;\n",
    "            elif isinstance(vi,complex):\n",
    "                t=complex;\n",
    "                break;\n",
    "            else:\n",
    "                raise TypeError(\"eda_cvec: tuple contains unsupported type %s\" % type(vi));\n",
    "        w = np.zeros((r,1),dtype=t);\n",
    "        w[:,0] = v;\n",
    "        return w;\n",
    "    else:\n",
    "        raise TypeError(\"eda_cvec: %s not supported\" % type(v));\n",
    "        \n",
    "# function to make a numpy N-by-1 column vector\n",
    "# c=eda_cvec(v1, v2, ...) from a list of several\n",
    "# array-like entities v1, v2, including a number\n",
    "# a list of numbers, a tuple of numbers, an N-by-0 np array\n",
    "# and a N-by-1 np array. The function\n",
    "# also insures that, if v is an np array, that\n",
    "# c is a copy, as contrasted to a view, of it\n",
    "# It promotes integers to floats, and integers\n",
    "# and floats to complex, by context.\n",
    "# This version concatenates many argments,\n",
    "# whereas c=eda_cvec1(v1) takes just one argument.\n",
    "# I recommend always using eda_cvec(v1, v2, ...)\n",
    "def eda_cvec(*argv):\n",
    "    t = int;\n",
    "    Nt = 0;\n",
    "    for a in argv:\n",
    "        v = eda_cvec1(a);\n",
    "        N,M = np.shape(v);\n",
    "        Nt = Nt + N;\n",
    "        if( N==0 ):\n",
    "            continue; # skip vector of zero length\n",
    "        if (t==int) and isinstance(v[0,0],float):\n",
    "            t=float;\n",
    "        elif isinstance(v[0,0],complex):\n",
    "            t=complex;\n",
    "    w = np.zeros((Nt,1),dtype=t);\n",
    "    Nt = 0;\n",
    "    for a in argv:\n",
    "        v = eda_cvec1(a);\n",
    "        N,M = np.shape(v);\n",
    "        w[Nt:Nt+N,0] = v[0:N,0];\n",
    "        Nt = Nt + N;\n",
    "    return w;\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "rows  110430 cols  2\n",
      "[[  0.042 -17.85 ]\n",
      " [  0.083 -18.42 ]\n",
      " [  0.125 -18.94 ]\n",
      " [  0.167 -19.29 ]]\n"
     ]
    }
   ],
   "source": [
    "# edapy_02_01, print numeric values of a few data\n",
    "# Black Rock Forest temperature data used as an example\n",
    "\n",
    "# load data from text file into matrix D\n",
    "D = np.genfromtxt('../Data/brf_temp.txt', delimiter='\\t');\n",
    "t = eda_cvec( D[:,0] ); # time in col 0\n",
    "d = eda_cvec( D[:,1] ); # temperature data in col 1\n",
    "[N, M]=D.shape;  # D has rows N and columns\n",
    "\n",
    "print('rows ',N,'cols ',M);\n",
    "\n",
    "# print first few rows\n",
    "print( D[1:5,:] );"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 27924 characters]\n",
      "text/plain": [
       "<Figure size 576x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# edapy_02_02: initial plot of the data\n",
    "# Black Rock Forest temperature data used as an example\n",
    "\n",
    "# load data from text file into matrix D\n",
    "D = np.genfromtxt('../Data/brf_temp.txt', delimiter='\\t');\n",
    "t = eda_cvec( D[:,0] ); # time in col 0\n",
    "d = eda_cvec( D[:,1] ); # temperature data in col 1\n",
    "[N, M]=D.shape;  # determine rows N and columns M of D\n",
    "\n",
    "# print size and frst few rows\n",
    "print('rows ',N,'cols ',M);\n",
    "\n",
    "# print first few rows\n",
    "print( D[1:5,:] );\n",
    "\n",
    "# load data from text file\n",
    "D = np.genfromtxt('../Data/brf_temp.txt', delimiter='\\t');\n",
    "[N, M]=D.shape;\n",
    "t = eda_cvec( D[:,0] );\n",
    "d = eda_cvec( D[:,1] );\n",
    "\n",
    "#  plot the data\n",
    "fig1 = plt.figure(1,figsize=(8,4)); # open figure 1\n",
    "ax1 = plt.subplot(1,1,1); # single plot figure\n",
    "plt.axis([0, 5000, -60, 40])  # plot bounds\n",
    "plt.plot(t,d,'k-');           # plot data\n",
    "plt.xlabel('Time (days)');    # horizontal axis label\n",
    "plt.ylabel('Teperature (C)'); # vertical axis label\n",
    "plt.title('Black Rock Forest weather station'); # title\n",
    "plt.show(); # display finished figure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 35616 characters]\n",
      "text/plain": [
       "<Figure size 576x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 23780 characters]\n",
      "text/plain": [
       "<Figure size 576x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# edapy_02_03, enlarged plot of part of the data\n",
    "# Black Rock Forest temperature data used as an example\n",
    "\n",
    "# load data from text file\n",
    "D = np.genfromtxt('../Data/brf_temp.txt', delimiter='\\t')\n",
    "[N, M]=D.shape;\n",
    "t = eda_cvec( D[:,0] );\n",
    "d = eda_cvec( D[:,1] );\n",
    "\n",
    "tc = 1560.0;        # center tine in seconds\n",
    "Dt = t[2,0]-t[1,0]; # sampling interval\n",
    "wh = 25.0;          # half width of window in days\n",
    "iwh = floor(wh/Dt); # half width of window in samples   \n",
    "r = np.where(t>tc); # find first index of t greater than tc\n",
    "k = r[0];\n",
    "i = k[0];\n",
    "\n",
    "# close-up plot, centered at tc, 2*w2 wide\n",
    "fig1 = plt.figure(1,figsize=(8,4));\n",
    "ax3 = plt.subplot(1,1,1);\n",
    "# note: plt.axis() omitted, so plot auto-scales\n",
    "plt.plot(t[i-iwh:i+iwh,0],d[i-iwh:i+iwh,0],'k-');\n",
    "plt.plot(t[i-iwh:i+iwh,0],d[i-iwh:i+iwh,0],'k*');\n",
    "plt.xlabel('time after Jan 1, 1997 (days)');\n",
    "plt.ylabel('teperature (C)');\n",
    "plt.title('Close-up: Black Rock Forest Temp');\n",
    "plt.show();\n",
    "\n",
    "# whole data figure with all data in grey, selected data in black\n",
    "fig2 = plt.figure(2,figsize=(8,4)); \n",
    "ax2 = plt.subplot(1,1,1);\n",
    "plt.axis([0, 5000, -60, 40])\n",
    "plt.plot(t,d,'k-',linewidth=2,color=(0.8,0.8,0.8));\n",
    "plt.plot(t[i-iwh:i+iwh,0],d[i-iwh:i+iwh,0],'k-');\n",
    "plt.xlabel('Time (days)');\n",
    "plt.ylabel('Teperature (C)');\n",
    "plt.title('Black Rock Forest weather station');\n",
    "plt.show;\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "click at 352.82 with data d[8262]=2.33\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 24160 characters]\n",
      "text/plain": [
       "<Figure size 576x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 23700 characters]\n",
      "text/plain": [
       "<Figure size 576x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# edapy_02_04: enlarged plot via mouse-click\n",
    "# Black Rock Forest temperature data used as an example\n",
    "\n",
    "#  interactive graphics from withing the Jupyter Notebook is a bit dicey\n",
    "# if METHOD=1 fails, try METHOD=2 (or vice-versa)\n",
    "METHOD=1;\n",
    "\n",
    "# get current backend\n",
    "if( METHOD==1):\n",
    "    bkend = matplotlib.get_backend();\n",
    "\n",
    "# load data from text file\n",
    "D = np.genfromtxt('../Data/brf_temp.txt', delimiter='\\t')\n",
    "[N, M]=D.shape;\n",
    "t = eda_cvec( D[:,0] );\n",
    "d = eda_cvec( D[:,1] );\n",
    "\n",
    "#  use interactive plotting\n",
    "if( METHOD==1 ):\n",
    "    matplotlib.use('Qt5Agg');\n",
    "else:\n",
    "    %matplotlib qt5\n",
    "fig1 = plt.figure(1,figsize=(8,4)); # interactive figure\n",
    "ax1 = plt.subplot(1,1,1);\n",
    "plt.axis([0, 5000, -60, 40])\n",
    "plt.plot(t,d,'k-');\n",
    "plt.xlabel('Time (days)');\n",
    "plt.ylabel('Teperature (C)');\n",
    "plt.title('Black Rock Forest weather station');\n",
    "plt.show();\n",
    "Dt = t[2,0]-t[1,0]; # sampling interval\n",
    "wh = 25.0;          # half width of window in days\n",
    "iwh = floor(wh/Dt); # half width of window in samples\n",
    "v=plt.ginput(1); # get (t,d) of 1 mouse click\n",
    "plt.close();     # close interactive figure\n",
    "tc=v[0][0];      # time of click\n",
    "dc=v[0][1];      # temperature of click\n",
    "r = np.where(t>tc); # find first index of t greater than tc\n",
    "k = r[0];\n",
    "i = k[0];\n",
    "print('click at %.2f with data d[%d]=%.2f' % (tc, i, d[i,0]) );\n",
    "\n",
    "# switch back to standard plotting\n",
    "if( METHOD==1 ):\n",
    "    matplotlib.use(bkend);\n",
    "else:\n",
    "    %matplotlib inline\n",
    "\n",
    "fig2 = plt.figure(2,figsize=(8,4)); # close-up figure\n",
    "ax3 = plt.subplot(1,1,1);\n",
    "plt.plot(t[i-iwh:i+iwh,0],d[i-iwh:i+iwh,0],'k-');\n",
    "plt.plot(t[i-iwh:i+iwh,0],d[i-iwh:i+iwh,0],'k*');\n",
    "plt.xlabel('time after Jan 1, 1997 (days)');\n",
    "plt.ylabel('teperature (C)');\n",
    "plt.title('Enlarged Section');\n",
    "plt.show();\n",
    "\n",
    "# whole data figure with all data in grey, selected data in black\n",
    "fig3 = plt.figure(3,figsize=(8,4)); \n",
    "ax2 = plt.subplot(1,1,1);\n",
    "plt.axis([0, 5000, -60, 40])\n",
    "plt.plot(t,d,'k-',linewidth=2,color=(0.8,0.8,0.8));\n",
    "plt.plot(t[i-iwh:i+iwh,0],d[i-iwh:i+iwh,0],'k-');\n",
    "plt.xlabel('Time (days)');\n",
    "plt.ylabel('Teperature (C)');\n",
    "plt.title('Black Rock Forest weather station');\n",
    "plt.show;\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 20196 characters]\n",
      "text/plain": [
       "<Figure size 576x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# edapy_02_05: histogram of data\n",
    "# Black Rock Forest temperature data used as an example\n",
    "\n",
    "# load data from text file\n",
    "D = np.genfromtxt('../Data/brf_temp.txt', delimiter='\\t')\n",
    "[N, M]=D.shape;\n",
    "t = eda_cvec( D[:,0] );\n",
    "d = eda_cvec( D[:,1] );\n",
    "\n",
    "# histogram\n",
    "Lh = 100;            # number of bins in histogram\n",
    "dmin = np.min(d);    # minimum bin\n",
    "dmax = np.max(d);    # maximum bin\n",
    "c, e = np.histogram(d,Lh,(dmin,dmax));  # create histogram\n",
    "Nc = len(c);  # lengths of counts\n",
    "Ne = len(e);  # length of edges\n",
    "counts = eda_cvec(c);            # vector of counts\n",
    "edges = eda_cvec( e[0:Ne-1] );  # vector of edges\n",
    "centers = edges + 0.5*(edges[1,0]-edges[0,0]); # centers\n",
    "\n",
    "# plot histogram\n",
    "fig1 = plt.figure(1,figsize=(8,4));\n",
    "ax1 = plt.subplot(1,1,1);\n",
    "\n",
    "plt.axis([-60, 40, 0, np.max(c)]);\n",
    "plt.plot(centers,counts,'k-');\n",
    "plt.xlabel('temperature (C)');\n",
    "plt.ylabel('counts');\n",
    "plt.title('Histogram of Temperature Measurements');\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 24168 characters]\n",
      "text/plain": [
       "<Figure size 576x576 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# edapy_02_06: plot histogram two ways\n",
    "# Black Rock Forest temperature data used as an example\n",
    "\n",
    "# load data from text file\n",
    "D = np.genfromtxt('../Data/brf_temp.txt', delimiter='\\t')\n",
    "[N, M]=D.shape;\n",
    "t = eda_cvec( np.copy( D[:,0] ) );\n",
    "d = eda_cvec( np.copy( D[:,1] ) );\n",
    "\n",
    "# histogram\n",
    "Lh = 100;         # number of bins\n",
    "dmin = np.min(d); # minimum bin\n",
    "dmax = np.max(d); # maximum bin\n",
    "c, e = np.histogram(d,Lh,(dmin,dmax)); # compute histogram\n",
    "Nc = len(c); # length of couts\n",
    "Ne = len(e); # length of edges\n",
    "counts = eda_cvec(c);          # vector of counts\n",
    "edges = eda_cvec(e[0:Ne-1]);  # vector of edges\n",
    "centers = edges + 0.5*(edges[1,0]-edges[0,0]); # centers\n",
    "\n",
    "# plot\n",
    "fig1 = plt.figure(1,figsize=(8,8));\n",
    "\n",
    "# subplot 1; traditional line plot\n",
    "ax1 = plt.subplot(1, 2, 1);\n",
    "plt.axis([0, np.max(counts), -60, 40]);\n",
    "plt.gca().invert_yaxis();  # y-axis positive down\n",
    "plt.plot(counts,centers,'k-');\n",
    "plt.xlabel('counts');\n",
    "plt.ylabel('temperature (C)');\n",
    "\n",
    "# subplot 1; grey-scale image\n",
    "ax2 = plt.subplot(1, 2, 2);\n",
    "mycmap = cm.gray;             # grey-scale colormap\n",
    "mycmap = mycmap.reversed();   # reverse intensities, so black=max\n",
    "plt.axis([0.0, 1.0, -60, 40]);\n",
    "myvmax = np.max(count);\n",
    "myextent=(0.0,1.0,np.max(centers),np.min(centers));\n",
    "plt.gca().invert_yaxis();  # y-axis positive down\n",
    "plt.imshow( count, cmap=mycmap, vmin=0, vmax=myvmax, extent=myextent, aspect='auto' );\n",
    "plt.colorbar(label='counts'); # draw color bar\n",
    "plt.show();\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 30596 characters]\n",
      "text/plain": [
       "<Figure size 576x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# edapy_02_07, moving-window histogram of data\n",
    "\n",
    "# load temperature data from text fole\n",
    "D = np.genfromtxt('../Data/brf_temp.txt', delimiter='\\t')\n",
    "[N, M]=D.shape;\n",
    "t = eda_cvec(D[:,0]);\n",
    "d = eda_cvec(D[:,1]);\n",
    "Dt = 1.0/24.0;\n",
    "\n",
    "# There are Lw histograms, each offset apart\n",
    "toffset = 30.0;   # offset of 30 days\n",
    "offset = floor(toffset/Dt); # offset in samples\n",
    "Lw = floor( N/offset ); # number of windows\n",
    "Lh = 100;  # number of bins in each histogram\n",
    "tend = toffset*Lw; # last time\n",
    "\n",
    "# an array to hold the histograms\n",
    "Dw = np.zeros((Lh,Lw));\n",
    "\n",
    "# histograms\n",
    "dmin = -60;\n",
    "dmax = 40;\n",
    "for i in range(Lw): # loop over time windows\n",
    "    j=(i-1)*offset; # index of first data in window\n",
    "    k=j+offset;     # index of last data in window\n",
    "    c, e = np.histogram(d[j:k,0],Lh,(dmin,dmax));\n",
    "    Nc = len(c);\n",
    "    counts = eda_cvec(c);\n",
    "    Dw[0:Lh,i:i+1]=counts;\n",
    "    \n",
    "# plot\n",
    "fig1 = plt.figure(1,figsize=(8,4));\n",
    "ax1 = plt.subplot(1,1,1);\n",
    "plt.rcParams['xtick.bottom'] = False;\n",
    "plt.rcParams['xtick.labelbottom'] = False;\n",
    "plt.rcParams['xtick.top'] = True;\n",
    "plt.rcParams['xtick.labeltop'] = True;\n",
    "plt.xlabel('time (days)');\n",
    "plt.ylabel('temperature (C)');\n",
    "ax1.xaxis.set_label_position('top');\n",
    "plt.axis([0, tend, dmin, dmax]);\n",
    "plt.gca().invert_yaxis();\n",
    "\n",
    "# make a greyscale colormap that is never completely white\n",
    "bw = np.zeros((256,4));\n",
    "v = eda_cvec( 0.9*(256 - np.linspace( 0, 255, 256 ))/255 );\n",
    "w = np.ones((256,1));\n",
    "bw[0:256,0:1] = v;\n",
    "bw[0:256,1:2] = v;\n",
    "bw[0:256,2:3] = v;\n",
    "bw[0:256,3:4] = w;\n",
    "bwcmap = ListedColormap(bw);\n",
    "\n",
    "# alternatively, could have ued colormap\n",
    "# bwcmap = cm.gray;\n",
    "# bwcmap = bwcmap.reversed();\n",
    "\n",
    "plt.imshow( Dw, cmap=bwcmap, vmin=0, vmax=np.max(Dw), aspect='auto', extent=(0,tend,dmin,dmax), origin='lower' );\n",
    "plt.colorbar(label='counts');\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 37912 characters]\n",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# edapy_02_08: superimposed plots with different line types\n",
    "\n",
    "# create simple data\n",
    "Dt=1.0;\n",
    "N=511;\n",
    "tmin=0;\n",
    "tmax=tmin+Dt*(N-1);\n",
    "t = eda_cvec( np.linspace( tmin, tmax, N ) );\n",
    "d1=np.sin(pi*t/tmax);\n",
    "d2=np.sin(2.0*pi*t/tmax);\n",
    "\n",
    "# plot the sample data\n",
    "fig1 = plt.figure(1,figsize=(8,8));\n",
    "ax1 = plt.subplot(1,1,1);\n",
    "plt.axis([tmin, tmax, -1.1, 1.1]);\n",
    "plt.plot(t,d1,'k-');\n",
    "plt.plot(t,d2,'k:');\n",
    "plt.xlabel('time (s)');\n",
    "plt.ylabel('data');\n",
    "plt.title('data consisting of sine waves');\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 37644 characters]\n",
      "text/plain": [
       "<Figure size 720x576 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# edapy_02_09: sideways side-by-side plots\n",
    "\n",
    "# time axis with N points, starting at tmin, with increment Dt\n",
    "Dt=1;\n",
    "N=511;\n",
    "tmin=0;\n",
    "tmax=tmin+Dt*(N-1);\n",
    "t = eda_cvec( np.linspace( tmin, tmax, N ) );\n",
    "\n",
    "#  trig functions for sample data\n",
    "d1=np.sin(pi*t/tmax);\n",
    "d2=np.sin(2*pi*t/tmax);\n",
    "\n",
    "# plot the sample data\n",
    "fig1 = plt.figure(1,figsize=(10,8));\n",
    "# change the position of the tics\n",
    "plt.rcParams['xtick.bottom'] = False;\n",
    "plt.rcParams['xtick.labelbottom'] = False;\n",
    "plt.rcParams['xtick.top'] = True;\n",
    "plt.rcParams['xtick.labeltop'] = True;\n",
    "\n",
    "ax1 = plt.subplot(1,2,1); # subplot 1 of d1\n",
    "ax1.xaxis.set_label_position('top');\n",
    "plt.axis([-1.1, 1.1, tmin, tmax, ]);\n",
    "ax1.invert_yaxis();\n",
    "plt.plot(d1, t,'k-'); # plot data d1\n",
    "plt.plot([0,0],[tmin,tmax],'k:'); # plot zero line\n",
    "plt.xlabel('d1');\n",
    "plt.ylabel('time (s)');\n",
    "\n",
    "\n",
    "ax2 = plt.subplot(1,2,2); # subplot 2 of d2\n",
    "ax2.xaxis.set_label_position('top');\n",
    "plt.axis([-1.1, 1.1, tmin, tmax]);\n",
    "ax2.invert_yaxis();\n",
    "plt.plot(d2, t,'k-'); # plot data d2\n",
    "plt.plot([0,0],[tmin,tmax],'k:');  # plot zero line\n",
    "plt.xlabel('d2');\n",
    "plt.ylabel('time (s)');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 8188 characters]\n",
      "text/plain": [
       "<Figure size 576x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 8280 characters]\n",
      "text/plain": [
       "<Figure size 576x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# edapy_02_10: colormaps\n",
    "# note you must: import matplotlib.cm as cm\n",
    "\n",
    "# test image\n",
    "N = 11;\n",
    "D = np.random.uniform( low=0.0, high=1.0, size=(N,N) );\n",
    "\n",
    "# plot with Python's gray colormap\n",
    "fig1 = plt.figure(1,figsize=(8,4));\n",
    "ax1 = plt.subplot(1,1,1);\n",
    "# Python-provided grey-scale colormap\n",
    "graycm = cm.get_cmap('gray', 256);\n",
    "plt.imshow( D, cmap=graycm, vmin=0.0, vmax=1.0 );\n",
    "plt.title('image plotted with Python''s grey colormap');\n",
    "plt.show();\n",
    "\n",
    "# alternaitive is to build a user-defined colormap\n",
    "n=256;\n",
    "bw = np.zeros((n,4));\n",
    "# note that this colormap is reversed; white is low\n",
    "v = eda_cvec( 0.9*(n - np.linspace( 0.0, n-1, n ))/((n-1)) );\n",
    "w = np.ones((n,1));\n",
    "bw[0:n,0:1] = v;\n",
    "bw[0:n,1:2] = v;\n",
    "bw[0:n,2:3] = v;\n",
    "bw[0:n,3:4] = w;\n",
    "usercmap = ListedColormap(bw);\n",
    "\n",
    "# plot\n",
    "fig2 = plt.figure(2,figsize=(8,4));\n",
    "ax1 = plt.subplot(1,1,1);\n",
    "plt.imshow( D, cmap=usercmap, vmin=0.0, vmax=1.0 );\n",
    "plt.title('image plotted with user-defined colormap');\n",
    "plt.show();\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 137,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 31004 characters]\n",
      "text/plain": [
       "<Figure size 576x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# edapy_02_11, example of plotting an image\n",
    "# based on edapy_02_08: running histogram of temperature\n",
    "\n",
    "# load data from text fild\n",
    "D = np.genfromtxt('../Data/brf_temp.txt', delimiter='\\t')\n",
    "[N, M]=D.shape;\n",
    "t = eda_cvec(D[:,0]);\n",
    "d = eda_cvec(D[:,1]);\n",
    "\n",
    "# moving window histogram\n",
    "offset = 1000;\n",
    "Lh = 100;\n",
    "Lw = floor( N/offset );\n",
    "Dw = np.zeros((Lh,Lw));\n",
    "tend = offset*(Lw-1); # last time\n",
    "\n",
    "# make image of histograms\n",
    "dmin = -60;\n",
    "dmax = 40;\n",
    "for i in range(Lw):\n",
    "    j=(i-1)*offset;\n",
    "    k=j+offset;\n",
    "    c, e = np.histogram(d[j:k,0],Lh,(dmin,dmax));\n",
    "    Nc = len(c);\n",
    "    counts = eda_cvec(c);\n",
    "    Dw[:,i]=counts.ravel();\n",
    "    \n",
    "# greyscale colormap\n",
    "bw = np.zeros((256,4));\n",
    "v = eda_cvec( 0.9*(256 - np.linspace( 0, 255, 256 ))/255 );\n",
    "w = np.ones((256,1));\n",
    "bw[0:n,0:1] = v;\n",
    "bw[0:n,1:2] = v;\n",
    "bw[0:n,2:3] = v;\n",
    "bw[0:n,3:4] = w;\n",
    "bwcmap = ListedColormap(bw);\n",
    "\n",
    "# plot image; note image in matrix Dw\n",
    "fig1 = plt.figure(1,figsize=(8,4)); # open figure\n",
    "ax1 = plt.subplot(1,1,1);       # one plot in figure\n",
    "plt.rcParams['xtick.bottom'] = False;  # no ticks on bottom\n",
    "plt.rcParams['xtick.labelbottom'] = False; # no tick labels\n",
    "                                           # on bottom\n",
    "plt.rcParams['xtick.top'] = True; # ticks on top\n",
    "plt.rcParams['xtick.labeltop'] = True; # tick labels on top\n",
    "ax1.xaxis.set_label_position('top'); # axis label on top\n",
    "vlo = 0.0;          # low value of color map\n",
    "vhi = np.max(Dw);   # high value of color map\n",
    "ex=(0,tend,dmin,dmax); # corners of image plot here\n",
    "plt.imshow(Dw, cmap=bwcmap, vmin=vlo, vmax=vhi,\n",
    "  aspect='auto', extent=ex, origin='lower' ); # plot matrix Dw\n",
    "ax1.invert_yaxis();  # vertical axis increases downward\n",
    "plt.colorbar(label='counts');  # plot colrbar with label\n",
    "plt.xlabel('time (days)');  # label horizontal axis\n",
    "plt.ylabel('temperature (C)'); # label verical axis\n",
    "plt.title('moving window histogram'); # title\n",
    "plt.show(); # plot finished figure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 149,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 67144 characters]\n",
      "text/plain": [
       "<Figure size 1152x576 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# edapy_02_12: plot and histgram of rates of change\n",
    "# using the Neuse River Hydrograph data as an example\n",
    "\n",
    "# load data from text fole\n",
    "D = np.genfromtxt('../Data/neuse.txt', delimiter='\\t')\n",
    "[N, M]=D.shape;\n",
    "\n",
    "# time in column 1, discharge in colun 2\n",
    "t = eda_cvec(D[:,0]);\n",
    "d = eda_cvec(D[:,1]);\n",
    "\n",
    "# min, max times for plotting\n",
    "tmin=500;\n",
    "tmax = 1000;\n",
    "\n",
    "# compute rate\n",
    "dddt = np.divide( d[1:N,0:1]-d[0:N-1,0:1],\n",
    "                        t[1:N,0:1]-t[0:N-1,0:1] );\n",
    "\n",
    "# plot the sample data\n",
    "fig1 = plt.figure(1,figsize=(16,8));\n",
    "plt.rcParams['xtick.bottom'] = False;\n",
    "plt.rcParams['xtick.labelbottom'] = False;\n",
    "plt.rcParams['xtick.top'] = True;\n",
    "plt.rcParams['xtick.labeltop'] = True;\n",
    "\n",
    "# subplot 1 of data\n",
    "ax1 = plt.subplot(1,3,1);\n",
    "ax1.xaxis.set_label_position('top');\n",
    "plt.axis([ 0, 20000, tmin, tmax,]);\n",
    "plt.plot(d,t,'k-');\n",
    "ax1.invert_yaxis();\n",
    "plt.xlabel('discharge d');\n",
    "plt.ylabel('time (days)');\n",
    "\n",
    "# subplot 2 of derivative\n",
    "ax2 = plt.subplot(1,3,2);\n",
    "ax2.xaxis.set_label_position('top');\n",
    "plt.axis([-5000, 5000, tmin, tmax, ]);\n",
    "ax2.invert_yaxis();\n",
    "plt.plot(dddt,t[0:N-1],'k-');\n",
    "plt.xlabel('derivative dd/dt');\n",
    "plt.ylabel('time (days)');\n",
    "\n",
    "# histogram\n",
    "Lh = 1000;\n",
    "dmin = np.min(dddt);\n",
    "dmax = np.max(dddt);\n",
    "c, e = np.histogram(dddt,Lh,(dmin,dmax));\n",
    "Dbin = (dmin-dmax)/(Lh-1);\n",
    "Nc = len(c);\n",
    "Ne = len(e);\n",
    "count = eda_cvec(c);\n",
    "edges = eda_cvec(e[0:Ne-1]);\n",
    "centers = edges + Dbin/2.0;\n",
    "\n",
    "# subplot 3 of histogram\n",
    "ax3 = plt.subplot(1, 3, 3);\n",
    "ax3.xaxis.set_label_position('top');\n",
    "plt.axis([0, np.max(c), -500, 500]);\n",
    "ax3.invert_yaxis();\n",
    "plt.plot(count,centers,'k-');\n",
    "plt.plot([0,500],[0,0],'k:');\n",
    "plt.xlabel('counts');\n",
    "plt.ylabel('dddt');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 64888 characters]\n",
      "text/plain": [
       "<Figure size 1008x864 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# edapy_02_13: segregating data based on the sign of rate\n",
    "# using the Neuse River Hydrograph data as an example\n",
    "\n",
    "# load data from text fole\n",
    "D = np.genfromtxt('../Data/neuse.txt', delimiter='\\t')\n",
    "[N, M]=D.shape;\n",
    "\n",
    "# time in column 1, discharge in colun 2\n",
    "t = D[0:N,0:1];\n",
    "d = D[0:N,1:2];\n",
    "\n",
    "# derivative\n",
    "dddt = np.divide( d[1:N,0:1]-d[0:N-1,0:1],\n",
    "                        t[1:N,0:1]-t[0:N-1,0:1] );\n",
    "\n",
    "# divide (d, dddt) pairs into two groups depending\n",
    "# upon whether dddt is positive or negativ\n",
    "r = np.where( dddt>=0 );\n",
    "pos = r[0];\n",
    "r = np.where( dddt<0 );\n",
    "neg = r[0];\n",
    "dpos=d[pos,0:1];\n",
    "dddtpos=dddt[pos,0:1];\n",
    "dneg=d[neg,0:1];\n",
    "dddtneg=dddt[neg,0:1];\n",
    "\n",
    "# plot\n",
    "fig1 = plt.figure(1,figsize=(14,12));\n",
    "plt.rcParams['xtick.bottom'] = True;\n",
    "plt.rcParams['xtick.labelbottom'] = True;\n",
    "plt.rcParams['xtick.top'] = False;\n",
    "plt.rcParams['xtick.labeltop'] = False;\n",
    "\n",
    "# d vs dddt scatter plot for positive rates\n",
    "ax1 = plt.subplot(2,1,1);\n",
    "plt.axis([ 0, 10000, 0, 5000]);\n",
    "plt.plot(dpos,dddtpos,'k.');\n",
    "plt.xlabel('dischagre (cfs)');\n",
    "plt.ylabel('d/dt of discharge');\n",
    "plt.title('positive rates');\n",
    "\n",
    "# d vs dddt scatter plot for negative rates\n",
    "ax2 = plt.subplot(2,1,2);\n",
    "plt.axis([ 0, 10000, -5000, 0]);\n",
    "plt.plot(dneg,dddtneg,'k.');\n",
    "plt.xlabel('dischagre (cfs)');\n",
    "plt.ylabel('d/dt of discharge');\n",
    "plt.title('negative rates');\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 166,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 276852 characters]\n",
      "text/plain": [
       "<Figure size 1152x2016 with 28 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# edapy_02_14: scatter plots of Atlantic Rock chemical data\n",
    "\n",
    "# load data from text fole\n",
    "D = np.genfromtxt('../Data/rocks.txt', delimiter='\\t')\n",
    "[N, M]=D.shape;\n",
    "\n",
    "# note the use of a list to hold names\n",
    "names = [ 'SiO2', 'TiO2', 'Al203', 'FeO-total',\n",
    "           'MgO', 'CaO', 'Na2O', 'K2O'];\n",
    "Ns = len(names);\n",
    "\n",
    "# Np: number of permulations of Ns things taken two at a time\n",
    "Np = Ns*Ns-(Ns*(Ns+1)/2);\n",
    "M = 4; # plots in a row\n",
    "Nr = ceil(Np/M); # number of rows\n",
    "i=2;\n",
    "\n",
    "fig1 = plt.figure(1,figsize=(4*M,4*Nr));\n",
    "sp = 0;  # counts plots\n",
    "for i in range(0,Ns-1):\n",
    "    for j in range(i+1,Ns):\n",
    "           sp = sp+1; # increment counter\n",
    "           ax1 = plt.subplot(Nr, M, sp);\n",
    "           plt.plot(D[:,i],D[:,j],'k.');\n",
    "           plt.xlabel(names[i]);\n",
    "           plt.ylabel(names[j]);\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 176,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Bill\n",
      "C:/bill\n",
      " \n",
      "number of names: 4\n",
      "name 0: lithium\n",
      "name 1: calcium\n",
      "name 2: sodium\n",
      "name 3: potassium\n",
      " \n",
      "element 2\n",
      "a=4.72\n",
      "row 2 column 4\n",
      "(2,4): 4.721\n",
      " \n",
      "name 0 is lithium\n",
      "name 1 is calcium\n",
      "name 2 is sodium\n",
      "name 3 is potassium\n"
     ]
    }
   ],
   "source": [
    "# edapy_02_15, character strings\n",
    "\n",
    "# Part 1.  String variables\n",
    "my_name = \"Bill\";\n",
    "my_path = \"C:/bill\";\n",
    "print(my_name);\n",
    "print(my_path);\n",
    "print(\" \");\n",
    "\n",
    "# Part 2.  Using lists to hold strings\n",
    "names = [ 'lithium', 'calcium', 'sodium', 'potassium' ];\n",
    "Ns = len(names);\n",
    "\n",
    "print(\"number of names:\", Ns)\n",
    "print(\"name 0:\", names[0] );\n",
    "print(\"name 1:\", names[1] );\n",
    "print(\"name 2:\", names[2] );\n",
    "print(\"name 3:\", names[3] );\n",
    "print(\" \");\n",
    "\n",
    "# Part 3. Formatting strings using the % operator\n",
    "\n",
    "# prints \"element 2\"\n",
    "i=2;\n",
    "mystring1 = \"element %d\" % i;\n",
    "print(mystring1);\n",
    "\n",
    "# prints \"a=4.72\"\n",
    "a=4.7213;\n",
    "mystring2 = \"a=%.2f\" % a;\n",
    "print(mystring2);\n",
    "\n",
    "# prints \"row 2 column 4\"\n",
    "i=2;\n",
    "j=4;\n",
    "mystring3 = \"row %d column %d\" % (i, j);\n",
    "print(mystring3);\n",
    "\n",
    "# prints \"(2,4): 4.721\"\n",
    "i=2;\n",
    "j=4;\n",
    "a=4.7213;\n",
    "mystring4 = \"(%d,%d): %.3f\" % (i, j, a);\n",
    "print(mystring4);\n",
    "print(\" \");\n",
    "\n",
    "# Part 4. Combining lists and formats\n",
    "for i in range(Ns):\n",
    "    print(\"name %d is %s\" % (i, names[i]));\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
