Assignment_3.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "fed53dcc",
   "metadata": {},
   "source": [
    "# Assignment 3: When will summer Arctic sea ice disappear?\n",
    "In Part 1 we will fit a linear model to Arctic sea ice extent data. This will be a relatively straightforward task because all you're asked to do is re-use a few lines of code from Assignment 2 and then interpret the results. \n",
    "\n",
    "Sources. This part of the assignment uses data from NSIDC: https://nsidc.org/data/nsidc-0051."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "433b5a68",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the data\n",
    "import pickle\n",
    "import pandas as pd\n",
    "\n",
    "with open('northern-hemisphere-sea-ice-nsidc.pickle', 'rb') as file:\n",
    "    sea_ice_data = pickle.load(file)\n",
    "    sea_ice_df = pd.DataFrame(sea_ice_data)\n",
    "    sea_ice_df['years'] = 1970+sea_ice_df['date'].astype(int)/1e9/86400/365\n",
    "    \n",
    "# Extract just the minima\n",
    "from scipy.signal import find_peaks\n",
    "inds = find_peaks(-sea_ice_df['extent'],prominence=1)[0]\n",
    "min_dates = sea_ice_df.iloc[inds]['years'] # use this variable as the first input to the linear regression\n",
    "min_extents = sea_ice_df.iloc[inds]['extent'] # use this variable as the second input to the linear regression\n",
    "\n",
    "# Plot the data\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "fig,ax=plt.subplots(figsize=(16,9))\n",
    "plt.plot(sea_ice_df['years'],sea_ice_df['extent'],'.',label='Data')\n",
    "plt.plot(sea_ice_df.iloc[inds]['years'],sea_ice_df.iloc[inds]['extent'],'o',markersize=10,label='Summertime minima')\n",
    "plt.plot((1970,2030),(0,0),'-r',label='zero sea ice')\n",
    "plt.ylim([-1,17.5])\n",
    "plt.xlim([1978,2020])\n",
    "plt.ylabel('Sea ice extent, million sq km')\n",
    "plt.legend()\n",
    "plt.grid()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "efaddff2",
   "metadata": {},
   "source": [
    "### Fitting the linear model.\n",
    "In this part, you want to fit a linear model between to the time series of minimum extents.  The individual data points consist of pairs of (time, extent) as stored int he variables ```min_dates, min_extents```.\n",
    "\n",
    "In the cell below, reuse the code from Assignment 2 to fit a straight line to the data.  This will consist of three lines of code that import the linregres function, define the variable ```output```, and then print the output.\n",
    "#### Write code here:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "74fe05fe",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7ff9b8df",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import linregress\n",
    "output = linregress( min_dates, min_extents )\n",
    "output"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6205fc00",
   "metadata": {},
   "source": [
    "### Interpreting the results.\n",
    "Once you've correctly carried out the linear regression in the cell above, you'll be left with a ```LinregressResult``` that has properties like ```slope```, ```intercept```, etc.  The linear model describes ```extent = slope*year + intercept```.\n",
    "\n",
    "### What year does this linear model predict that the Arctic will be, on average, sea ice free?\n",
    "Solving the above equation, this can be found by evaluating the expression ``` - output.intercept/output.slope```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9c964bc7",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "48f9ead9",
   "metadata": {},
   "source": [
    "### The following code calculates confidence intervals. \n",
    "You don't need to modify this code, but you do need to extract the upper bound and lower bounds for analysis. Hint: look at the variable names."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1209c563",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate upper and lower bounds\n",
    "import numpy as np\n",
    "import scipy.stats as stats\n",
    "\n",
    "model_years = np.linspace(1970, 2100, 100)\n",
    "expected_model = output.intercept + output.slope * model_years\n",
    "\n",
    "# Confidence interval\n",
    "alpha = 0.05 # significance level\n",
    "tails = 2\n",
    "n = len(min_dates)\n",
    "dof = n - 2\n",
    "t_critical = stats.t.ppf(1 - (alpha / tails), dof)\n",
    "x = min_dates\n",
    "ci = 100*t_critical * output.stderr * np.sqrt(1 / n + (model_years - np.mean(x))**2 / np.sum((x - np.mean(x))**2))\n",
    "\n",
    "upper_bound = expected_model + ci\n",
    "lower_bound = expected_model - ci\n",
    "\n",
    "upper_bound_year = round(model_years[np.where(upper_bound<0)].min())\n",
    "lower_bound_year = round(model_years[np.where(lower_bound<0)].min())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "040048cd",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "plt.subplots(figsize=(16,9))\n",
    "plt.plot(min_dates, min_extents,'o')\n",
    "\n",
    "# Line of best fit\n",
    "plt.plot(model_years, expected_model, \n",
    "         label=f'Line of Best Fit, R² = {-output.rvalue:.2f}')\n",
    "\n",
    "plt.fill_between(\n",
    "    model_years, upper_bound, lower_bound, facecolor='#b9b9b9', zorder=0,\n",
    "    label=r'95% Confidence Interval'\n",
    ")\n",
    "\n",
    "plt.plot((1970,2200),(0,0),'-')\n",
    "plt.plot(expected_year,0,'o',label='Expected first year with no sea ice')\n",
    "plt.plot(lower_bound_year,0,'o',label='95% Confidence lower bound')\n",
    "plt.plot(upper_bound_year,0,'o',label='95% Confidence upper bound')\n",
    "plt.ylim([-1,10])\n",
    "plt.xlim([1978,2100])\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3590c1bf",
   "metadata": {},
   "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.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}