geol3910_pythonintro_globalchange_ANSWERKEY.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# GEOL 3910: Understanding Global Change through Isotope Data and Python Data Operations\n",
    "\n",
    "## *ANSWER KEY* \n",
    "\n",
    "Authors: Amanda Alexander, Tristan Caro, Kevin Rozmiarek\n",
    "\n",
    "Due for GEOL 5042 November 10, 2020"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Learning Goals\n",
    "1. introduction to pandas data frames, numpy in python\n",
    "2. import and use of scipy\n",
    "3. fitting functions for datasets with least-squares method\n",
    "4. describe anthropogenic causes of incresed carbon emissions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Tip:** Before we get started, go ahead and click 'Run' -> 'Run all cells' from the navigation tool bar up top. Later, when you are ready to run any cell you have modified, simply do <code>shift+enter</code> in the cell box or click 'Run' -> 'Run selected cells'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Table of Contents\n",
    "1. **[Stable Isotope Fractionation](#StableIsotopeFractionation1)**\n",
    "\n",
    "    1 [Background](#StableIsotopeFractionation2)\n",
    "    \n",
    "    2 [Concept Questions](#StableIsotopeFractionation3)\n",
    "    \n",
    "\n",
    "2. **[Data Structures and Visualization](#DataStructures1)**\n",
    "\n",
    "    1 [Introduction](#DataStructures2)\n",
    "    \n",
    "    2 [Example](#DataStructures3)\n",
    "    \n",
    "    3 [Practice Questions](#DataStructures4)\n",
    "    \n",
    "\n",
    "3. **[Linear Least Squares Regression Method](#LeastSquares1)**\n",
    "\n",
    "    1 [Background](#LeastSquares2)\n",
    "    \n",
    "    2 [Practice Questions](#LeastSquares3)\n",
    "    \n",
    "    3 [Example with SciPy](#LeastSquares4)\n",
    "    \n",
    "\n",
    "4. **[Summary and big picture example: Stable Isotopes of the Carbon Cycle](#BigPicture1)**\n",
    "\n",
    "    1 [Summary Questions](#BigPicture2)\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a name=\"StableIsotopeFractionation1\"></a>\n",
    "## 1. Stable Isotope Fractionation "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a name=\"StableIsotopeFractionation2\"></a>\n",
    "### 1.1 Background\n",
    "All elements consist of protons, neutrons, and electrons. Elements are defined by their *atomic number*, the number of protons contained within their nucleus, such as Oxygen, whose atomic number is **8**. Often, elements can exist with the same atomic number, but a different number of neutrons. An atom with the same atomic number, but a differing number of neutrons, is defined as an **isotope**. Often, isotopes are radioactive and cause atomic decay. Carbon 14 $^{14}C$, for example, is Carbon 12 with two extra neutrons, leading it to have an atomic mass of 14. $^{14}C$ is radioactive and decays into Nitrogen 14 via [Beta Decay](https://www.esrl.noaa.gov/gmd/ccgg/isotopes/decay.html#:~:text=C%20decays%20by%20a%20process,rather%20than%20a%20carbon%20atom.). Atoms with differing numbers of neutrons that are *not* radioactive are called stable isotopes. $^{13}C$, for example, is a stable isotope of carbon.\n",
    "\n",
    "Because stable isotopes of a given element have differing weights, they are often termed \"light\" or \"heavy.\" This weight difference is often small (just the weight of one or two neutrons), but can be significant enough to cause changes in the reaction characteristics of molecules, especially when a reaction is mediated by biological enzymes. These differing reaction rates often lead to one isotope being preferentially reacted when compared to the other - this is known as **stable isotope fractionation**. Often times, but not always, biological systems \"prefer\" the light isotope of a given element, due to it having faster reaction kinetics than its heavy counterpart. In shorthand, scientists will often say \"life prefers light.\" The result of this is that one can often search for the source of a compound of interest by looking at its isotopic composition *i.e.* does the isotopic composition of our compound of interest match that of a biological system or that of an abiotic reaction?\n",
    "\n",
    "### Commonly Studied Stable Isotopes\n",
    "A few commonly studied stable isotopes in biogeochemical systems are:\n",
    "- $^2H$ : Deuterium, stable isotope of Hydrogen $^1H$\n",
    "- $^{13}C$ : Carbon 13, stable isotope of Carbon $^{12}C$\n",
    "- $^{15}N$ : Nitrogen 15, stable isotope of Nitrogen $^{14}N$\n",
    "- $^{17}O$ : Oxygen 17, stable isotope of Oxygen $^{16}O$\n",
    "- $^{18}O$ : Oxygen 18, stable isotope of Oxygen  $^{16}O$\n",
    "\n",
    "### Delta Notation\n",
    "Stable isotope fractionation is often written out using \"delta\" ($\\delta$) notation. Delta notation is a ratio of ratios, where the ratios are those of one isotope to another. In other words:\n",
    "\n",
    "$$\n",
    "\\delta = \\frac{R_{sample} - R_{standard}}{R_{standard}}\n",
    "$$\n",
    "\n",
    "where $R_{sample}$ is the isotope ratio of the sample and  $R_{standard}$ is the isotope ratio of an international standard. In the case of Carbon, a common international standard is VPDB (Vienna Pee Dee Belemnite).\n",
    "\n",
    "Delta notation is often expressed in \"permil\". Similar to percent, delta values can be converted to permil by multiplying by 1000.\n",
    "\n",
    "$$\n",
    "\\delta \\space in \\space permil: = \\frac{R_{sample} - R_{standard}}{R_{standard}} * 1000\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.2 Concept questions\n",
    "1. Hydrogen isotopes $^1H$ and $^2H$ often fractionate more strongly in biological systems than those of Carbon $^{12}C$ and $^{13}C$. Why do you suspect this is the case?\n",
    "\n",
    "> Answer: Think about the relative weights of these elements. The stable isotopes of Carbon have atomic masses of 13 and 12, a difference of roughly 8%. Deuterium (2H) is twice as heavy as light hydrogen (1H). Because an extra neutron adds a greater amount of *relative* mass, deuterium often fractionates much more strongly than heavier elements such as Carbon.\n",
    "\n",
    "2. Consider the made-up element *Bozium* that has two isotopes of heavy and light. What is the delta value of a sample whose istotope ratio (heavy/light) is .25, compared to an international Bozium standard that has an isotope ratio of .1?\n",
    "\n",
    "$$\n",
    "\\delta = \\frac{R_{sample} - R_{standard}}{R_{standard}} \\\\\n",
    "\\delta = \\frac{0.25 - 0.1}{0.1} \\\\\n",
    "\\delta = 1.5\n",
    "$$\n",
    "\n",
    "> Answer: 1.5. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a name=\"DataStructures1\"></a>\n",
    "## 2. Data Structures and Visualization\n",
    "\n",
    "We will be using python to read, modify, and visualize our data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a name=\"DataStructures2\"></a>\n",
    "### 2.1 Introduction to Pandas Dataframes – First Steps\n",
    "There are three packages, or libraries, that are important for data analysis in python. These are \"numpy\", \"pandas\", and \"matplotlib\".\n",
    "\n",
    "You can import these into a python script using the following commands:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "pd, np, and plt are shorthand that allow us to call functions from these packages without referring to their full names."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a name=\"DataStructures3\"></a>\n",
    "### 2.2 Introduction to Pandas Dataframes – Example\n",
    "For brevity, we will only go through basic Pandas and Matplotlib syntax. For a more in-depth discussion on numpy, pandas, and matplotlib:\n",
    "- Numpy: https://numpy.org/devdocs/user/quickstart.html\n",
    "- Pandas: https://pandas.pydata.org/pandas-docs/stable/user_guide/10min.html\n",
    "- Matplotlib: https://matplotlib.org/tutorials/introductory/pyplot.html"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Pandas\n",
    "There are two main components of data in the `pandas` universe: series and Data Frames. A series is a list of numbers -- like a column. A Data Frame is multiple series, a collection of columns.\n",
    "\n",
    "For example, we can define some data using a python dictionary:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = {\n",
    "        'customer': [1,2,3,4], # a series\n",
    "        'bananas': [5,12,2,4], # a series\n",
    "        'pears': [5,6,7,8] # a series\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And tell pandas to turn it into a `DataFrame`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "inventory = pd.DataFrame(data)\n",
    "inventory"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The leftmost column is the index, and the two series we defined above are placed side-by-side. Note that in python, data is indexed starting at `0`: the first element of the list is in the \"0th\" position.\n",
    "\n",
    "We can now plot this data using matplotlib. Matplotlib can plot directly from pandas dataframes as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "inventory.plot(kind = 'scatter', x = 'customer', y = 'bananas')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a name=\"DataStructures4\"></a>\n",
    "### 2.3 Practice Questions \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Using a python dictionary, write the temperature forecast for every day next week. HINT: Days of the week should be written as strings in quotes, within a python list e.g. `[\"Monday\", \"Tuesday\", ...]`\n",
    "\n",
    "> Answer:\n",
    "```\n",
    "weather = {\n",
    "    \"day\": [\"Monday\", \"Tuesday\", \"Wednesday\", \"Thursday\", \"Friday\", \"Saturday\", \"Sunday\"],\n",
    "    \"temp (F)\": [70, 69, 31, 66, 55, 89, 87]\n",
    "}\n",
    "```\n",
    "\n",
    "2. Plot the forecast as a scatterplot using matplotlib\n",
    "> See code snippet below\n",
    "\n",
    "3. Import the data that we will use for the exercise and inspect it. What do the series consist of? How is the data structured? HINT: use the command `co2ppm_Data =  pd.read_csv('data/co2_mlo_surface-flask_1_ccgg_event.txt',delim_whitespace=True)` to import the data.\n",
    "\n",
    "> The dataframe is 9931 rows by 27 columns. Each row represents an *observation* and each column represents a *variable*. Some data is categorical, meaning that it can fit into discrete bins e.g. year, sample method, etc. Categorical data is defined by *categories*. Some data is continuous, meaning it sits on a numerical scale. Continuous data us numerical e.g. sample latitude, longitude, and elevation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Question 2: Plotting Weather Data\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "weather = {\n",
    "    \"day\": [\"Monday\", \"Tuesday\", \"Wednesday\", \"Thursday\", \"Friday\", \"Saturday\", \"Sunday\"],\n",
    "    \"temp (F)\": [70, 69, 31, 66, 55, 89, 87]\n",
    "}\n",
    "weather_df = pd.DataFrame(weather)\n",
    "weather_df\n",
    "\n",
    "# Plot\n",
    "weather_df.plot(kind = 'line', x = 'day', y = 'temp (F)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Question 3: Reading in real data\n",
    "co2ppm_Data =  pd.read_csv('data/co2_mlo_surface-flask_1_ccgg_event.txt',delim_whitespace=True)\n",
    "co2ppm_Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a name=\"LeastSquares1\"></a>\n",
    "## 3. Least Squares"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a name=\"LeastSquares2\"></a>\n",
    "### 3.1 Background"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider the data that you read in and plotted in the previous section. You could probably imagine drawing a line by eye through the majority of the points. But it is more accurate to calculate the line of best fit using the least squares method. The least squares regression method is an important application for fitting data. In this exercise, we will focus on the linear least squares. \n",
    "\n",
    "Recall the equation of a line:\n",
    "$$y = mx+b$$\n",
    "where $y$ is the dependent variable, $m$ is the slope, $x$ is the independent variable and $b$ is the y-intercept.\n",
    "\n",
    "Now we want to compute the line that fits our data via the following steps:\n",
    "\n",
    "1. For each (x,y) point calculate x2 and xy\n",
    "\n",
    "2. Sum all x, y, x2 and xy ($\\sum(x)$, $\\sum(y)$, $\\sum(xy)$ and $\\sum(x^2)$)\n",
    "\n",
    "3. Calculate Slope (m) as the following:\n",
    "\n",
    "$$m = \\frac{N\\sum{(xy)} − \\sum{x}\\sum{y}}{N \\sum{(x^2)} − (\\sum{x^{2}})}$$\n",
    "\n",
    "(Note: **N** is the number of points.)\n",
    "\n",
    "4. Calculate the y-intercept (b):\n",
    "\n",
    "$$b = \\frac{\\sum{y} − m\\sum{x}}{N}$$\n",
    "\n",
    "5. Assemble the equation of a line\n",
    "$$ y = mx + b $$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An important caution with this method is to watch out for outliers. A way to visualize this method is to look at the scatter plot of your data, and then imagine a spring connects each data point to the line of best fit, so any outliers or strange values will pull the line towards it. A spring naturally wants to be less stretched and so it turns out that minimizing the overall energy in the springs is actually equivalent to fitting a regression line using the method of least squares.\n",
    "\n",
    "<img src=\"lsm_img.png\" alt=\"A visualization of the errors (residuals) with LSM\" width=\"500\" height=\"600\">\n",
    "\n",
    "This is because the least squares method works by making the total of the square of the errors as small as possible. When we square each of the errors and then sum them all up, that total is as small as possible -- in a perfect world, the sum would be zero. Neat, huh?\n",
    "\n",
    "Note: these errors are also sometimes called *residuals*. To infer how well a line fits data: the better the line fit, the smaller the residuals. \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a name=\"LeastSquares3\"></a>\n",
    "### 3.2 Practice with least squares"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Making up some data\n",
    "np.random.seed(12345678) #THIS CAN BE CHANGED TO THE PREVIOUS DATA STRUCTURE EXAMPLES\n",
    "x = np.random.random(10)\n",
    "y = 1.6*x + np.random.random(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Using the data we just made up above, follow the five steps of the least squares method to assemble the equation of a line. Enter code in the boxes provided. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#1. for each point compute x^2 and xy:\n",
    "x2  = x**2\n",
    "xy  = x*y "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#2. sum all x, y, x2 and xy -- let's use the numpy sum() function\n",
    "sum_X  = np.sum(x)\n",
    "sum_Y  = np.sum(y)\n",
    "sum_X2 = np.sum(x2)\n",
    "sum_XY = np.sum(xy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#3. calculate slope\n",
    "# recall N is the number of points, to get this use the len() function\n",
    "N = len(x)\n",
    "\n",
    "# it might be easiest to break this into a few parts\n",
    "# first, let's work the numerator:\n",
    "numerator   = N * sum_XY - (sum_X*sum_Y)\n",
    "\n",
    "# next, the denomenator:\n",
    "denomenator = N * sum_X2 - sum_X**2\n",
    "\n",
    "# finally, the fraction as a whole:\n",
    "slope = numerator/denomenator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 4. calculate y-intercept\n",
    "y_int = (sum_Y*slope - sum_X) / N"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 5. assemble the equation in the form y = mx + b\n",
    "line_of_fit = slope*x + y_int"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. Plot the made up data as a scatter plot and add the line of best fit as a RED line. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#create basic scatterplot\n",
    "plt.scatter(x, y)\n",
    "\n",
    "#then we add the line of fit, colored red\n",
    "plt.plot(x,line_of_fit, color='r')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a name=\"LeastSquares4\"></a>\n",
    "### 3.3 Scipy example\n",
    "\n",
    "So now you've built your own least squares method and used it to fit to your data and you've plotted it and that was awesome but implimenting it was kinda a pain wasn't it? Now that you've gone through this process, you can appreciate the value of programs that do that heavy lifting for you.\n",
    "\n",
    "[Scipy](https://www.scipy.org/) is a self-described \"ecosystem\" of packages and libraries that have a great number of the tools that you will use as a scientist/geologist. In the software world, we would call this a **stack** or series of interdependent programs that build on themselves to perform tasks. Scipy is open source meaning that it is maintained by passionate community memebers invested in it's success. We are going to utilze Scipy to do what we did above in a few lines."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Import the library\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np #Can be removed because it will be used previously\n",
    "from scipy import stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Making up some data\n",
    "np.random.seed(12345678) #THIS CAN BE CHANGED TO THE PREVIOUS DATA STRUCTURE EXAMPLES\n",
    "x = np.random.random(10)\n",
    "y = 1.6*x + np.random.random(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Doing all that work in ONE line\n",
    "slope, intercept, r_value, p_value, std_err = stats.linregress(x, y)\n",
    "print(\"slope: %f    intercept: %f\" % (slope, intercept))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"R-squared: %f\" % r_value**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#Using our new plotting skills, lets plot our original data vs our fit\n",
    "plt.plot(x, y, 'o', label='original data')\n",
    "plt.plot(x, intercept + slope*x, 'r', label='fitted line')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pretty easy huh?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's look at a big picture example..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a name=\"BigPicture1\"></a>\n",
    "## 4. Summary and big picture example: Stable Isotopes of the Carbon Cycle\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a name=\"BigPicture2\"></a>\n",
    "### 4.1 Summary questions\n",
    "\n",
    "So a lot of you have probably heard of this global warming thing by now. We release carbon dioxide and other **greenhouse gases** into the atmosphere and cause solar radiation reflected off the surface of the Earth to remain trapped as heat. Lets use a combination of all of our sweat new tools to infer something about what is going on in this **carbon cycle** and where it is going.\n",
    "\n",
    "We are going to pull from the NOAA Global Monitoring Divison's flask network. This network represents dozens of different measurements of atmospheric gases over dozens of locations around the world. For now, we are going to looks at the carbon dioxide mole fraction of air (part per million to be exact) at Mauna Loa, an active volcano in Hawai'i. You'll see, quite depressingly, that it continues to go upwards. The periodic component of this plot is the seasonal cycle and the total plant biomass' ability to fix carbon out of the atmosphere. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note: you must run the following box twice for jupyter notebook to render the plot correctly :)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Using pandas, we are going to read in the text file delimited by whitespace\n",
    "co2ppm_Data =  pd.read_csv('data/co2_mlo_surface-flask_1_ccgg_event.txt',delim_whitespace=True)\n",
    "\n",
    "#Lets move the data into a form that can be plotted\n",
    "co2ppm_Dates_col = [co2ppm_Data['sample_year'],co2ppm_Data['sample_month'],co2ppm_Data['sample_day']]\n",
    "co2ppm_Dates = pd.concat(co2ppm_Dates_col, axis=1)\n",
    "co2ppm_Dates.columns = ['year', 'month','day']\n",
    "co2ppm_DateTime = pd.to_datetime(co2ppm_Dates)\n",
    "\n",
    "#The CO2 value\n",
    "co2ppm_co2 = co2ppm_Data['analysis_value']\n",
    "\n",
    "#Plot!\n",
    "plt.plot(co2ppm_DateTime, co2ppm_co2, 'o')\n",
    "plt.rcParams['figure.figsize'] = [20,10]\n",
    "plt.rcParams.update({'font.size':22})\n",
    "plt.ylim(300, 430)\n",
    "plt.xlabel('Date')\n",
    "plt.ylabel('CO2 (ppm)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Drawing from your newfound knowledge of stable isotopes, let's talk about the stable isotopes of carbon. In general, faster reactions fractionate more than slower ones. That means that quick reactions with carbon in their products will drastically prefer carbon 12 to carbon 13. Combustion, or a process of burning fuel, is one of these quick processes.\n",
    "\n",
    "**Question:** How would you expect the isotopic value of the atmosphere to change if we (as in us pesky humans) produced more combustion CO2 than other sources?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Answer:** The isotope value would most likely become more depleted in 13C. The isotope value would get more negative."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below is the isotope record at Mauna Loa. This just so happens to produced at CU's Insitute of Arctic and Alpine Research."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Same as before, just with a different data file\n",
    "co2c13_Data =  pd.read_csv('data/co2c13_mlo_surface-flask_1_sil_event.txt',delim_whitespace=True)\n",
    "\n",
    "co2c13_Dates_col = [co2c13_Data['sample_year'],co2c13_Data['sample_month'],co2c13_Data['sample_day']]\n",
    "co2c13_Dates = pd.concat(co2c13_Dates_col, axis=1)\n",
    "co2c13_Dates.columns = ['year', 'month','day']\n",
    "co2c13_DateTime = pd.to_datetime(co2c13_Dates)\n",
    "\n",
    "co2c13_co2 = co2c13_Data['analysis_value']\n",
    "\n",
    "\n",
    "plt.plot(co2c13_DateTime, co2c13_co2, 'o')\n",
    "plt.rcParams['figure.figsize'] = [20,10]\n",
    "plt.rcParams.update({'font.size':22})\n",
    "plt.ylim(-9, -7.4)\n",
    "plt.xlabel('Date')\n",
    "plt.ylabel('CO2 13C')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question:** Speculate about why this graph trend downward based on your knowledge of isotopes and the sources of carbon. What does carbon isotopes having a seasonal cycle imply about plants?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Answer:** It implies that there is a fractionation associated with photosynthesis because isotopes changes with the seasons (and by extension, plant biomass)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now lets use our awesome skills of linear regression to make a prediction of where the isotope signal may be going."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use scipy to fit a linear regression of the isotope data and plot it against the data (revisit section 3.2.2 for help).\n",
    "\n",
    "#This line should help\n",
    "import matplotlib.dates as mdates\n",
    "x = mdates.date2num(co2c13_DateTime)\n",
    "\n",
    "# You can use this to go back for plotting -> your_var_name = mdates.num2date(points)\n",
    "\n",
    "#Use your new regression to predict what the value will be in 2020. What about 2024?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This isotope data set isn't up to date for a reason! Let's find the data ourselves. Go to the data visualizer for [NOAA GMD](https://www.esrl.noaa.gov/gmd/dv/iadv/), find Mauna Loa, hit \"Carbon Cycle Gases\" and then \"Time Series\". Select the Carbon-13/Carbon-12 parameter and hit submit.\n",
    "\n",
    "**Question:** How was your prediction for 2020? Speculate on how you would fit to include the seasonal pattern to the dataset?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Answer:** It looks pretty accurate! I would fit a trig function to the plot. Some sort of fit to a function like A + Bx + Csin(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question:** What does this data set imply about human, *anthropogenic*, emissions? Understand that interpreting the isotope record of the atmosphere is an ongoing reasearch topic so don't be afraid of a wrong answer!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Answer:** The short answer is that it seems likely that the atmosphere is becoming more and more composed of anthropogenic sources of carbon."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}