Part_2_Calibration.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Recalibration of SMP density coefficients\n",
    "*Josh King, Environment and Climate Change Canada, 2020*\n",
    "\n",
    "The matched snow pit and SMP measurements from Part 1 are used to recalibrate the bilinear regression model of [Proksch et al. (2015)](https://agupubs.onlinelibrary.wiley.com/doi/full/10.1002/2014JF003266). We apply a K-Folds approach to minimize sampling bias while evaluating the OLS regression. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import community packages\n",
    "import string\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "plt.rcParams[\"font.family\"] = \"Times New Roman\"\n",
    "\n",
    "import matplotlib.cm as cm\n",
    "import pickle\n",
    "\n",
    "from scipy import stats\n",
    "from statsmodels.formula.api import ols\n",
    "from sklearn.model_selection import KFold, StratifiedShuffleSplit\n",
    "\n",
    "# Seed to replicate the paper result exactly\n",
    "RANDOM_SEED = 2019\n",
    "\n",
    "# Load comparison result from Part 1\n",
    "result = pd.read_pickle('./output/smp_pit_filtered.pkl')\n",
    "result['force_log'] = np.log(result['force_median'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## K-Folds OLS method\n",
    "\n",
    "OLS regression with 10 folds to minimize sampling bias.\n",
    "Model coefficients and skill are evaluated as the mean of all folds. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "k_fold = KFold(n_splits = 10, shuffle = True, random_state = RANDOM_SEED)\n",
    "rmse = []; error = []; r = []; params = None\n",
    "\n",
    "# Split the dataset into 10 roughly equal groups, \n",
    "# train on all but one test group\n",
    "for train_idx, test_idx in k_fold.split(result):\n",
    "    train = result.iloc[train_idx]\n",
    "    test = result.iloc[test_idx]\n",
    "    \n",
    "    model_rho = ols(\"RHO ~ force_log + force_log * l\", train).fit()\n",
    "    predict_rho = model_rho.predict(exog=dict(force_log=test['force_log'], l=test['l']))\n",
    "    rmse = np.append(rmse, np.sqrt(np.mean((predict_rho - test['RHO'])**2)))\n",
    "    r = np.append(r,np.corrcoef(predict_rho, test['RHO'])[1][0])\n",
    "    error = np.append(error, predict_rho - test['RHO'])\n",
    "    \n",
    "    if params is None:\n",
    "        params = model_rho.params.values\n",
    "    else:\n",
    "        params =  np.vstack((params, model_rho.params.values))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# K20a evaluation stats (from k-folds method)\n",
    "# Metrics presented as mean of all fold permutations\n",
    "\n",
    "k20a_rmse = rmse.mean()\n",
    "k20a_bias = error.mean()\n",
    "k20a_r2 = r.mean()**2\n",
    "\n",
    "print('K20a recalibration evaluation')\n",
    "print('N: %i' % len(result))\n",
    "print('RMSE: %0.1f' % k20a_rmse)\n",
    "print('RMSE percent: %0.2f' % np.round(k20a_rmse/result.RHO.mean(),2))\n",
    "print('bias: %0.1f' % k20a_bias)\n",
    "print('r^2: %0.2f' % k20a_r2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model_k20a_coeff = [np.round(params[:,0].mean(),2), np.round(params[:,1].mean(),2),\n",
    "              np.round(params[:,3].mean(),2), np.round(params[:,2].mean(),2)]\n",
    "var_coeffs = [np.round(params[:,0].std(),2), np.round(params[:,1].std(),2),\n",
    "              np.round(params[:,3].std(),2), np.round(params[:,2].std(),2)]\n",
    "\n",
    "# Unbiased coeffs\n",
    "print(model_k20a_coeff)\n",
    "\n",
    "# Save coeffs\n",
    "filename = './output/density_k20a_coeffs.pkl'\n",
    "pickle.dump(model_k20a_coeff, open(filename, 'wb'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Apply the new coeffs to estimate density\n",
    "k20a_rho = model_k20a_coeff[0]+(model_k20a_coeff[1]*result['force_log'])+ \\\n",
    "           (model_k20a_coeff[2]*result['force_log']*result['l'])+ \\\n",
    "           model_k20a_coeff[3]*result['l']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## K-folds OLS with outliers removed\n",
    "\n",
    "Outliers were defined as SMP/Cutter comparisons where error > than the 95th quantile in the K20a recalibration.  \n",
    "We justify this in the paper in the context of the matching procedure which cannot be assumed to be perfect.  \n",
    "Note that this removes a small number of comparisons(n=10) which are not isolated to any one profile."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Remove outliers\n",
    "result_lim = result.copy()\n",
    "result_lim['f_l'] = (result_lim['l'])*result_lim['force_log']\n",
    "result_lim['abs_error'] = np.abs(k20a_rho - result_lim['RHO'])\n",
    "q_95 = result_lim['abs_error'].quantile(0.95)\n",
    "result_lim = result_lim[result_lim['abs_error'] < q_95]\n",
    "n_removed  = len(result) - len(result_lim)\n",
    "\n",
    "print('Error threshold: %i kg m^-3' % q_95)\n",
    "print('Data points removed: %i' % n_removed)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Figure 5 with caption\n",
    "<img src=\"./output/figures/Fig05_RegressionTerms_lowres.png\" alt=\"Figure 5\" style=\"width: 800px;\"/>\n",
    "\n",
    "##### Comparison of the SMP regression parameters and corresponding snow density observations. Parameters include log-transformed median force (ln(F ̃)), microstructure length scale (L) and an interaction term (f ̃L). Relationships are separated by ice surface environment.'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Code to construct Figure 5\n",
    "f, (ax1, ax2, ax3) = plt.subplots(1, 3, sharey=True, figsize=(15,5))\n",
    "\n",
    "axs = [ax1, ax2, ax3]\n",
    "for n, ax in enumerate(axs):\n",
    "    ax.text(0.02, 0.92, '('+string.ascii_lowercase[n]+')', transform=ax.transAxes, \n",
    "            size=16, weight='bold')\n",
    "\n",
    "axis_value_size = 14\n",
    "axis_label_szie = 14\n",
    "point_size = 15\n",
    "\n",
    "ax1.tick_params(axis='both', which='major', labelsize=axis_value_size)\n",
    "ax2.tick_params(axis='both', which='major', labelsize=axis_value_size)\n",
    "ax3.tick_params(axis='both', which='major', labelsize=axis_value_size)\n",
    "\n",
    "kws_myi = dict(color = 'deepskyblue', s = point_size)\n",
    "kws_fyi = dict(color = 'black', s = point_size)\n",
    "\n",
    "ax1.scatter(result_lim[result_lim['ice_type'] == 'f']['force_log'], \n",
    "            result_lim[result_lim['ice_type'] == 'f']['RHO'], \n",
    "            label='FYI', **kws_fyi)\n",
    "ax1.scatter(result_lim[result_lim['ice_type'] == 'm']['force_log'], \n",
    "            result_lim[result_lim['ice_type'] == 'm']['RHO'], \n",
    "            label='MYI', **kws_myi)\n",
    "\n",
    "\n",
    "ax2.scatter(result_lim[result_lim['ice_type'] == 'f']['l'], \n",
    "            result_lim[result_lim['ice_type'] == 'f']['RHO'], \n",
    "            label='FYI', **kws_fyi)\n",
    "ax2.scatter(result_lim[result_lim['ice_type'] == 'm']['l'], \n",
    "            result_lim[result_lim['ice_type'] == 'm']['RHO'], \n",
    "            label='MYI', **kws_myi)\n",
    "\n",
    "ax3.scatter(result_lim[result_lim['ice_type'] == 'f']['f_l'], \n",
    "            result_lim[result_lim['ice_type'] == 'f']['RHO'], \n",
    "            label='FYI', **kws_fyi)\n",
    "ax3.scatter(result_lim[result_lim['ice_type'] == 'm']['f_l'], \n",
    "            result_lim[result_lim['ice_type'] == 'm']['RHO'], \n",
    "            label='MYI', **kws_myi)\n",
    "\n",
    "ax1.set_ylabel('Snow pit density [kg m$\\mathregular{^{-3}}$]', fontsize = axis_label_szie)\n",
    "ax1.set_xlabel(r'$ln(\\tilde{F})$ [N]', fontsize = axis_label_szie)\n",
    "ax2.set_xlabel(r'$L$ [mm]', fontsize = axis_label_szie)\n",
    "ax3.set_xlabel(r'$ln(\\tilde{F}) \\cdot L$ [N mm]', fontsize = axis_label_szie)\n",
    "\n",
    "ax1.set_xlim(-4, 4)\n",
    "ax2.set_xlim(0, 1.5)\n",
    "ax3.set_xlim(-3, 3)\n",
    "\n",
    "ax3.legend(fontsize = 12)\n",
    "ax3.legend(handletextpad=0, fontsize=12)\n",
    "\n",
    "\n",
    "f.savefig('./output/figures/Fig05_RegressionTerms_lowres.png', format='png')\n",
    "f.savefig('./output/figures/Fig05_RegressionTerms_production.pdf', format='pdf', dpi = 300)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Correlation between snow pit observed density and median force.L for all comparisons\n",
    "print(np.round(np.corrcoef(result_lim.RHO, result_lim.force_log)[0][1],2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Correlation between snow pit observed density and median force, microstructural length scale and the interaction term\n",
    "result_lim.groupby(['ice_type'])[['RHO','force_log','l','f_l']].corr()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Correlation between snow pit observed density and median force, microstructural length scale and the interaction term\n",
    "result_lim.groupby(['campaign','ice_type'])[['RHO','force_log','l','f_l']].corr()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "k_fold = KFold(10, True, RANDOM_SEED)\n",
    "rmse = []\n",
    "error = []\n",
    "r = []\n",
    "params = None\n",
    "\n",
    "for train_idx, test_idx in k_fold.split(result_lim):\n",
    "    train = result_lim.iloc[train_idx]\n",
    "    test = result_lim.iloc[test_idx]\n",
    "    \n",
    "    model_rho = ols(\"RHO ~ force_log + force_log * l + l\", train).fit()\n",
    "    predict_rho = model_rho.predict(exog=dict(force_log=test['force_log'], l=test['l']))\n",
    "    rmse = np.append(rmse, np.sqrt(np.mean((predict_rho - test['RHO'])**2)))\n",
    "    r = np.append(r,np.corrcoef(predict_rho, test['RHO'])[1][0])\n",
    "    error = np.append(error, predict_rho - test['RHO'])\n",
    "    \n",
    "    if params is None:\n",
    "        params = model_rho.params.values\n",
    "    else:\n",
    "        params =  np.vstack((params, model_rho.params.values))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# K20a evaluation stats (from kfolds method)\n",
    "k20b_rmse = rmse.mean()\n",
    "k20b_bias = error.mean()\n",
    "k20b_r2 = r.mean()**2\n",
    "\n",
    "print('K20b recalibration evaluation')\n",
    "print('N: %i' % len(result_lim))\n",
    "print('RMSE: %0.1f' % k20b_rmse)\n",
    "print('RMSE percent: %0.2f' % np.round(k20b_rmse/result_lim.RHO.mean(),2))\n",
    "print('bias: %0.1f' % k20b_bias)\n",
    "print('r^2: %0.2f' % k20b_r2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "model_k20b_coeff = [np.round(params[:,0].mean(),2), np.round(params[:,1].mean(),2),\n",
    "              np.round(params[:,3].mean(),2), np.round(params[:,2].mean(),2)]\n",
    "var_coeffs = [np.round(params[:,0].std(),2), np.round(params[:,1].std(),2),\n",
    "              np.round(params[:,3].std(),2), np.round(params[:,2].std(),2)]\n",
    "\n",
    "# Unbiased coeffs\n",
    "print(model_k20b_coeff)\n",
    "\n",
    "# Save coeffs\n",
    "filename = './output/density_k20b_coeffs.pkl'\n",
    "pickle.dump(model_k20b_coeff, open(filename, 'wb'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "k20b_rho = model_k20b_coeff[0]+(model_k20b_coeff[1]*result_lim['force_log'])+ \\\n",
    "           (model_k20b_coeff[2]*result_lim['force_log']*result_lim['l'])+ \\\n",
    "           model_k20b_coeff[3]*result_lim['l']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Figure 4 with caption\n",
    "<img src=\"./output/figures/Fig04_ModelEval_lowres.png\" alt=\"Figure 4\" style=\"width: 600px;\"/>\n",
    "\n",
    "##### Evaluation of the SMP density model parametrization of Proksch et al. (2015) (P15; Top) and recalibrated coefficients for snow on sea ice (K19b; Bottom). In both cases the model is evaluated against manual density cutter measurements of snow density."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Code to construct Figure 4\n",
    "# P15 statistics\n",
    "p2015_rmse = np.sqrt(np.mean((result['mean_samp']-result['RHO'])**2))\n",
    "p2015_bias = (result['error']).mean()\n",
    "p2015_r2 = np.ma.corrcoef(result['mean_samp'],result['RHO'])[0, 1]**2\n",
    "p2015_n = len(result['mean_samp'])\n",
    "p2015_p = stats.pearsonr(result['mean_samp'],result['RHO'])[1]\n",
    "\n",
    "\n",
    "# Fig constants\n",
    "axis_value_size = 14\n",
    "axis_label_size = 14\n",
    "rho_bin_size = 20 #in kg m-3\n",
    "line_start_a = 100\n",
    "line_end_a = 700\n",
    "line_end_b = 550\n",
    "point_size = 8\n",
    "common_bin = np.arange(line_start_a,line_end_a, rho_bin_size)\n",
    "\n",
    "\n",
    "f, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(12,10))\n",
    "f.subplots_adjust(hspace=.25)\n",
    "\n",
    "axs = [ax1, ax2, ax3, ax4]\n",
    "for n, ax in enumerate(axs):\n",
    "    ax.text(0.02, 0.9, '('+string.ascii_lowercase[n]+')', transform=ax.transAxes, \n",
    "            size=20, weight='bold')\n",
    "\n",
    "ax1.set_xlim(line_start_a,line_end_a)\n",
    "ax1.set_ylim(line_start_a,line_end_a)\n",
    "ax2.set_xlim(line_start_a,line_end_a)\n",
    "ax3.set_xlim(line_start_a,line_end_b)\n",
    "ax3.set_ylim(line_start_a,line_end_b)\n",
    "ax4.set_xlim(line_start_a,line_end_b)\n",
    "\n",
    "ax3.set_yticks(np.arange(line_start_a, line_end_b+1, 100.0))\n",
    "ax1.tick_params(axis='both', which='major', labelsize=axis_value_size)\n",
    "ax2.tick_params(axis='both', which='major', labelsize=axis_value_size)\n",
    "ax3.tick_params(axis='both', which='major', labelsize=axis_value_size)\n",
    "ax4.tick_params(axis='both', which='major', labelsize=axis_value_size)\n",
    "\n",
    "\n",
    "ax1.scatter(result['RHO'], result['mean_samp'], \n",
    "            s = point_size, color ='black', zorder = 1000)\n",
    "ax1.plot([line_start_a, line_end_a], [line_start_a, line_end_a], \n",
    "         'k-', color = 'gray' ,alpha= 0.8, zorder = 500)\n",
    "\n",
    "hist_kws = dict(density=True, bins=common_bin, histtype=\"stepfilled\", linewidth=1.25)\n",
    "ax2.ticklabel_format(axis='y',style='sci', scilimits=(1,5), useMathText=False)\n",
    "ax2.hist(result['RHO'], alpha = 1, edgecolor=\"black\", \n",
    "         color = 'grey', label = 'Pit', **hist_kws)\n",
    "ax2.hist(result['mean_samp'], alpha = 0.6, edgecolor=\"black\", \n",
    "         color = 'deepskyblue', label = 'SMP', **hist_kws)\n",
    "\n",
    "ax1.set_ylabel('P15 density [kg m$\\mathregular{^{-3}}$]', fontsize=axis_label_size)\n",
    "ax1.set_xlabel('Snow pit density [kg m$\\mathregular{^{-3}}$]', fontsize=axis_label_size)\n",
    "ax2.set_ylabel('Probability density', fontsize=axis_label_size)\n",
    "ax2.set_xlabel('Snow density [kg m$\\mathregular{^{-3}}$]', fontsize=axis_label_size)\n",
    "ax2.legend(edgecolor = 'black',  fontsize=axis_value_size)\n",
    "\n",
    "ax3.scatter(result_lim['RHO'], k20b_rho, s = point_size, color ='black',  zorder = 1000)\n",
    "ax3.plot([line_start_a, line_end_b], [line_start_a, line_end_b], \n",
    "         'k-', color = 'grey' ,alpha= 0.8,  zorder = 500)\n",
    "\n",
    "hist_kws = dict(density=True, bins = common_bin, histtype= \"stepfilled\", linewidth = 1.25)\n",
    "ax4.ticklabel_format(axis='y',style='sci', scilimits=(1,5), useMathText=False)\n",
    "\n",
    "ax4.hist(result_lim['RHO'], alpha = 1, edgecolor=\"black\", color = 'grey', label = 'Pit', **hist_kws)\n",
    "ax4.hist(k20b_rho, alpha = 0.6, edgecolor=\"black\", color = 'deepskyblue', label = 'SMP', **hist_kws)\n",
    "\n",
    "ax3.set_ylabel('K20b density [kg m$\\mathregular{^{-3}}$]',fontsize=axis_label_size)\n",
    "ax3.set_xlabel('Snow pit density [kg m$\\mathregular{^{-3}}$]',fontsize=axis_label_size)\n",
    "ax4.set_ylabel('Probability density',fontsize=axis_label_size)\n",
    "ax4.set_xlabel('Snow density [kg m$\\mathregular{^{-3}}$]',fontsize=axis_label_size)\n",
    "\n",
    "# Display stats\n",
    "ax1.text(550, 150,'N: %i \\nRMSE: %i \\nR$^2$: %0.2f'%(p2015_n, p2015_rmse, p2015_r2), fontsize=12)\n",
    "ax3.text(440, 130,'N: %i \\nRMSE: %i \\nR$^2$: %0.2f'%(len(result_lim), k20b_rmse, k20b_r2),  fontsize=12)\n",
    "\n",
    "f.savefig('./output/figures/Fig04_ModelEval_lowres.png', format='png')\n",
    "f.savefig('./output/figures/Fig04_ModelEval_production.pdf', format='pdf', dpi = 300)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Error metrics\n",
    "def rmse(data):\n",
    "    return np.sqrt(np.mean(data**2))\n",
    "\n",
    "result_lim['model_rho'] = k20b_rho\n",
    "result_lim['abs_error'] = np.abs(k20b_rho-result_lim['RHO']).values\n",
    "result_lim['error'] = (k20b_rho-result_lim['RHO']).values\n",
    "\n",
    "# Error by layer type\n",
    "np.round(result_lim.groupby('TYPE')['error'].apply(rmse)/result_lim.groupby('TYPE')['model_rho'].mean(),3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.round(result_lim.groupby('campaign')['error'].apply(rmse)/result_lim.groupby('campaign')['model_rho'].mean(),3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Overall error in %\n",
    "np.round((rmse(result_lim['error'])/result_lim['model_rho'].mean()),3)"
   ]
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}