notebooks/CO2_Flux_Decomposition/co2_flux_decomposition.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# CO$_{2}$ Flux Decomposition\n",
    "\n",
    "This is the working area to develop a script that takes time series from a given region and decomposes the CO2 flux changes from a given climate event (e.g. ENSO) into the relative components contributing to that change (U, SST, S, DIC, etc.)\n",
    "\n",
    "To improve this script, can compute sensitivities in 5-10 year chunks to reduce error due to changing mean state."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import xarray as xr\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "plt.style.use('seaborn-whitegrid')\n",
    "import seaborn as sns\n",
    "from esmtools.physics import stress_to_speed\n",
    "from esmtools.carbon import co2_sol, schmidt\n",
    "import cartopy.crs as ccrs\n",
    "import cartopy.feature as cfeature\n",
    "import esmtools as et"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Some constants that are pertinent for this.\n",
    "\"\"\"\n",
    "a = 6.972e-7 # Flux coefficient in s/m\n",
    "conversion = (60 * 60 * 24 * 365) / 1e3 # Applied to every term of the expansion to have proper units (see notes)\n",
    "S0 = 35 # Salinity constant"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "e = 'CalCS'\n",
    "x = 'NPGO'\n",
    "smooth = 12"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Load in all important ensemble means\n",
    "\n",
    "Will use this to compute sensitivities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def load_ensemble_mean(e, v):\n",
    "    \"\"\"\n",
    "    Loads in the ensemble mean time series for a given upwelling system (e)\n",
    "    and variable (v). I will use these to compute the sensitivities, which shouldn't\n",
    "    cause too much error since it's a much larger value than the residuals.\n",
    "    \"\"\"\n",
    "    filepath = ('/glade/p/work/rbrady/EBUS_BGC_Variability/' + v + '/' + e +\n",
    "               '/filtered_output/' + e.lower() + '-' + v + '-forced-signal-AW-chavez-800km.nc')\n",
    "    da = xr.open_dataarray(filepath)\n",
    "    da.name = v\n",
    "    return da"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "variables = ['FG_CO2', 'SST', 'SALT', 'DIC', 'ALK', 'TAUX',\n",
    "             'TAUY', 'TAUX2', 'TAUY2', 'sDIC', 'sALK', 'U',\n",
    "             'DCO2STAR', 'RHO', 'pCO2SURF']\n",
    "forced = xr.Dataset()\n",
    "for i, name in enumerate(variables):\n",
    "    da = load_ensemble_mean(e, name)\n",
    "    forced = xr.merge([forced, da])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "forced = forced.rename({'pCO2SURF': 'pCO2'})\n",
    "means = forced.mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Load in the regression results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def load_regression(e, x, y, smooth):\n",
    "    \"\"\"\n",
    "    Loads in the area-weighted regression results for a given upwelling system\n",
    "    (e), predictor variable (x), and dependent variable (y). Need to add functionality\n",
    "    for lagging. Currently assumes values of 0 for each of those features.\n",
    "    \"\"\"\n",
    "    filepath = ('/glade/p/work/rbrady/EBUS_BGC_Variability/area_weighted_regional_regressions/' +\n",
    "                e + '/' + y + '/' + x + '/')\n",
    "    if smooth == 0:\n",
    "        filename = (x + '.' + y + '.' + e + '.unsmoothed.area_weighted_regional_regression.lag0.nc')\n",
    "    else:\n",
    "        filename = (x + '.' + y + '.' + e + '.smoothed' + str(smooth) + '.area_weighted_regional_regression.lag0.nc')\n",
    "    ds = xr.open_dataset(filepath + filename)\n",
    "    return ds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "delta_FG_CO2 = load_regression(e, x, 'FG_CO2', smooth)\n",
    "delta_U = load_regression(e, x, 'U', smooth)\n",
    "delta_SST = load_regression(e, x, 'SST', smooth)\n",
    "delta_SALT = load_regression(e, x, 'SALT', smooth)\n",
    "delta_sDIC = load_regression(e, x, 'sDIC', smooth)\n",
    "delta_sALK = load_regression(e, x, 'sALK', smooth)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "tolerance = 0.1\n",
    "np.where((delta_FG_CO2['p'] <= tolerance) & \n",
    "         (delta_U['p'] <= tolerance) &\n",
    "         (delta_SST['p'] <= tolerance) &\n",
    "         (delta_SALT['p'] <= tolerance) &\n",
    "         (delta_sDIC['p'] <= tolerance) &\n",
    "         (delta_sALK['p'] <= tolerance))\n",
    "\"\"\"\n",
    "# Seems impossible to get all 34 members to be significant in all terms. So maybe just\n",
    "# test significance of FG_CO2?\n",
    "tolerance = 0.1\n",
    "delta_U = delta_U.where(delta_FG_CO2['p'] <= tolerance, drop=True)\n",
    "delta_SST = delta_SST.where(delta_FG_CO2['p'] <= tolerance, drop=True)\n",
    "delta_SALT = delta_SALT.where(delta_FG_CO2['p'] <= tolerance, drop=True)\n",
    "delta_sDIC = delta_sDIC.where(delta_FG_CO2['p'] <= tolerance, drop=True)\n",
    "delta_sALK = delta_sALK.where(delta_FG_CO2['p'] <= tolerance, drop=True)\n",
    "delta_FG_CO2 = delta_FG_CO2.where(delta_FG_CO2['p'] <= tolerance, drop=True)\n",
    "N = len(delta_FG_CO2['m'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compute Terms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Repetitive terms\n",
    "Sc = schmidt(forced['SST'])\n",
    "ff = co2_sol(forced['SST'], forced['SALT'])\n",
    "gamma_DIC = (3*means['ALK']*means['DIC'] - 2*means['DIC']**2) / \\\n",
    "            ((2*means['DIC']-means['ALK'])*(means['ALK']-means['DIC']))\n",
    "gamma_ALK = (-means['ALK']**2)/((2*means['DIC']-means['ALK'])*(means['ALK']-means['DIC']))\n",
    "\n",
    "\n",
    "\"\"\" \n",
    "The conversion term converts to the proper mol/m2/yr for every unit. So that we actually have each term\n",
    "summing to a mol/m2/yr perturbation from some climate event.\n",
    "\n",
    "Can run a np.where command on the original delta datasets to drop non-significant results.\n",
    "\"\"\"\n",
    "# SENSITIVITIES\n",
    "dFdpCO2 = means['RHO']*ff*a*np.sqrt(660/Sc)*means['U']**2\n",
    "dFdU = -2*a*np.sqrt(660/Sc)*means['U']*means['DCO2STAR']\n",
    "dFdSST = dFdpCO2*0.0423*means['pCO2']\n",
    "dFdSALT = dFdpCO2*(means['pCO2']/means['SALT'])\n",
    "dpCO2dDIC = (means['pCO2']*gamma_DIC)/means['DIC']\n",
    "dpCO2dALK = (means['pCO2']*gamma_ALK)/means['ALK']\n",
    "\n",
    "# FULL TERM\n",
    "F = delta_FG_CO2['m']\n",
    "U = dFdU*delta_U['m']*conversion\n",
    "SST = dFdSST*delta_SST['m']*conversion\n",
    "SALT = dFdSALT*delta_SALT['m']*conversion\n",
    "sDIC = dFdpCO2*dpCO2dDIC*(means['SALT']/S0)*delta_sDIC['m']*conversion\n",
    "sALK = dFdpCO2*dpCO2dALK*(means['SALT']/S0)*delta_sALK['m']*conversion\n",
    "\n",
    "# Freshwater\n",
    "term1 = dFdpCO2*dpCO2dDIC*(means['sDIC']/S0)*delta_SALT['m']*conversion\n",
    "term2 = dFdpCO2*dpCO2dALK*(means['sALK']/S0)*delta_SALT['m']*conversion\n",
    "FW = term1 + term2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "U: -0.01 + 0.0\n",
      "SST: -0.12 + 0.02\n",
      "SALT: -0.01 + 0.0\n",
      "sDIC: 0.11 + 0.04\n",
      "sALK: -0.01 + 0.02\n",
      "FW: -0.01 + 0.0\n"
     ]
    }
   ],
   "source": [
    "variables = ['U','SST','SALT','sDIC','sALK','FW']\n",
    "\n",
    "for i in variables:\n",
    "    newvar = eval(i)\n",
    "    ensmean = np.round(newvar.mean().values, 2)\n",
    "    spread = np.round(newvar.std().values, 2)\n",
    "    print(i + \": \" + str(ensmean) + \" + \" + str(spread))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Approx : -0.03 + 0.01\n",
      "True : -0.1 + 0.01\n"
     ]
    }
   ],
   "source": [
    "# Approx. sum\n",
    "taylor_sum = U + SST + SALT + sDIC + sALK\n",
    "ensmean = np.round(taylor_sum.mean().values, 2)\n",
    "spread = np.round(taylor_sum.std().values, 2)\n",
    "print(\"Approx : \" + str(ensmean) + \" + \" + str(spread))\n",
    "\n",
    "# Modeled\n",
    "ensmean = np.round(F.mean().values, 2)\n",
    "spread = np.round(F.std().values, 2)\n",
    "print(\"True : \" + str(ensmean) + \" + \" + str(spread))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sensitivities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-1.7912809155897704e-06\n",
      "1.553407664155826e-05\n",
      "1.1329994197344304e-05\n",
      "2.528269010672599e-06\n",
      "-2.1038629363655715e-06\n"
     ]
    }
   ],
   "source": [
    "variables = ['U','SST','SALT','DIC','ALK']\n",
    "for i in variables:\n",
    "    if i not in ['U','SST','SALT']:\n",
    "        temp = eval('dpCO2d' + i)\n",
    "        newvar = dFdpCO2 * temp\n",
    "    else:\n",
    "        newvar = eval('dFd' + i)\n",
    "    val = newvar.mean().values\n",
    "    print(str(val))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Regressions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.12195940676331439\n",
      "-0.24250648354388704\n",
      "-0.02038709660145013\n",
      "1.5407108338167197\n",
      "0.12688637454996382\n"
     ]
    }
   ],
   "source": [
    "variables = ['delta_U','delta_SST','delta_SALT','delta_sDIC','delta_sALK']\n",
    "for i in variables:\n",
    "    newvar = eval(i)\n",
    "    val = newvar['m'].mean().values\n",
    "    print(str(val))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Sensitivity vs Delta\n",
    "Quick vis to look at for different systems."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "d1 = {'U': [means['U'].values],\n",
    "     'DCO2STAR': [means['DCO2STAR'].values],\n",
    "     'Schmidt': [np.sqrt(660/Sc).values]}\n",
    "d2 = {'Sensitivity(e-3)': [dFdU.values*1000],\n",
    "      'Regression': [delta_U['m'].mean().values]}\n",
    "df1 = pd.DataFrame(data=d1)\n",
    "df2 = pd.DataFrame(data=d2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,1,'HumCS NINO3')"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAd8AAADOCAYAAAB//QjSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xt4TNf+x/F3LkJcg0riVEQRmkrcW/qjnI5cEBGS0HoUdalWKdrSntNqNCl1qFZDFXGJ9ke11dCUaatECOW4n0ZKL5ymEj+mLkXdMoz5/eGxT+ckNCQmOvm8nsfz2Huvtee792T2d9aavfZys9vtdkRERMRp3Ms6ABERkfJGyVdERMTJlHxFREScTMlXRETEyZR8RUREnEzJV0RExMmUfEVERJxMyfc2MJlMbNmyxWHdihUr6Nev3217zZ9++onRo0fTrl072rRpQ3R0NKmpqdhsNgCWL19O165dadWqFf/zP//D8OHDOXv2bJH7GjBgAKGhoRw5csRYt2XLFkwmk7H8+2NcsWIFTZs2ZcGCBQ776dSpE9u2bTOWDxw4wFNPPUWbNm1o1aoVAwYMYPfu3cb2kydP8uijj9KuXTvatm3LI488wq5du0p+ckRE7jBKvi7g0KFD9O3bl7p167Jq1Sp27dpFcnIyOTk5nDt3ju3btzNjxgzeeust9uzZw+eff063bt1uuM/KlSvz7rvvFjsGHx8f5s+ff92EfujQIfr160fTpk3JyMhg06ZNhIeHM3ToUPbs2QNAlSpVeP3119m6dSs7duzgiSeeYMSIEVy+fLn4J0NE5E9AybeMNG3alJ9//tlY/tvf/saMGTMA2LZtG506dWL+/Pk8+OCDdOzYkXXr1rFx40YiIyN54IEHmDt3rlF35syZtGrVir///e/4+voC0LBhQ958802qV6/O3r17admyJffddx9wNVH27t2bqlWrXje+AQMGsHr1aocYb6Rhw4a0atWKxYsXF7l91qxZtGzZkmeffRYfHx+qVq3KwIED6dmzJ9OnTwegYsWKNGzYEHd3d+x2O+7u7pw+fZrTp08XKwYRkT8LJd871PHjxykoKCArK4vRo0czYcIEPvvsM9LS0li6dCmzZ88mLy8PgK1btxIZGXndfbVo0YLNmzczc+ZMdu3ahdVq/cPX9/Pzo2/fvsyaNavYMY8ZM4bFixdz6tSpQtu2bNlC165dC63v1q0bu3fv5sKFC8a66OhomjdvzogRI+jTpw+1a9cudgwiIn8GnmUdgKsaOXIkHh4exvKlS5eMlmdxeHp6MmLECDw8POjevTuvvPIKAwcOpGrVqgQFBREUFMT3339PQEAAp06dok6dOtfdV9u2bZk1axbLli3j/fffx2az0bdvX1544QWHGP/bk08+SXh4OD/++GOxYg4ODqZDhw7Mnz+f8ePHO2z79ddfi4yxTp06XLlyhTNnzuDt7Q3AqlWrKCgoYO3atVy6dKlYry0i8meilu9tMnv2bHbu3Gn8mzhx4k3V9/HxMRJjpUqVABxagBUrVuTcuXNG2WPHjt1wf507d2bu3Lls376dd999l5UrV7J8+fIb1qlVqxaPPfYYM2fOLHbco0ePZtmyZYXiqVmzZpExHjt2DHd3d6pXr+6wvmLFivTo0YOUlBS+++67Yr++iMifgZJvGfH29nboav2j5HkjDz74IF999VWxyrq7u/Pggw/Svn37YrVohw4dyrZt28jJySnW/hs1akRERATz5s0rFOOXX35ZqPwXX3xBy5YtjVbvf7t8+bLRvS4i4iqUfMvIvffey+rVq7HZbGRlZbFjx45b3tfo0aPZs2cPU6dONZL4zz//zLhx4zhz5gzr1q3DbDZz+vRp7HY72dnZbN++nRYtWvzhvqtXr87gwYNZuHBhseMZOXIkaWlp/Pbbb8a6UaNGsWfPHmbMmMGpU6c4e/Ys//u//0t6ejrjxo0D4F//+hc7d+7EarVy8eJFUlJSOH78OM2bN7/JMyIicmdT8i0jL7/8MpmZmbRt25ZVq1YRFhZ2y/uqX78+H374IYcPH6ZHjx60adOGZ555hpCQEKpUqUKNGjX4+OOPiYiIoHXr1owfP56hQ4fSs2fPYu1/4MCBuLsX/08lICCAmJgYzp8/b6xr0KABH3zwAd999x0mk4mHHnqIr776igULFtCmTRsArFYrSUlJtGvXjk6dOpGVlUVKSgp+fn43d0JERO5wbna73V7WQYiIiJQnavmKiIg4Wakk38WLFxMVFUWPHj147rnnKCgoKI3disgtyMrKIjIykvDwcFJSUq5b7ssvv6Rp06bs3bvXWDdv3jzCw8OJjIxk06ZNzghXpFwqcfK1WCy8//77pKWlGTcQmc3m0ohNRG6SzWYjKSmJBQsWYDabWb16NQcOHChU7toNb7+/6e7AgQOYzWbMZjMLFiwgMTHReDa4iJSuUmn52mw2Ll68yOXLl7l48aLxiEMRca7s7GwCAwMJCAjAy8uLqKgoMjIyCpVLTk5m2LBhVKxY0ViXkZFBVFQUXl5eBAQEEBgYSHZ2tjPDFyk3SvyEKz8/P4YMGcLDDz9MxYoV6dChAx07dnQok5ub6zCmVUSKFhwcXKL6FosFf39/Y9nPz69QAt23bx9Hjx7l4YcfZtGiRQ51f98S9vPzw2KxFHoNfZ5FiudGn+cSJ9/Tp0+TkZFBRkYG1apVY8yYMaSnpxMTE2OUadCgQUlfRkSKoajBC25ubsb/r1y5wpQpU5gyZcpN171Gn2eRkitxt/OWLVuoV68etWrVokKFCkRERBhTxImIc/n7+3P06FFj2WKxOPwMdO7cOX744QcGDhyIyWTiX//6FyNGjGDv3r1/WFdESk+Jk+9f/vIXvvnmGy5cuIDdbmfr1q00atSoNGITkZsUGhpKbm4ueXl5WK1WzGYzJpPJ2F6tWjW2bdvG+vXrWb9+PS1btmTOnDmEhoZiMpkwm81YrVby8vLIzc3V08VEbpMSdzu3aNGCyMhIevfujaenJ8HBwTzyyCOlEZuI3CRPT08SEhIYNmwYNpuNuLg4goKCSE5OJiQkhC5duly3blBQEN26daN79+54eHiQkJBww1mvROTW6QlXIiIiTqYnXImIiDiZkq+IiIiTKfmKiIg4mZKviIiIkyn5ioiIOJmSr4iIiJMp+YqIiDiZkq+IiIiTKfmKiIg4mZKviIiIkyn5ioiIOJmSr4iIiJMp+Yq4mKysLCIjIwkPDyclJaXQ9mXLlhEdHU1MTAz9+vXjwIEDAOTn59O8eXNiYmKIiYkhISHB2aGLlBua1UjEhdhsNiIjI0lNTcXPz4/4+HjeeustGjdubJQ5e/YsVatWBSAjI4MPPviAhQsXkp+fz1NPPcXq1avLKnyRckMtXxEXkp2dTWBgIAEBAXh5eREVFUVGRoZDmWuJF+DChQu4ubk5O0yRcs+zrAMQkdJjsVjw9/c3lv38/MjOzi5UbunSpaSmpnLp0iXee+89Y31+fj69evWiatWqjB07lrZt2xaqm5uby4ULF27PAYi4kODg4OtuU/IVcSFF/YpUVMu2f//+9O/fn1WrVjFnzhymTp2Kr68vmZmZ1KxZk5ycHEaOHInZbHZoKQM0aNDgdoUvUm6o21nEhfj7+3P06FFj2WKx4Ovre93yUVFRrFu3DgAvLy9q1qwJQEhICPXr1+enn366vQGLlFOlknzPnDnD6NGj6dq1K926dWPPnj2lsVsRuUmhoaHk5uaSl5eH1WrFbDZjMpkcyuTm5hr/37BhA4GBgQCcPHkSm80GQF5eHrm5uQQEBDgtdpHypFS6nSdPnsxDDz3EzJkzsVqtXLx4sTR2KyI3ydPTk4SEBIYNG4bNZiMuLo6goCCSk5MJCQmhS5cuLFmyhK1bt+Lp6Un16tWZOnUqADt27GDmzJl4eHjg4eFBYmIiPj4+ZXxEIq6pxEONzp49S8+ePcnIyNBdkyIiIsVQ4pZvXl4etWrV4u9//zvfffcdzZo14+WXX6Zy5cpGGd0dWXwff/xxWYdwW/Tt27esQ/hTuNHdkSLiOkqcfC9fvsy+fft45ZVXaNGiBZMmTSIlJYWxY8caZXR3pCipiIj8R4lvuPL398ff358WLVoA0LVrV/bt21fiwERERFxViZNvnTp18Pf359///jcAW7dupVGjRiUOTERExFWVyt3Or7zyCuPGjePSpUsEBAQwZcqU0titiIiISyqV5BscHMyKFStKY1ciIiIuT0+4EhERcTIlXxERESdT8hUREXEyJV8REREnU/IVERFxMiVfEReTlZVFZGQk4eHhpKSkFNq+bNkyoqOjiYmJoV+/fhw4cMDYNm/ePMLDw4mMjGTTpk3ODFukXCmVoUYicmew2WwkJSWRmpqKn58f8fHxmEwmGjdubJSJjo6mX79+AGRkZDBlyhQWLlzIgQMHMJvNmM1mLBYLgwcPZs2aNXh4eJTV4Yi4LLV8RVxIdnY2gYGBBAQE4OXlRVRUFBkZGQ5lqlatavz/woULxmxkGRkZREVF4eXlRUBAAIGBgWRnZzs1fpHyQi1fERdisVjw9/c3lv38/IpMoEuXLiU1NZVLly7x3nvvGXWvPaP9Wl2LxVKormYpEymeG00oo+Qr4kKKmp67qHm2+/fvT//+/Vm1ahVz5sxh6tSpxa6rWcpESk7dziIuxN/fn6NHjxrLFosFX1/f65aPiopi3bp1t1RXRG6dkq+ICwkNDSU3N5e8vDysVitmsxmTyeRQJjc31/j/hg0bCAwMBMBkMmE2m7FareTl5ZGbm0vz5s2dGb5IuaFuZxEX4unpSUJCAsOGDcNmsxEXF0dQUBDJycmEhITQpUsXlixZwtatW/H09KR69epMnToVgKCgILp160b37t3x8PAgISFBdzqL3CZu9qJ+6JEyk5iYWNYh3BYTJ04s6xBERO4Y6nYWERFxMiVfERERJ1PyFRERcTIlXxERESdT8hUREXGyUku+NpuNXr168eSTT5bWLkVERFxSqSXf999/n0aNGpXW7kRERFxWqSTfo0ePsmHDBuLj40tjdyIiIi6tVJ5w9frrrzN+/HjOnTtX5HbNgiL79+8v6xD+FG40C4qIuI4SJ9/MzExq1apFSEgI27ZtK7KMZkERJRURkf8ocfLdvXs369evJysri4KCAs6ePcu4ceOYPn16acQnIiLickqcfJ9//nmef/55ALZt28aiRYuUeEXKUFZWFpMnT+bKlSv06dOH4cOHO2xPTU1l+fLleHh4UKtWLV5//XXuvvtu4GoPRZMmTQCoW7cuc+fOdXr8IuWBZjUScSE2m42kpCRSU1Px8/MjPj4ek8lE48aNjTLBwcGkpaXh7e3NBx98wBtvvMHbb78NQKVKlUhPTy+r8EXKjVJ9yEa7du2YN29eae5SRG5CdnY2gYGBBAQE4OXlRVRUFBkZGQ5l2rdvj7e3NwAtW7bk6NGjZRGqSLmmlq+IC7FYLPj7+xvLfn5+ZGdnX7f8J598QqdOnYzlgoICYmNj8fT0ZPjw4YSFhRWqo9ELIsVzoxtNlXxFXEhR03O7ubkVWTY9PZ2cnByWLFlirMvMzMTPz4+8vDwGDRpEkyZNqF+/vkM9jV4QKTk921nEhfj7+zt0I1ssFnx9fQuV27JlC3PnzmXOnDl4eXkZ6/38/AAICAjggQceYN++fbc/aJFySMlXxIWEhoaSm5tLXl4eVqsVs9mMyWRyKLNv3z4SEhKYM2cOtWvXNtafPn0aq9UKwMmTJ9m9e7fDjVoiUnrU7SziQjw9PUlISGDYsGHYbDbi4uIICgoiOTmZkJAQunTpwrRp0zh//jxjxowB/jOk6ODBg0ycOBE3NzfsdjtPPPGEkq/IbeJmL+pHIikziYmJZR3CbTFx4sSyDkFE5I6hbmcREREnU/IVERFxMiVfERERJ1PyFRERcTIlXxERESdT8hUREXEyJV8REREnU/IVERFxMiVfERERJ1PyFRERcTIlXxERESdT8hUREXGyEiffI0eOMGDAALp160ZUVBTvvfdeacQlIrcoKyuLyMhIwsPDSUlJKbQ9NTWV7t27Ex0dzaBBgzh8+LCxbeXKlURERBAREcHKlSudGbZIuVLiKQU9PDz429/+RrNmzTh79ixxcXF06NBBU5GJlAGbzUZSUhKpqan4+fkRHx+PyWRy+DwGBweTlpaGt7c3H3zwAW+88QZvv/02p06d4p133iEtLQ03NzdiY2MxmUzUqFGjDI9IxDWVuOXr6+tLs2bNAKhatSoNGzbEYrGUODARuXnZ2dkEBgYSEBCAl5cXUVFRZGRkOJRp37493t7eALRs2ZKjR48CsHnzZjp06ICPjw81atSgQ4cObNq0yenHIFIelLjl+3v5+fns37+fFi1aOKzPzc3lwoULpflS8iezf//+sg7hTyE4OLhE9S0WC/7+/sayn58f2dnZ1y3/ySef0KlTp+vWLeqLtD7PIsVzo89zqSXfc+fOMXr0aF566SWqVq3qsK1Bgwal9TLyJ1XSpCLFY7fbC61zc3Mrsmx6ejo5OTksWbLkpurq8yxScqVyt/OlS5cYPXo00dHRRERElMYuReQW+Pv7G93IcLU16+vrW6jcli1bmDt3LnPmzMHLy+um6opIyZU4+drtdl5++WUaNmzI4MGDSyMmEblFoaGh5ObmkpeXh9VqxWw2YzKZHMrs27ePhIQE5syZQ+3atY31HTt2ZPPmzZw+fZrTp0+zefNmOnbs6OxDECkXStztvGvXLtLT02nSpAkxMTEAPPfcc3Tu3LnEwYnIzfH09CQhIYFhw4Zhs9mIi4sjKCiI5ORkQkJC6NKlC9OmTeP8+fOMGTMGgLp16zJ37lx8fHx4+umniY+PB2DkyJH4+PiU5eGIuCw3e1E/9EiZSUxMLOsQbouJEyeWdQgiIncMPeFKRETEyZR8RUREnEzJV0RExMmUfEVERJxMyVdERMTJlHxFREScTMlXRETEyUp1YoVboXGtIiJS3qjlKyIi4mRKviIiIk6m5CsiIuJkSr4iLiYrK4vIyEjCw8NJSUkptH3Hjh307t2b++67jy+//NJhW3BwMDExMcTExPDUU085K2SRcqfMb7gSkdJjs9lISkoiNTUVPz8/4uPjMZlMNG7c2ChTt25dpkyZwqJFiwrVr1SpEunp6c4MWaRcUvIVcSHZ2dkEBgYSEBAAQFRUFBkZGQ7Jt169egC4u6vjS6SsKPmKuBCLxYK/v7+x7OfnR3Z2drHrFxQUEBsbi6enJ8OHDycsLKxQmdzcXC5cuFAq8Yq4suDg4OtuU/IVcSFFTc/t5uZW7PqZmZn4+fmRl5fHoEGDaNKkCfXr13co06BBg5KGKVLuqd9JxIX4+/tz9OhRY9liseDr61vs+n5+fgAEBATwwAMPsG/fvlKPUUSUfEVcSmhoKLm5ueTl5WG1WjGbzZhMpmLVPX36NFarFYCTJ0+ye/duh9+KRaT0lEq3c1ZWFpMnT+bKlSv06dOH4cOHl8ZuReQmeXp6kpCQwLBhw7DZbMTFxREUFERycjIhISF06dKF7OxsRo0axZkzZ8jMzGTWrFmYzWYOHjzIxIkTcXNzw26388QTTyj5itwmJU6+xRnaICLO07lzZzp37uywbsyYMcb/mzdvTlZWVqF6rVu3ZtWqVbc9PhEphW7n3w9t8PLyMoY2iIiISNFK3PItztCGGw1N6Nu3b0lDuCPt37//lurpfPzHxx9/fBsiuTNc732+0dAEEXEdJU6+xRnaoKEJIo6UZEXKtxJ3O5d0aIOIiEh5U+LkW5KhDSIiIuVRibudrze0QURERIpWKuN8ixraICIiIkXTE65EREScTMlXRETEyZR8RUREnEzJV0RExMmUfEVERJxMyVfExWRlZREZGUl4eDgpKSmFtu/YsYPevXtz33338eWXXzpsW7lyJREREURERLBy5UpnhSxS7pTKUCMRuTMUZ5axunXrMmXKFBYtWuRQ99SpU7zzzjukpaXh5uZGbGwsJpOJGjVqOPswRFyeWr4iLqQ4s4zVq1ePe++9F3d3x4//5s2b6dChAz4+PtSoUYMOHTqwadMmZ4YvUm6o5SviQoozy9jN1LVYLIXK3WiWMhH5jxtNoKLkK+JCijPLWEnrapYykZJT8hVxISWZZczf35/t27c71H3ggQduOZbExMRbritXTZw4saxDkNtEv/mKuJCSzDLWsWNHNm/ezOnTpzl9+jSbN2+mY8eOtzlikfJJLV8RF3K9WcaSk5MJCQmhS5cuZGdnM2rUKM6cOUNmZiazZs3CbDbj4+PD008/TXx8PAAjR47Ex8enjI9IxDW52Yv6oUfkDuDK3ZbloTvRld8/ZykPfyfllbqdRUREnEzJV0RExMmUfEVERJxMyVdERMTJSnS389SpU8nMzKRChQrUr1+fKVOmUL169dKKTURExCWVqOXboUMHVq9ezapVq2jQoAHz5s0rrbhERERcVomSb8eOHfH0vNp4btmypcOTdURERKRopfaQjbS0NLp161bkNj2IXcTR/v37i1x/owexi4jr+MPk+/jjj3P8+PFC68eOHUtYWBgAc+bMwcPDg549exa5Dz2IXcSRkqxI+faHyXfx4sU33L5y5Uo2bNjA4sWLiz17ioiISHlWom7nrKws5s+fz5IlS/D29i6tmEQAPVpPRFxXiW64eu211zh37hyDBw8mJiaGhISE0opLRG5RVlYWkZGRhIeHk5KSUmi71Wpl7NixhIeH06dPH/Lz8wHIz8+nefPmxMTE6PMscpuVqOW7du3a0opDREqBzWYjKSmJ1NRU/Pz8iI+Px2Qy0bhxY6PM8uXLqV69OmvXrsVsNjN9+nTefvttAOrXr096enpZhS9SbugJVyIuJDs7m8DAQAICAvDy8iIqKoqMjAyHMuvXr6d3794AREZGsnXrVjS5mYhzaT5fERdisVjw9/c3lv38/MjOzi5Upm7dusDV+X+rVavGr7/+Clzteu7VqxdVq1Zl7NixtG3bttBraOig81xvSJr8OdxoVIOSr4gLKaoF+9+jEK5XxtfXl8zMTGrWrElOTg4jR47EbDZTtWpVh7IaOug8GpLmutTtLOJC/P39HZ40Z7FY8PX1LVTmyJEjAFy+fJnffvsNHx8fvLy8qFmzJgAhISHUr1+fn376yXnBi5QjSr4iLiQ0NJTc3Fzy8vKwWq2YzWZMJpNDGZPJxMqVKwFYs2YN7du3x83NjZMnT2Kz2QDIy8sjNzeXgIAApx+DSHmgbmcRF+Lp6UlCQgLDhg3DZrMRFxdHUFAQycnJhISE0KVLF+Lj4xk/fjzh4eHUqFGDGTNmALBjxw5mzpyJh4cHHh4eJCYm4uPjU8ZHJOKa3Oy6zVFEboPExMSyDuFPTw+acV3qdhYREXEyJV8REREnU/IVERFxMiVfERERJ1PyFRERcTIlXxERESfTOF8RkXJEQ8BKrjSGgKnlKyIi4mRKviIiIk6m5CsiIuJkSr4iIiJOVirJd+HChTRt2pSTJ0+Wxu5EpASysrKIjIwkPDyclJSUQtutVitjx44lPDycPn36kJ+fb2ybN28e4eHhREZGsmnTJmeGLVKulDj5HjlyhC1btvCXv/ylNOIRkRKw2WwkJSWxYMECzGYzq1ev5sCBAw5lli9fTvXq1Vm7di2PP/4406dPB+DAgQOYzWbMZjMLFiwgMTHRmGJQREpXiZPvlClTGD9+PG5ubqURj4iUQHZ2NoGBgQQEBODl5UVUVBQZGRkOZdavX0/v3r0BiIyMZOvWrdjtdjIyMoiKisLLy4uAgAACAwPJzs4ui8MQcXklGuebkZGBr68v9957b2nFIyIlYLFY8Pf3N5b9/PwKJVCLxULdunWBq/P/VqtWjV9//RWLxUKLFi0c6loslluORdPh3Zn0vtwZ/jD5Pv744xw/frzQ+rFjxzJv3jwWLVp0WwITkZtX1PTc/90rdb0yxakrIqXjD5Pv4sWLi1z//fffk5+fT0xMDABHjx4lNjaW5cuXU6dOnVINUkSKx9/fn6NHjxrLFosFX1/fQmWOHDmCv78/ly9f5rfffsPHx6dYdUWkdNzyb75NmzZl69atrF+/nvXr1+Pv78+KFSuUeEXKUGhoKLm5ueTl5WG1WjGbzZhMJocyJpOJlStXArBmzRrat2+Pm5sbJpMJs9mM1WolLy+P3NxcmjdvXhaHIeLy9GxnERfi6elJQkICw4YNw2azERcXR1BQEMnJyYSEhNClSxfi4+MZP3484eHh1KhRgxkzZgAQFBREt27d6N69Ox4eHiQkJODh4VHGRyTimtzsRf3QcweZM2cOq1evxt3dHXd3d5KSkhxuCrlmxYoV5OTkkJCQcEuvs3fvXtLT05kwYUKhbSaTiU8++QRPT09WrVpF//79b+k17gT5+fk89dRTrF692lg3a9YsKleuzNChQ8swsqIFBwfTpEkTLl++jIeHB71792bQoEG4u1/ttMnOzmbq1KkcP34cNzc32rRpw4QJE/D29mbdunUkJydz6dIlKlSowJgxYwgLCwNg6tSpZGZmUqFCBerXr8+UKVOoXr06Fy5cYMKECfzwww/Y7XaqVavG9OnTefrppwE4fvw47u7u1KpVC7g6bMfLy4u1a9cyatQoPv/8cxo1agRcPdfdu3fnnnvu4dKlS4SEhDB58mQqVKhQBmfyz+Pae26z2ahXrx7Tpk2jevXqZR2W4YknnuDNN9+8o2IqieJeY2/Fo48+yocffkh+fj579uwhOjoauPH19pply5bh7e1Nr169WLFiBR06dMDPz++Gr2e32xk0aBDvvvsuVatWLVaM2dnZvPLKK0b9Z555hvDwcKxWK4MHD+a9997D0/M2tFPtd7Ddu3fb+/btay8oKLDb7Xb7iRMn7EePHi2ybFpamj0xMfG2xPHwww/bT5w4Yc/Ly7NHRUXdltdwlqKOYebMmfYFCxaUUUQ31rJlS+P/x48ftw8aNMienJxst9vt9mPHjtn/+te/2nfv3m232+32K1eu2L/44gv7sWPH7Pv377eHhYXZDx06ZLfb7fZDhw7Zw8LC7Pv377fb7Xb7pk2b7JcuXbLb7Xb7tGnT7NOmTbPb7Xb73Llz7a+MeL3CAAANEElEQVS//rrxmgcPHjT+/uz265+r0aNH2/v162efOXOmse735/ry5cv2AQMG2NPT00t+Ulzc79/zF154wf7uu++Wyn6vvd/yHzdzjS2Jf/7zn/bhw4ffcv3HHnvMnp2d/YflMjMz7ZMnT76pfZ8/f97427BYLPb27dsby7Nmzbptn9k7+vGSx44do2bNmnh5eQFQq1YtY+jEo48+Ss+ePYmPj+fs2bMA/PLLLwwdOpSIiAimTZtm7KdVq1a88cYbxMbG8vjjj5Odnc2AAQPo0qWLMQZy27ZtPPnkkwD8+uuvDBkyhF69epGQkGDcBfrmm29y6NAhYmJimDp1qjNPhQC1a9fmtddeY+nSpdjtdpYuXUqvXr1o1aoVcPXO3K5du3LXXXexcOFCnnzySQICAgAICAhg+PDhLFy4EICOHTsa32Zbtmxp3Gh07Ngxh2/XDRs2NP7+rufcuXPs3r2byZMnYzabiyzj4eFB8+bNSzR0pzxq2bKlwzlbsGABcXFxREdHM3PmTGP97Nmz6dq1K4MHD+a5554z3ucBAwbw1ltv8dhjj/H+++9z8uRJnnnmGeLi4oiLi2PXrl0AbN++nZiYGGJiYujVqxdnz57ll19+oX///sTExNCjRw927twJXO0Ju/Y0v9TUVHr06EGPHj2Mm1Pz8/Pp1q0bEyZMICoqiiFDhnDx4kVnnK6bdr1rbE5ODo899hixsbEMHTqUX375Bbh6Pt944w3i4+OJjIw0zsmPP/5IfHw8MTExREdHk5ubC2B8Nt9880127txJTEwMixcvNq63V65cwWQycebMGSOm8PBwjh8/zqxZs1i4cCFffvklOTk5jBs3jpiYGDZs2MDIkSON8l9//TWjRo0CYNWqVXTp0sXYlp6ebsSVkJBQ5ENjvL29jWtBQUGBwx3+YWFhrFq1qsTnuSh3dPLt0KEDR44cITIykldffZXt27djtVp59tlneemll/jss89YvHgxlSpVAmD//v28/fbbrFq1ii+++IIjR44AcP78eR544AFWrFhBlSpVePvtt1m0aBGzZ892+ABfM3v2bFq3bs2nn36KyWTi//7v/wB4/vnnqV+/Punp6bz44ovOOxFiCAgI4MqVK5w4cYIff/yRZs2aFVnuwIEDhISEOKwLDQ0t9LQngLS0NDp16gRAXFwc8+fP55FHHmHGjBnGReRG1q1bx0MPPcQ999yDj48P3377baEyBQUFfPPNNzz00EPFOEqBq0/r2rp1q3HD2ObNm/n555/55JNPSE9P59tvv2XHjh3s3buXr776ik8//ZRZs2aRk5PjsJ8zZ86wZMkShgwZwuTJkxk0aBBpaWnMmjXL6PZctGgRCQkJpKens3TpUipVqsTq1avp2LEj6enppKenF3qeQU5ODitWrODjjz/mo48+Yvny5ezbtw+An3/+mf79+2M2m6lWrRpr1qxxwhm7eUVdYy9dusSkSZOYOXMmK1asIC4uzrgvAK6+L5988gkvvfQS77zzDgAffvghAwcOJD09nbS0NIex5nD12tm2bVvS09N5/PHHjfXu7u6YTCbWrl0LwDfffMPdd9/NXXfdZZTp2rUrISEhTJ8+nfT0dDp37szBgweNL0ArVqwgNjYWgN27dxvXhIMHD/LFF1+wbNky0tPTcXd3v24i/eabb4iKiqJnz54kJiYayTgoKIi9e/eW5BRf1x19w1WVKlVYsWIFO3fuZNu2bTz77LM89dRT1KlTx7gL8/f9+g8++CDVqlUDoFGjRhw+fJi6detSoUIF4+LapEkTvLy8qFChAk2aNOHw4cOFXnfHjh3GH9Vf//pXatSocbsP1WmuN27zzzSe016M2xSKKmO32wsd55w5c/Dw8KBnz57A1d8b161bx9dff82WLVuIj4/no48+Mn7HLYrZbGbQoEEAdO/endWrVxsXgGs9JT///DORkZF6IE0xXLx4kZiYGA4fPkyzZs3o0KEDcLWF8/XXX9OrVy/g6pfq3Nxczp07R5cuXYwv4Q8//LDD/rp37278f8uWLQ5fwM6ePcvZs2dp3bo1//jHP4iOjiYiIoIqVaoQGhrKSy+9xOXLlwkLCyM4ONhhv7t27SIsLIzKlSsDV1tsO3fuxGQyUa9ePaN8s2bNirzO3AmKusaOGDGCH374gcGDBwNw5coVh1Es4eHhgONxtWzZkrlz53L06FEiIiJo0KBBsWPo3r07s2fPJi4uDrPZ7PB+FcXNzY2YmBg+++wzYmNj2bNnj9ETeerUKSMnbN26lZycHOLj44Grf1e1a9cucp8tWrTAbDZz8OBBXnzxRTp16kTFihXx8PCgQoUKnD17tti/IRfXHZ184Wp3Xbt27WjXrh1NmjRh6dKl100Uv+8e9PDwMLoYKlSoYNRxd3c3yrm7u5e7Z9f6+Phw+vRph3WnT5+mXr16ZRTRzcnLy8PDw4PatWvTuHFjvv32W+Mmqt9r3LgxOTk5Dsnu22+/dUiiK1euZMOGDSxevNjhb6pKlSpEREQQERGBu7s7GzduvG7y/fXXX/nnP//Jjz/+iJubGzabDTc3N1544QUAo6fkl19+YcCAAWRkZDh0i0lhlSpVIj09nd9++40nn3ySpUuXMnDgQOx2O8OHD+fRRx91KH+9ZxFc4+3tbfz/ypUrfPTRR0aivmb48OF07tyZjRs30rdvX1JTU7n//vtZsmQJGzdu5IUXXmDo0KFG4ocbfwn872tRQUFBcQ69TBR1jQ0KCuKjjz4qsnxR18/o6GhatGjBhg0bGDp0KJMmTeLBBx8s1uu3atWKQ4cOcfLkSdatW8eIESP+sE5sbCwjRozAy8uLrl27Gi1VT09Prly5gru7O3a7nd69e/P888871F27dq3RuJo0aRKhoaHGtkaNGuHt7c0PP/xgrLdarVSsWLFYx3Iz7uhu53//+98O3X779++nUaNG/PLLL8Yj886ePcvly5dL9XXvv/9+o3ti48aNRrKqUqUK586dK9XXcrYqVapQp04dtm7dClz9prhp0ybatGlTxpH9sZMnTzJx4kT69++Pm5sbjz32GJ9++inffPONUSY9PZ1jx44xdOhQUlJSjBl78vPzmTdvHkOGDAGuzvwzf/585syZ43Bx3rVrl/F+W61WDhw4cMNJQ9asWUOvXr3IzMxk/fr1bNy4kXr16hm/JV7j6+vLuHHjipxlSIpWrVo1JkyYwKJFi7h06RIdO3YkLS3N+AxaLBZOnDhB69atyczMpKCggHPnzrFhw4br7rNjx44sWbLEWN6/fz9wtYeiadOmDB8+nJCQEH766ScOHz5M7dq16du3L3FxcYV+Trj//vtZt24dFy5c4Pz586xbt462bduW/om4ja53jT158iR79uwB4NKlS/z444833E9eXh4BAQEMHDgQk8nE999/77D9RtdONzc3wsLCmDJlCo0aNaJmzZqFyvx3fT8/P3x9fZkzZ47R5Qxwzz33kJeXB1ztCV2zZg0nTpwArl7rDh8+THh4uPFTQmhoKHl5eUYOOXz4MD/99BN33303cPXLda1atW7LCIU7uuV7/vx5Jk2axJkzZ/Dw8CAwMJCkpCRiY2OZNGkSFy9epFKlSqSmppbq644cOZLnn3+e3r17c//99xsX35o1a9K6dWt69OjBQw899Kf93XfatGkkJibyj3/8A7h6vPXr1y/jqIp2rQvy2lCjmJgYozvsrrvu4q233mLq1KmcOHECd3d37r//fiIiIqhTpw7jxo1jxIgRxlCj8ePHG12Br732mjGUAK52OyUlJZGXl8err74KXG0lde7cmcjIyOvGZzabeeKJJxzWRUREsGrVqkLrw8LCmDVrFjt37vzTXaTLyn333ce9996L2WymV69eHDx40Gj5Vq5cmTfeeIPmzZtjMpno2bMnd999NyEhIcbPT//t5ZdfJikpiejoaGw2G23btiUpKYn33nuPbdu24e7uTuPGjenUqRNms5mFCxfi6elJ5cqVC91k2axZM2JjY+nTpw8A8fHx3HfffQ5TNN7prneNfeSRR5g0aRK//fYbNpuNQYMGERQUdN39fP7553z22Wd4enpy1113OdwQBVcfynTt553Y2NhCXfjdu3cnPj7euCb9t969ezNx4kQqVapk9FxER0dz8uRJGjdubJTr3Lkz27dvJzAwkMaNGzN27FiGDBnClStXqFChAgkJCUZivWbXrl3Mnz8fT09P3N3defXVV42hhNu2baNz5843dU6L644f5ysi8kfOnTtHlSpVuHDhAv379+e111677s144hqSkpIIDg42vvzA1REvL774Yqk1yEaNGsVzzz1Hw4YNS2V/v3dHt3xFRIojISGBAwcOUFBQQO/evZV4XVxsbCze3t787W9/c1jv6+tLnz59SuUGKavVSlhY2G1JvKCWr4iIiNPd0TdciYiIuCIlXxERESdT8hUREXEyJV8REREnU/IVERFxsv8H6mh4vQL5b1UAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x216 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.set_style('whitegrid')\n",
    "order = list(['Schmidt','U','DCO2STAR'])\n",
    "f, ax = plt.subplots(ncols=2, figsize=(8,3))\n",
    "bars1 = sns.barplot(data=df1, ax=ax[0], order=order, color='gray')\n",
    "bars2 = sns.barplot(data=df2, ax=ax[1], color='gray')\n",
    "\n",
    "ax[0].set_ylim([-4,8])\n",
    "ax[1].set_ylim([0,0.4])\n",
    "sns.despine(left=True,bottom=True)\n",
    "\n",
    "ax[0].set_title(e + \" \" + x)\n",
    "\n",
    "#et.vis.savefig('cancs_nao_U', extension='.png', dpi=300)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def setup_dataframe(v, da):\n",
    "    \"\"\"\n",
    "    Creates a DataFrame with term values given a variable name (v) and the\n",
    "    DataArray containing the expansion results (da).\n",
    "    \"\"\"\n",
    "    da.name = v\n",
    "    df = da.to_dataframe()\n",
    "    del df['z_t']\n",
    "    return df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "variables = ['U','SST','SALT','sDIC','sALK','FW']\n",
    "df = pd.DataFrame()\n",
    "for i, name in enumerate(variables):\n",
    "    temp_df = setup_dataframe(name, eval(name))\n",
    "    df = pd.concat([df, temp_df], axis=1)\n",
    "df['Total'] = df.sum(axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# For putting the variance explained.\n",
    "avg_r = delta_FG_CO2['r'].median()\n",
    "avg_r = avg_r.values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set_style('whitegrid')\n",
    "if F.mean() > 0:\n",
    "    clist = ['r', '#a9a9a9', '#a9a9a9', '#a9a9a9', '#a9a9a9', '#a9a9a9', '#a9a9a9']\n",
    "else:\n",
    "    clist = ['b', '#a9a9a9', '#a9a9a9', '#a9a9a9', '#a9a9a9', '#a9a9a9', '#a9a9a9']\n",
    "order = list(['Total', 'U', 'SST', 'SALT', 'sDIC', 'sALK', 'FW'])\n",
    "f, ax = plt.subplots(figsize=(8,4))\n",
    "bars = sns.barplot(data=df, ax=ax, ci='sd', capsize=0.1, errwidth=1,\n",
    "                   palette=clist, order=order)\n",
    "sns.despine(left=True, bottom=True)\n",
    "ax.set_title(e + ' 1$\\sigma$' + ' ' + x + ' Response (N = ' + str(N) + ')', fontsize=14)\n",
    "ax.set_ylabel('Contribution to F$_{\\mathrm{CO}_{2}}$ Change' + '\\n' + '[mol/m2/yr]',\n",
    "              fontsize=14)\n",
    "\n",
    "# Add true regression\n",
    "ax.errorbar(0, F.mean(), yerr=F.std(), fmt='o', capsize=3, color='k', markersize=6)\n",
    "\n",
    "# Add variance explained.\n",
    "plt.text(5.5, 0.12, 'R = ' + str(avg_r.round(2)), fontsize=12, fontweight='bold',\n",
    "        bbox={'edgecolor':'w', 'facecolor':'w', 'alpha':0.2})\n",
    "\n",
    "# Fontsizes\n",
    "ax.tick_params(labelsize=13)\n",
    "\n",
    "# Save figure.\n",
    "#et.vis.savefig('CanCS_nao', extension='.eps', dpi=300)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Signs of Sensitivity and Regression Terms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mean(x):\n",
    "    return x['m'].mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def color_negative_red(val):\n",
    "    \"\"\"\n",
    "    Takes a scalar and returns a string with\n",
    "    the css property `'color: red'` for negative\n",
    "    strings, black otherwise.\n",
    "    \"\"\"\n",
    "    color = 'red' if val == '-' else 'green'\n",
    "    return 'color: %s' % color"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set up list of terms\n",
    "sens_list = np.asarray([dFdU, dFdSST, dFdSALT, (dFdpCO2*dpCO2dDIC), (dFdpCO2*dpCO2dALK)])\n",
    "reg_list = np.asarray([mean(delta_U), mean(delta_SST), mean(delta_SALT), mean(delta_sDIC),\n",
    "                       mean(delta_sALK)])\n",
    "term_list = np.asarray([U.mean(),SST.mean(),SALT.mean(),sDIC.mean(),sALK.mean()])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Replace with + and -\n",
    "sens = []\n",
    "reg = []\n",
    "term = []\n",
    "for s, r, t in zip(sens_list, reg_list, term_list):\n",
    "    if s < 0:\n",
    "        sens.append('-')\n",
    "    else:\n",
    "        sens.append('+')\n",
    "    if r < 0:\n",
    "        reg.append('-')\n",
    "    else:\n",
    "        reg.append('+')\n",
    "    if t < 0:\n",
    "        term.append('-')\n",
    "    else:\n",
    "        term.append('+')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.DataFrame({'Sensitivity': sens, '$\\Delta$': reg, 'Term': term},\n",
    "                  index=['U','SST','SALT','sDIC','sALK'])\n",
    "df = df[['Sensitivity','$\\Delta$','Term']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s = df.style.applymap(color_negative_red)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "means['DCO2STAR']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Degrees or Sigma?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds = xr.open_dataset('/glade/p/work/rbrady/cesmLE_CVDP/processed/cvdp_detrended_BGC.nc')\n",
    "#ds = xr.open_mfdataset('/glade/p/work/rbrady/EBUS_BGC_Variability/NPGO/*.nc', concat_dim='ensemble')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds['nao']"
   ]
  }
 ],
 "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}