historic_analysis.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from glob import glob\n",
    "from netCDF4 import Dataset\n",
    "import matplotlib.pyplot as plt\n",
    "plt.style.use('fivethirtyeight')\n",
    "from skimage.graph import route_through_array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load Sea/Land fixed data to create a mask\n",
    "nc_f = 'data/sftlf_ARC-22_CCCma-CanESM2_rcp45_r1i1p1_CCCma-CanRCM4_r2_fx.nc'\n",
    "nc = Dataset(nc_f, 'r')\n",
    "#nc.variables.keys()\n",
    "sftlf = nc['sftlf'][0, :, :]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get historical transit years\n",
    "hist_transits = pd.read_html('http://www.nauticapedia.ca/Articles/NWP_Fulltransits.php')[0]\n",
    "# skip where transit was only by an icebreaker, stress easier travels\n",
    "no_break = hist_transits[~hist_transits['Vessel Type'].isna()]\n",
    "no_break = no_break[~no_break['Vessel Type'].str.contains('break')]\n",
    "no_break_transit_years = no_break.Year.str.split().str.get(0)\n",
    "no_break_transit_years = no_break_transit_years[no_break_transit_years.str.isalnum()].astype(np.int32).unique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# files downloaded from http://climate-modelling.canada.ca/climatemodeldata/canrcm/CanRCM4/ARC-22_CCCma-CanESM2_historical/day/atmos/sic/index.shtml\n",
    "#hist_files = glob('data/historic/day/1/*.nc')\n",
    "hist_files = glob('data/historic/day/*.nc')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "# process years\n",
    "\n",
    "weights = []\n",
    "indices_len = []\n",
    "times = []\n",
    "\n",
    "for hist_file in hist_files:\n",
    "    \n",
    "    nc = Dataset(hist_file, 'r')\n",
    "    print(hist_file, nc['time'].shape)\n",
    "    times.extend(list(nc['time'][:]))\n",
    "    \n",
    "    for i in range(nc.variables['time'].shape[0]):\n",
    "        sic = nc['sic'][i, :, :]\n",
    "        sic[np.where(sftlf!=0)] = -np.inf\n",
    "        indices, weight = route_through_array(sic, [0,116], [265,231], geometric=True, fully_connected=True)\n",
    "        weights.append(weight)\n",
    "        indices_len.append(len(indices))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.DataFrame([weights, indices_len, times]).T\n",
    "df.columns = ['weights', 'num_index', 'time']\n",
    "df['date'] = pd.to_datetime(df.time, unit='D', origin='1949-12-01')\n",
    "df['year'] = df['date'].dt.year\n",
    "df['dayofyear'] = df.date.dt.dayofyear\n",
    "df['no_break_transit'] = df.year.isin(no_break_transit_years)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.sort_values('date', inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(10,8))\n",
    "ax.scatter(df[df.no_break_transit==False].date, df[df.no_break_transit==False].weights, label='Year no transit')\n",
    "ax.scatter(df[df.no_break_transit==True].date, df[df.no_break_transit==True].weights, label='Year transit')\n",
    "ax.set_ylim(0.0001, 50000)\n",
    "ax.set_yscale('log')\n",
    "ax.set_ylabel('MCP weight', labelpad=10, fontsize=20)\n",
    "ax.set_xlabel('Days', labelpad=10, fontsize=20)\n",
    "plt.legend()\n",
    "plt.savefig('plots/transit_years.png', dpi=72, pad_inches=0, bbox_inches=\"tight\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import ks_2samp, ttest_ind"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "test_res = []\n",
    "for x in np.arange(0, 60, 0.5):\n",
    "    res = []\n",
    "    for k,g in df.groupby('year'):\n",
    "        res.append([k, (g.weights<x).sum() / len(g), g.no_break_transit.iloc[0]])\n",
    "    res = pd.DataFrame(res, columns=['year', 'days_per_year_under', 'no_break_transit'])\n",
    "    t = res.loc[res.no_break_transit==True, 'days_per_year_under']\n",
    "    f = res.loc[res.no_break_transit==False, 'days_per_year_under']\n",
    "    tt_stat, tt_pval = ttest_ind(t, f)\n",
    "    ks_stat, ks_pval = ks_2samp(t,f)\n",
    "    test_res.append([x, tt_stat, tt_pval, ks_stat, ks_pval])\n",
    "test_res = pd.DataFrame(test_res, columns=['thresh', 'tt_stat', 'tt_pval', 'ks_stat', 'ks_pval'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "day_res = []\n",
    "for x in np.arange(0, 60, 0.5):\n",
    "    res = []\n",
    "    for k,g in df.groupby('year'):\n",
    "        res.append([k, (g.weights<x).sum(), g.no_break_transit.iloc[0]])\n",
    "    res = pd.DataFrame(res, columns=['year', 'days_under', 'no_break_transit'])\n",
    "    t = res.loc[res.no_break_transit==True, 'days_under']\n",
    "    f = res.loc[res.no_break_transit==False, 'days_under']\n",
    "    day_res.append([x, np.median(t), np.median(f)])\n",
    "day_res = np.array(day_res)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(8,6))\n",
    "ax.plot(test_res.thresh, test_res['tt_pval'], label='t-test')\n",
    "ax.plot(test_res.thresh, test_res['ks_pval'], label='ks-test')\n",
    "ax.set_ylabel('P value', labelpad=10, fontsize=20)\n",
    "ax.set_xlabel('Threshold', labelpad=10, fontsize=20)\n",
    "ax.set_yscale('log')\n",
    "plt.legend()\n",
    "plt.savefig('plots/threshold_test.png', dpi=72, pad_inches=0, bbox_inches=\"tight\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(8,6))\n",
    "ax.plot(day_res[:, 0], day_res[:,1], label='transit')\n",
    "ax.plot(day_res[:, 0], day_res[:,2], label='non-transit')\n",
    "ax.set_ylabel('Median Days Under Threshold', labelpad=10, fontsize=20)\n",
    "ax.set_xlabel('Threshold', labelpad=10, fontsize=20)\n",
    "plt.legend(loc=\"lower right\")\n",
    "plt.savefig('plots/threshold_days.png', dpi=72, pad_inches=0, bbox_inches=\"tight\")\n",
    "plt.show()"
   ]
  }
 ],
 "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}