forecast_analysis.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import requests\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from glob import glob\n",
    "from netCDF4 import Dataset\n",
    "from skimage.graph import route_through_array"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Process Projections"
   ]
  },
  {
   "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": [
    "# go through each rcp\n",
    "rcps = [\"rcp45\", \"rcp85\"]\n",
    "\n",
    "weights = []\n",
    "indices_len = []\n",
    "times = []\n",
    "rcp_track = []\n",
    "\n",
    "for rcp in rcps:\n",
    "    \n",
    "    # get filenames\n",
    "    site = \"http://climate-modelling.canada.ca/climatemodeldata/canrcm/CanRCM4/ARC-22_CCCma-CanESM2_{}/day/atmos/sic/index.shtml\".format(rcp)\n",
    "    file_df = pd.read_html(site)[0]\n",
    "    \n",
    "    # for through files\n",
    "    for filename in file_df['File name']:\n",
    "        print(filename)\n",
    "        \n",
    "        # check if file already downloaded\n",
    "        if not os.path.isfile('data/future/{}/{}'.format(rcp, filename)):\n",
    "            # download file\n",
    "            print('Downloading')\n",
    "            download_base = \"http://climate-modelling.canada.ca/cgi-bin/data_get/get_www_nc?path=\"\n",
    "            path = \"/data3/CORDEX/output/CCCma/CanRCM4/ARC-22_CCCma-CanESM2_{}/day/atmos/sic/r1i1p1/{}\".format(rcp, filename)\n",
    "            link = download_base + path\n",
    "            response = requests.get(link, stream=True)\n",
    "            response.raise_for_status()\n",
    "            with open('data/future/{}/{}'.format(rcp, filename), 'wb') as f:\n",
    "                for block in response.iter_content(1024):\n",
    "                    f.write(block)\n",
    "            \n",
    "        print('Processing')\n",
    "        nc = Dataset('data/future/{}/{}'.format(rcp, filename), 'r')\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",
    "            rcp_track.append(rcp)\n",
    "            indices_len.append(len(indices))\n",
    "                    \n",
    "df = pd.DataFrame([weights, indices_len, times, rcp_track]).T\n",
    "df.columns = ['weights', 'num_index', 'time', 'rcp']\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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.to_csv('forecast_results.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Analyze Results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.style.use('fivethirtyeight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv('data/forecast_results.csv', index_col=[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res_year = df.groupby(['year', 'rcp']).weights.apply(lambda x : (x<5).sum()).reset_index()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res_year.sort_values('year', inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(10,8))\n",
    "ax.axhline(32, linestyle='--', c='black')\n",
    "ax.plot(res_year.loc[res_year.rcp=='rcp45', 'year'], res_year.loc[res_year.rcp=='rcp45', 'weights'], marker='o', linestyle='-', label='RCP4.5')\n",
    "ax.plot(res_year.loc[res_year.rcp=='rcp85', 'year'], res_year.loc[res_year.rcp=='rcp85', 'weights'], marker='o', linestyle='-', label='RCP8.5')\n",
    "ax.set_xticks(range(2000,2100,10))\n",
    "ax.set_ylabel('Days Under Threshold (5)', labelpad=10, fontsize=20)\n",
    "ax.set_xlabel('Year', labelpad=10, fontsize=20)\n",
    "plt.legend()\n",
    "#plt.savefig('plots/forecast.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
}