hovmoller_joint.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now combining the hovmoller plots for the two different sets of time periods into a single joint figure, following the previous template used in banded_and_local_changes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from bootstrap import *\n",
    "from chr_matplotlib import remappedColorMap\n",
    "from jd_matplotlib_tools import add_double_calendar_axis, mark_time_periods\n",
    "import matplotlib.gridspec as gridspec\n",
    "import matplotlib.pyplot as plt\n",
    "import netCDF4 as nc\n",
    "import numpy as np\n",
    "import numpy.ma as ma\n",
    "import os\n",
    "from shutil import copy\n",
    "import time\n",
    "\n",
    "%load_ext autotime"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "time: 4.36 ms\n"
     ]
    }
   ],
   "source": [
    "ferretpath = \"/Users/Siwen/Desktop/ferret/bin/\"\n",
    "ff = \"Pchina_type_notaiwan.nc\"\n",
    "filename = ferretpath + ff\n",
    "myf = nc.Dataset(filename, 'r')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "time: 51.5 ms\n"
     ]
    }
   ],
   "source": [
    "P_total =  myf.variables['Pchina_total'][:]\n",
    "P_band = myf.variables['Pchina_band'][:]\n",
    "P_local = myf.variables['Pchina_local'][:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "time: 21.7 ms\n"
     ]
    }
   ],
   "source": [
    "#years of interest - start and end date, with index starting at 1 out of habit.\n",
    "p1_yr_b = [1,30]\n",
    "p1_yr_e = [29,43]\n",
    "p2_yr_b = [30,44]\n",
    "p2_yr_e = [57,57]\n",
    "\n",
    "P_total_diff = []\n",
    "P_band_diff = []\n",
    "P_local_diff = []\n",
    "\n",
    "#cycle through all time periods of interest and find changes in each type of rainfall b/w periods\n",
    "for p1_b, p1_e, p2_b, p2_e, i in zip(p1_yr_b, p1_yr_e, p2_yr_b, p2_yr_e, np.arange(len(p1_yr_b))): \n",
    "    \n",
    "    P_total_diff.append(np.mean(P_total[:,:,p2_b-1:p2_e], axis=2) - \\\n",
    "                         np.mean(P_total[:,:,p1_b-1:p1_e], axis=2))\n",
    "    P_band_diff.append(np.mean(P_band[:,:,p2_b-1:p2_e], axis=2) \\\n",
    "                       - np.mean(P_band[:,:,p1_b-1:p1_e], axis=2))\n",
    "    P_local_diff.append(np.mean(P_local[:,:,p2_b-1:p2_e], axis=2) \\\n",
    "                        - np.mean(P_local[:,:,p1_b-1:p1_e], axis=2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "time: 7.08 ms\n"
     ]
    }
   ],
   "source": [
    "def smooth(P,dayrange,latrange):\n",
    "    Psmooth=np.zeros(P.shape)\n",
    "    \n",
    "    ll=int((dayrange-1)/2)\n",
    "    yy=int((latrange-1)/2)\n",
    "    \n",
    "    days = P.shape[0]\n",
    "    lats = P.shape[1]\n",
    "    \n",
    "    for d in range(days):\n",
    "        for j in range(lats):\n",
    "            \n",
    "            Psample = P.take(range(d-ll,d+ll+1), mode='wrap', axis=0)\n",
    "            Ps = Psample[:,max(0,j-yy):min(days,j+yy+1)]\n",
    "            Psmooth[d,j] = np.mean(np.mean(Ps))\n",
    "    \n",
    "    return Psmooth"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "time: 4.54 s\n"
     ]
    }
   ],
   "source": [
    "## SMOOTHING ##\n",
    "#smooth out the decadal precipitation means\n",
    "#but, no latitudinal smoothing since rainfall already effectively smoothed in latitude\n",
    "daysmth = 15\n",
    "latsmth = 0\n",
    "\n",
    "#one line! thank you python!\n",
    "P_total_diff_smth = [ smooth(x,daysmth,latsmth) for x in P_total_diff ]\n",
    "P_band_diff_smth = [ smooth(x,daysmth,latsmth) for x in P_band_diff ]\n",
    "P_local_diff_smth = [ smooth(x,daysmth,latsmth) for x in P_local_diff ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Data/chinarain_diff_pval_notaiwan_5179_80107_total_2day_2000iter.nc successfully found\n",
      "Data/chinarain_diff_pval_notaiwan_5179_80107_band_2day_2000iter.nc successfully found\n",
      "Data/chinarain_diff_pval_notaiwan_5179_80107_local_2day_2000iter.nc successfully found\n",
      "Data/chinarain_diff_pval_notaiwan_8093_94107_total_2day_2000iter.nc successfully found\n",
      "Data/chinarain_diff_pval_notaiwan_8093_94107_band_2day_2000iter.nc successfully found\n",
      "Data/chinarain_diff_pval_notaiwan_8093_94107_local_2day_2000iter.nc successfully found\n",
      "time: 59.1 ms\n"
     ]
    }
   ],
   "source": [
    "## LOAD COMPUTED P-VALUES ##\n",
    "\n",
    "#which set of sim values are we trying to access?\n",
    "pval_root = 'Data/chinarain_diff_pval_notaiwan_'\n",
    "blklen = 2\n",
    "niter = 2000\n",
    "\n",
    "P_total_diff_pval = []\n",
    "P_band_diff_pval = []\n",
    "P_local_diff_pval = []\n",
    "\n",
    "#cycle through all time periods of interest and find changes in each type of rainfall b/w periods\n",
    "for p1_b, p1_e, p2_b, p2_e, i in zip(p1_yr_b, p1_yr_e, p2_yr_b, p2_yr_e, np.arange(len(p1_yr_b))): \n",
    "    \n",
    "    #LOAD P-VALUES FOR TOTAL RAINFALL\n",
    "    raintype = 'total'\n",
    "    yrs = str(p1_b+50) + str(p1_e+50) + '_' + str(p2_b+50) + str(p2_e+50) + '_'\n",
    "    pval_filename = pval_root + yrs + raintype + '_' + str(blklen) + 'day_' \\\n",
    "        + str(niter) + 'iter.nc'\n",
    "    \n",
    "    if os.path.isfile(pval_filename):\n",
    "        print(pval_filename + ' successfully found')\n",
    "        myf = nc.Dataset(pval_filename, 'r')\n",
    "        P_total_diff_pval.append(myf.variables['pval'][:])\n",
    "        myf.close()\n",
    "        \n",
    "    else:\n",
    "        print(pval_filename + ' does not exist')\n",
    "    \n",
    "    \n",
    "    ##LOAD P-VALUES FOR BANDED RAINFALL\n",
    "    raintype = 'band'\n",
    "    yrs = str(p1_b+50) + str(p1_e+50) + '_' + str(p2_b+50) + str(p2_e+50) + '_'\n",
    "    pval_filename = pval_root + yrs + raintype + '_' + str(blklen) + 'day_' \\\n",
    "        + str(niter) + 'iter.nc'\n",
    "    \n",
    "    if os.path.isfile(pval_filename):\n",
    "        print(pval_filename + ' successfully found')\n",
    "        myf = nc.Dataset(pval_filename, 'r')\n",
    "        P_band_diff_pval.append(myf.variables['pval'][:])\n",
    "        myf.close()\n",
    "        \n",
    "    else:\n",
    "        print(pval_filename + ' does not exist')\n",
    "            \n",
    "    ##LOAD P-VALUES FOR LOCAL RAINFALL\n",
    "    raintype = 'local'\n",
    "    yrs = str(p1_b+50) + str(p1_e+50) + '_' + str(p2_b+50) + str(p2_e+50) + '_'\n",
    "    pval_filename = pval_root + yrs + raintype + '_' + str(blklen) + 'day_' \\\n",
    "        + str(niter) + 'iter.nc'\n",
    "    \n",
    "    if os.path.isfile(pval_filename):\n",
    "        print(pval_filename + ' successfully found')\n",
    "        myf = nc.Dataset(pval_filename, 'r')\n",
    "        P_local_diff_pval.append(myf.variables['pval'][:])\n",
    "        myf.close()\n",
    "        \n",
    "    else:\n",
    "        print(pval_filename + ' does not exist')    \n",
    "        \n",
    "    #P_total_diff.append()\n",
    "    #P_band_diff.append(np.mean(P_band[:,:,p2_b-1:p2_e], axis=2) \\\n",
    "    #P_local_diff.append(np.mean(P_local[:,:,p2_b-1:p2_e], axis=2) \\"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "time: 16.1 ms\n"
     ]
    }
   ],
   "source": [
    "## CREATE MASKED ARRAYS - only show values that are statistically significant\n",
    "#create two masked arrays - one showing 95% confidence, one showing 99% confidence\n",
    "\n",
    "##two separate for loops: first one creates masks, second one produces masked rainfall.\n",
    "#reason for doing things this way is that masks can then be accessed independently at a later point in time.\n",
    "\n",
    "#not super-pythonic method, but at least it makes it easy to understand.\n",
    "\n",
    "mk_total = [] #storage format will be a LIST of 2-element tuples - each set of tuples is the 95% and 99% mask respectively\n",
    "mk_band = []\n",
    "mk_local = []\n",
    "\n",
    "#set up 95th and 99th percentile masks for each type of rainfall\n",
    "for P_total_diff, P_band_diff, P_local_diff, P_total_pval, P_band_pval, P_local_pval in \\\n",
    "    zip(P_total_diff_smth, P_band_diff_smth, P_local_diff_smth, P_total_diff_pval, P_band_diff_pval, P_local_diff_pval):\n",
    "    \n",
    "    mk_total_95 = abs(P_total_pval-.5) < .475 #95th percentile mask for rainfall changes\n",
    "    mk_total_99 = abs(P_total_pval-.5) < .495 #99th percentile mask for rainfall changes\n",
    "    mk_total.append([mk_total_95, mk_total_99])\n",
    "    \n",
    "    mk_band_95 = abs(P_band_pval-.5) < .475 #95th percentile mask for rainfall changes\n",
    "    mk_band_99 = abs(P_band_pval-.5) < .495 #99th percentile mask for rainfall changes\n",
    "    mk_band.append([mk_band_95, mk_band_99])\n",
    "    \n",
    "    mk_local_95 = abs(P_local_pval-.5) < .475 #95th percentile mask for rainfall changes\n",
    "    mk_local_99 = abs(P_local_pval-.5) < .495 #99th percentile mask for rainfall changes\n",
    "    mk_local.append([mk_local_95, mk_local_99])\n",
    "    \n",
    "#create masked arrays\n",
    "P_total_diff_mask = []\n",
    "P_band_diff_mask = []\n",
    "P_local_diff_mask = []\n",
    "\n",
    "#create masked arrays with significant rainfall changes only based on 95th percentile masks\n",
    "for P_total_diff, P_band_diff, P_local_diff, mask_total, mask_band, mask_local in \\\n",
    "    zip(P_total_diff_smth, P_band_diff_smth, P_local_diff_smth, mk_total, mk_band, mk_local):\n",
    "    \n",
    "    P_total_diff_mask.append(ma.array(P_total_diff, mask = mask_total[0]))\n",
    "    P_band_diff_mask.append(ma.array(P_band_diff, mask = mask_band[0]))\n",
    "    P_local_diff_mask.append(ma.array(P_local_diff, mask = mask_local[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "5179\n",
      "8007\n",
      "8093\n",
      "9407\n",
      "time: 51.7 ms\n"
     ]
    }
   ],
   "source": [
    "## LOAD 1980-2007 v 1951-1979 changes FROM NETCF ##\n",
    "\n",
    "## NETCDF FILENAMES\n",
    "freq_diff_file = \"RDA_freq_diff_tau_seasonal\"\n",
    "freq_diff_name = ferretpath + freq_diff_file + '.nc'\n",
    "\n",
    "int_diff_file = \"RDA_int_diff\"\n",
    "int_diff_name = ferretpath + int_diff_file + '.nc'\n",
    "\n",
    "int_pval_file = \"RDA_int_diff_bootstrap_2000iter_perm\"\n",
    "int_pval_name = ferretpath + int_pval_file + '.nc'\n",
    "\n",
    "\n",
    "\n",
    "#import from netCDF file and load 57-year climatology\n",
    "freq_diff_netcdf = nc.Dataset(freq_diff_name, \"r\")\n",
    "int_diff_netcdf = nc.Dataset(int_diff_name, \"r\")\n",
    "int_pval_netcdf = nc.Dataset(int_pval_name, \"r\")\n",
    "\n",
    "#limits of different time periods that we're comparing\n",
    "p1_beg = [1951,1980]\n",
    "p1_end = [1979,1993]\n",
    "p2_beg = [1980,1994]\n",
    "p2_end = [2007,2007]\n",
    "\n",
    "#in this case, we're just going to load the smoothed plots of hovmoller frequency\n",
    "freq_diff = []\n",
    "freq_pval = []\n",
    "int_diff = []\n",
    "int_pval = []\n",
    "\n",
    "for p1_b, p1_e, p2_b, p2_e in zip(p1_beg, p1_end, p2_beg, p2_end):\n",
    "\n",
    "    years_p1= str(p1_b-1900) + str(\"{0:0=2d}\".format((p1_e-1900)%100))\n",
    "    print(years_p1)\n",
    "    years_p2= str(p2_b-1900) + str(\"{0:0=2d}\".format((p2_e-1900)%100))\n",
    "    print(years_p2)\n",
    "    \n",
    "    freq_diff.append(freq_diff_netcdf['lat_diff_tau_' + years_p2 + '_' + years_p1][:])\n",
    "    freq_pval.append(freq_diff_netcdf['lat_pval_tau_' + years_p2 + '_' + years_p1][:])\n",
    "    int_diff.append(int_diff_netcdf['int_diff_' + years_p2 + '_' + years_p1][:])\n",
    "    int_pval.append(int_pval_netcdf['int_pval_' + years_p2 + '_' + years_p1][:])\n",
    "    \n",
    "#close it\n",
    "freq_diff_netcdf.close()\n",
    "int_diff_netcdf.close()\n",
    "int_pval_netcdf.close()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "time: 9.19 ms\n"
     ]
    }
   ],
   "source": [
    "## USE P-VALS TO CREATE MASKS APPLIED TO ONLY SHOW SIGNIFICANT DIFFERENCES\n",
    "## 1980-2007 versus 1951-1979 - clearly cuts out a lot of random noise!\n",
    "\n",
    "## frequency masking by analytic p-value\n",
    "freq_diff_mask_95 = []\n",
    "freq_diff_mask_99 = []\n",
    "\n",
    "for fq, pval in zip(freq_diff, freq_pval):\n",
    "    freq_diff_mask_95.append(ma.array(fq, mask = ~(abs(pval-.5) > .475)))\n",
    "    freq_diff_mask_99.append(ma.array(fq, mask = ~(abs(pval-.5) > .495)))\n",
    "\n",
    "#intensity masking by bootstrapped p-value\n",
    "int_diff_mask_95 = []\n",
    "int_diff_mask_99 = []   \n",
    "    \n",
    "for ints, pval in zip(int_diff, int_pval):\n",
    "    int_diff_mask_95.append(ma.array(ints, mask = ~(abs(pval-.5) > .475)))\n",
    "    int_diff_mask_99.append(ma.array(ints, mask = ~(abs(pval-.5) > .495)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In retrospect, the way that I did the rest of this code is not up to par, but hindsight is 20/20"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "time: 7.78 ms\n"
     ]
    }
   ],
   "source": [
    "## Set up variables that will be plotted in correct order\n",
    "P_plot_tuples = {'a':[P_total_diff_mask[0].T,P_band_diff_mask[0].T,100*freq_diff_mask_95[0][:,1:21].T,\\\n",
    "                   int_diff_mask_95[0][:,1:21].T],\n",
    "                 'b':[P_total_diff_mask[1].T,P_band_diff_mask[1].T,100*freq_diff_mask_95[1][:,1:21].T,\\\n",
    "                   int_diff_mask_95[1][:,1:21].T]}\n",
    "\n",
    "P_contour_tuples = {'a':[mk_total[0][1].T,mk_band[0][1].T,freq_diff_mask_99[0][:,1:21].mask.T,\\\n",
    "                         int_diff_mask_99[0][:,1:21].mask.T],\n",
    "                    'b':[mk_total[1][1].T,mk_band[1][1].T,freq_diff_mask_99[1][:,1:21].mask.T,\\\n",
    "                         int_diff_mask_99[1][:,1:21].mask.T]}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Row 1\n",
      "2.14082144474 5.03482786761 -2.22076927579 -2.96129301969\n",
      "Row 2\n",
      "1.72500510079 5.19934700677 -2.65746721782 -2.69964654647\n",
      "Row 3\n",
      "2.14942528736 3.52380952381 -3.09031198686 -2.66666666667\n",
      "Row 4\n",
      "15.4565939207 11.7522309632 -8.97392758901 -45.5776983787\n",
      "time: 8.07 ms\n"
     ]
    }
   ],
   "source": [
    "for i in range(4):\n",
    "    a = P_plot_tuples['a'][i].max()\n",
    "    b = P_plot_tuples['b'][i].max()\n",
    "    c = P_plot_tuples['a'][i].min()\n",
    "    d = P_plot_tuples['b'][i].min()\n",
    "    print('Row',i+1)\n",
    "    print(a,b,c,d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/Siwen/anaconda/lib/python3.6/site-packages/matplotlib/text.py:2138: UserWarning: You have used the `textcoords` kwarg, but not the `xytext` kwarg.  This can lead to surprising results.\n",
      "  warnings.warn(\"You have used the `textcoords` kwarg, but not \"\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "time: 1.83 s\n"
     ]
    }
   ],
   "source": [
    "## New version of Hovmoller plots of difference, based off of banded_and_local_changes\n",
    "axlims = [2,2,2,6]\n",
    "palettes = ['bwr','bwr','bwr','bwr']\n",
    "yticklabels = [\"20$\\!^\\circ\\!$N\",\"24$\\!^\\circ$\",\"28$\\!^\\circ$\",\"32$\\!^\\circ$\",\"36$\\!^\\circ$\",\"40$\\!^\\circ\\!$N\"]\n",
    "## plots a-d are at .25 degree latitudinal resolution, e-h are all at 1 degree rez\n",
    "ytickpos = [np.linspace(0,80,6),np.linspace(0,80,6),np.linspace(0,20,6),np.linspace(0,20,6)]\n",
    "\n",
    "## Define time periods of interest - these will be demarcated on the plots with vertical lines\n",
    "time_periods = [[60.5,120.5],[120.5,160.5],[160.5,200.5],[200.5,273.5],[273.5,320.5]]\n",
    "\n",
    "fig = plt.figure(figsize=(20,10), facecolor='w')\n",
    "#st = plt.suptitle(\"Changes in Rainfall by Type\", size = 32, y=1, x=.515)\n",
    "\n",
    "## create grids for left-hand and right-hand sides\n",
    "grid_left = gridspec.GridSpec(4, 1)\n",
    "grid_left.update(left=0.05, right=0.49, hspace = .6)\n",
    "grid_right = gridspec.GridSpec(4, 1)\n",
    "grid_right.update(left=0.54, right=0.98, hspace = .6)\n",
    "grids = {'a':grid_left, 'b':grid_right}\n",
    "\n",
    "## create storage dict for axes objects for each side; a is left, b is right\n",
    "axes = {'a':[], 'b':[]}\n",
    "plots = {'a':[], 'b':[]}\n",
    "\n",
    "## Define title for each set of subplots; first row a-b, 2nd c-d and so forth\n",
    "subplot_titles={'a':[\"a) Change in Total Rainfall (mm day$^{-1}$)\",\"c) Change in Frontal Rainfall (mm day$^{-1}$)\",\\\n",
    "                \"e) Change in Frontal Rain Event Frequency (%)\",\"g) Change in Event Daily Accumulation (mm day$^{-1}$)\"],\n",
    "                'b':[\"b) Change in Total Rainfall (mm day$^{-1}$)\",\"d) Change in Frontal Rainfall (mm day$^{-1}$)\",\\\n",
    "                \"f) Change in Frontal Rain Event Frequency (%)\",\"h) Change in Event Daily Accumulation (mm day$^{-1}$)\"]}\n",
    "\n",
    "## Define title for each half of the figure\n",
    "panel_titles = {'a':'1980-2007 v 1951-1979','b':'1994-2007 v 1980-1993'}\n",
    "\n",
    "## Biggest change from previous version is that instead of explicitly\n",
    "## enumerating each subplot, now line everything up at the beginning and\n",
    "## then iterate through each subplot with an identical loop\n",
    "for half in ['a','b']:\n",
    "    for row,(precip,contour,title,ytickposition,axlim,palette) in enumerate(zip(P_plot_tuples[half],P_contour_tuples[half],\\\n",
    "                                                  subplot_titles[half],ytickpos,axlims,palettes)):\n",
    "        ax = plt.subplot(grids[half][row,0])\n",
    "        axes[half].append(ax)\n",
    "        plot = ax.pcolormesh(precip, cmap = palette, vmin = -axlim, vmax = axlim); \n",
    "        plots[half].append(plot)\n",
    "        ax.contour(contour,1,colors = 'gray')\n",
    "        ax.set_title(title, fontsize=20, x=0, y = 1.14, horizontalalignment = 'left')\n",
    "\n",
    "        ## only mark latitude axis\n",
    "        if half == 'a':\n",
    "            ax.set_ylabel('Latitude ($\\!^\\circ\\!$N)', fontsize=12, labelpad=0)\n",
    "\n",
    "        ## swap side of latitude tickmarks to make room for colorbars\n",
    "        if half == 'b':\n",
    "            ax.yaxis.tick_right()\n",
    "        \n",
    "        ax.set_yticks(ytickposition)\n",
    "        ax.set_yticklabels(yticklabels)\n",
    "        add_double_calendar_axis()\n",
    "        mark_time_periods(time_periods, numbers = 'yes', box = 'no')\n",
    "\n",
    "    ## add title to this half - place right above middle panel in top row\n",
    "    pos = axes[half][0].get_position() # get the original position \n",
    "    xx = (pos.x0+pos.x1)/2\n",
    "    yy = pos.y1+.07\n",
    "    plt.figtext(xx,yy,panel_titles[half],horizontalalignment='center',size=24,weight='bold')\n",
    "    \n",
    "    \n",
    "#########################\n",
    "## COLOR BARS\n",
    "## Experimenting with a new layout for color bars - attempting to use 3 color bars, one for\n",
    "## top 4 panels, one for next 2, one for bottom 2\n",
    "\n",
    "cbar_axes = []\n",
    "\n",
    "## Color bar 1 - top 4 panels\n",
    "cbar_ax_1 = fig.add_axes([0.50, 0.5883, 0.01, 0.25])\n",
    "cbar_axes.append(cbar_ax_1)\n",
    "cb_ticks = [ int(x) for x in np.linspace(-axlims[0],axlims[0],5) ]\n",
    "cb1 = fig.colorbar(plots['a'][0],cax=cbar_ax_1, ticks=cb_ticks, orientation='vertical')\n",
    "cb_ticklabels = ['<-2',int(-1),int(0),int(1),'>2']\n",
    "cb1.ax.set_yticklabels(cb_ticklabels,size=18)  # horizontal colorbar\n",
    "cb1.ax.tick_params(length=0)\n",
    "\n",
    "subplot_pos = axes['a'][2].get_position().get_points()\n",
    "subplot_ylen = subplot_pos[1][1]-subplot_pos[0][1]\n",
    "\n",
    "## Color bar 2 - panels e and f\n",
    "cbar2_ybottom = axes['a'][2].get_position().get_points()[0][1]\n",
    "cbar_ax_2 = fig.add_axes([0.50, cbar2_ybottom, 0.01, subplot_ylen])\n",
    "cbar_axes.append(cbar_ax_2)\n",
    "cb_ticks = [ int(x) for x in np.linspace(-axlims[2],axlims[2],5) ]\n",
    "cb2 = fig.colorbar(plots['a'][2],cax=cbar_ax_2, ticks=cb_ticks, orientation='vertical')\n",
    "cb_ticklabels = ['<-2',int(-1),int(0),int(1),'>2']\n",
    "cb2.ax.set_yticklabels(cb_ticklabels,size=16)  # horizontal colorbar\n",
    "cb2.ax.tick_params(length=0)\n",
    "\n",
    "## Color bar 3 - bottom 2 plots\n",
    "cbar3_ybottom = axes['a'][3].get_position().get_points()[0][1]\n",
    "cbar_ax_3 = fig.add_axes([0.50, cbar3_ybottom, 0.01, subplot_ylen])\n",
    "cbar_axes.append(cbar_ax_3)\n",
    "cb_ticks = [ int(x) for x in np.linspace(-axlims[3],axlims[3],5) ]\n",
    "cb3 = fig.colorbar(plots['a'][3],cax=cbar_ax_3, ticks=cb_ticks, orientation='vertical')\n",
    "cb_ticklabels = ['<-6',int(-3),int(0),int(3),'>6']\n",
    "cb3.ax.set_yticklabels(cb_ticklabels,size=16)  # horizontal colorbar\n",
    "cb3.ax.tick_params(length=0)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting existing version of file in RDA manuscript directory...\n",
      "Overwriting existing version of file in RDA Python directory...\n",
      "time: 1min 23s\n"
     ]
    }
   ],
   "source": [
    "## SAVE COPIES OF FIGURE TO MANUSCRIPT AND ANALYSIS DIRECTORIES\n",
    "figure_name = \"hov_summary.pdf\"\n",
    "\n",
    "savepath = \"/Users/Siwen/Desktop/RDA/Figures/\" + figure_name\n",
    "\n",
    "if os.path.exists(savepath):\n",
    "    print(\"Overwriting existing version of file in RDA manuscript directory...\")\n",
    "    os.remove(savepath)\n",
    "    \n",
    "plt.savefig(savepath, bbox_inches='tight')\n",
    "    \n",
    "copypath = \"/Users/Siwen/RDA/Figures/\" + figure_name\n",
    "    \n",
    "if os.path.exists(copypath):\n",
    "    print(\"Overwriting existing version of file in RDA Python directory...\")\n",
    "    os.remove(copypath)\n",
    "    \n",
    "copy(savepath,copypath)\n",
    "\n",
    "plt.savefig(savepath, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}