ESCI451_Module_3A.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction to Python for Earth Scientists\n",
    "\n",
    "These notebooks have been developed by Calum Chamberlain, Finnigan Illsley-Kemp and John Townend at [Victoria University of Wellington-Te Herenga Waka](https://www.wgtn.ac.nz) for use by Earth Science graduate students. \n",
    "\n",
    "The notebooks cover material that we think will be of particular benefit to those students with little or no previous experience of computer-based data analysis. We presume very little background in command-line or code-based computing, and have compiled this material with an emphasis on general tasks that a grad student might encounter on a daily basis. \n",
    "\n",
    "In 2021, this material will be delivered at the start of Trimester 1 in conjunction with [ESCI451 Active Earth](https://www.wgtn.ac.nz/courses/esci/451/2021/offering?crn=32176). Space and pandemic alert levels permitting, interested students not enrolled in ESCI451 are encouraged to come along too but please contact Calum, Finn, or John first.\n",
    "\n",
    "| Notebook | Contents | Data |\n",
    "| --- | --- | --- |\n",
    "| [1A](ESCI451_Module_1A.ipynb) | Introduction to programming, Python, and Jupyter notebooks | - |\n",
    "| [1B](ESCI451_Module_1B.ipynb) | Basic data types and variables, getting data, and plotting with Matplotlib | Geodetic positions |\n",
    "| [2A](ESCI451_Module_2A.ipynb) | More complex plotting, introduction to Numpy | Geodetic positions; DFDP-2B temperatures |\n",
    "| [2B](ESCI451_Module_2B.ipynb) | Using Pandas to load, peruse and plot data | Earthquake catalogue  |\n",
    "| **[3A](ESCI451_Module_3A.ipynb)** | **Working with Pandas dataframes** | **Geochemical data set; earthquake catalogue** |\n",
    "| [3B](ESCI451_Module_3B.ipynb) | Simple time series analysis using Pandas | Historical temperature records |\n",
    "| [4A](ESCI451_Module_4A.ipynb) | Making maps with Cartopy | Earthquake catalogue |\n",
    "| [4B](ESCI451_Module_4B.ipynb) | Working with gridded data | DEMs and Ashfall data |\n",
    "\n",
    "The content may change in response to students' questions or current events. Each of the four modules has been designed to take about three hours, with a short break between each of the two parts.\n",
    "\n",
    "Remember to run the first cell below to set up plotting."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib widget"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# This notebook\n",
    "\n",
    "1. More with Pandas dataframes\n",
    "   - Loading and summarising another dataset\n",
    "   - Working with subplots\n",
    "   - Regression\n",
    "   - Dates and times and datetimes\n",
    "   - Querying a dataframe using a function\n",
    "   - Applying functions to a dataframe\n",
    "2. Extension: more sophisticated analysis of a big dataset\n",
    "   \n",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Loading and summarising another dataset\n",
    "\n",
    "Let's recap where we got to with Pandas in the previous notebook. We learnt how to obtain data from GeoNet and store it in .csv file and a Pandas dataframe, which we could then inspect and plot in various ways. \n",
    "\n",
    "We'll start this notebook by reviewing some of the ways in which we can inspect, describe and plot Pandas dataframes, and then focus on working with dates and times. The first step is to reload the necessary modules."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import requests\n",
    "import matplotlib.pyplot as plt\n",
    "import datetime\n",
    "%matplotlib widget"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's load a geochemical data set obtained by a PhD student, Elliot Swallow, for volcanic rocks from Huckleberry Ridge."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "        SiO2 (wt%)  TiO2 (wt%)  Al2O3 (wt%)  Fe2O3 (T) (wt%)  MnO (wt%)  \\\n",
      "Sample                                                                    \n",
      "YP114        72.98        0.26        13.56             2.96       0.05   \n",
      "YP307        76.50        0.13        11.85             1.70       0.04   \n",
      "YP359        74.56        0.17        12.83             2.06       0.04   \n",
      "YP363        76.61        0.12        12.06             1.73       0.04   \n",
      "YP414        76.08        0.14        12.14             1.79       0.04   \n",
      "...            ...         ...          ...              ...        ...   \n",
      "YP564        76.64        0.17        12.23             2.13       0.02   \n",
      "YP603        76.93        0.19        11.77             2.29       0.04   \n",
      "YP081        75.68        0.10        12.81             1.84       0.04   \n",
      "YP133        76.69        0.10        12.04             1.55       0.02   \n",
      "YP600        76.39        0.10        12.54             1.59       0.03   \n",
      "\n",
      "        MgO (wt%)  CaO (wt%)  Na2O (wt%)  K2O (wt%)  P2O5 (wt%)  ...  \\\n",
      "Sample                                                           ...   \n",
      "YP114        0.08       0.82        4.05       5.22       0.040  ...   \n",
      "YP307        0.08       1.15        2.87       5.67       0.020  ...   \n",
      "YP359        0.23       1.28        2.95       5.86       0.030  ...   \n",
      "YP363        0.04       0.69        2.78       5.94       0.010  ...   \n",
      "YP414        0.06       0.96        2.77       5.99       0.020  ...   \n",
      "...           ...        ...         ...        ...         ...  ...   \n",
      "YP564        0.01       0.37        3.30       4.89       0.030  ...   \n",
      "YP603        0.01       0.39        3.48       4.58       0.020  ...   \n",
      "YP081        0.14       0.68        3.34       5.35       0.021  ...   \n",
      "YP133        0.25       1.45        3.03       4.85       0.020  ...   \n",
      "YP600        0.06       1.01        2.84       5.17       0.010  ...   \n",
      "\n",
      "        Dy (ppm)  Ho (ppm)  Er (ppm)  Tm (ppm)  Yb (ppm)  Lu (ppm)  Hf (ppm)  \\\n",
      "Sample                                                                         \n",
      "YP114       11.2   2.20002       6.4      0.89       5.7   0.82957       9.8   \n",
      "YP307       10.1   1.90000       6.0      0.82       5.3   0.88000       6.0   \n",
      "YP359        9.8   1.90000       5.6      0.75       4.9   0.85000       6.8   \n",
      "YP363       10.5   2.10000       6.3      0.89       5.6   1.02000       6.6   \n",
      "YP414        8.9   1.73000       5.3      0.73       4.6   0.84000       6.4   \n",
      "...          ...       ...       ...       ...       ...       ...       ...   \n",
      "YP564        9.0   1.70000       4.9      0.64       4.1   0.60000       8.1   \n",
      "YP603        7.2   1.47000       4.4      0.62       4.2   0.62000       8.2   \n",
      "YP081       14.7   3.00000       8.9      1.29       8.3   1.20000       7.8   \n",
      "YP133       11.9   2.40000       7.4      1.06       6.7   1.19000       7.0   \n",
      "YP600       12.2   2.40000       7.4      1.05       6.6   1.12000       6.9   \n",
      "\n",
      "        Pb (ppm)  Th (ppm)  U (ppm)  \n",
      "Sample                               \n",
      "YP114       25.0  19.36672      2.2  \n",
      "YP307       40.0  23.00000      4.5  \n",
      "YP359       55.0  23.00000      3.8  \n",
      "YP363       38.0  23.00000      5.1  \n",
      "YP414       28.0  20.00000      4.0  \n",
      "...          ...       ...      ...  \n",
      "YP564       38.0  21.00000      6.3  \n",
      "YP603       35.0  16.90000      4.2  \n",
      "YP081       29.0  31.00000      7.5  \n",
      "YP133       42.0  25.00000      6.4  \n",
      "YP600       80.0  25.00000      6.0  \n",
      "\n",
      "[65 rows x 38 columns]\n"
     ]
    }
   ],
   "source": [
    "geochem = pd.read_csv(\n",
    "    \"data/Edited Swallow et al J Petrol data for plotting.csv\",\n",
    "    index_col=\"Sample\")\n",
    "print(geochem)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Printing a dataframe out gives us a rough look at the data and we've seen before how we can compute various statistics to explore the dataframe piece by piece. But there's another way of getting a nicely-formatted, informative summary of the whole dataframe and that is to use the `.describe()` command:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>SiO2 (wt%)</th>\n",
       "      <th>TiO2 (wt%)</th>\n",
       "      <th>Al2O3 (wt%)</th>\n",
       "      <th>Fe2O3 (T) (wt%)</th>\n",
       "      <th>MnO (wt%)</th>\n",
       "      <th>MgO (wt%)</th>\n",
       "      <th>CaO (wt%)</th>\n",
       "      <th>Na2O (wt%)</th>\n",
       "      <th>K2O (wt%)</th>\n",
       "      <th>P2O5 (wt%)</th>\n",
       "      <th>...</th>\n",
       "      <th>Dy (ppm)</th>\n",
       "      <th>Ho (ppm)</th>\n",
       "      <th>Er (ppm)</th>\n",
       "      <th>Tm (ppm)</th>\n",
       "      <th>Yb (ppm)</th>\n",
       "      <th>Lu (ppm)</th>\n",
       "      <th>Hf (ppm)</th>\n",
       "      <th>Pb (ppm)</th>\n",
       "      <th>Th (ppm)</th>\n",
       "      <th>U (ppm)</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>count</th>\n",
       "      <td>65.000000</td>\n",
       "      <td>65.000000</td>\n",
       "      <td>65.000000</td>\n",
       "      <td>65.000000</td>\n",
       "      <td>65.000000</td>\n",
       "      <td>65.000000</td>\n",
       "      <td>65.000000</td>\n",
       "      <td>65.000000</td>\n",
       "      <td>65.000000</td>\n",
       "      <td>65.000000</td>\n",
       "      <td>...</td>\n",
       "      <td>65.000000</td>\n",
       "      <td>65.000000</td>\n",
       "      <td>65.000000</td>\n",
       "      <td>65.000000</td>\n",
       "      <td>65.000000</td>\n",
       "      <td>65.000000</td>\n",
       "      <td>65.000000</td>\n",
       "      <td>65.000000</td>\n",
       "      <td>65.000000</td>\n",
       "      <td>65.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mean</th>\n",
       "      <td>74.183231</td>\n",
       "      <td>0.332462</td>\n",
       "      <td>12.286308</td>\n",
       "      <td>3.753692</td>\n",
       "      <td>0.060615</td>\n",
       "      <td>0.125385</td>\n",
       "      <td>1.287692</td>\n",
       "      <td>3.251846</td>\n",
       "      <td>4.600769</td>\n",
       "      <td>0.071015</td>\n",
       "      <td>...</td>\n",
       "      <td>9.144615</td>\n",
       "      <td>1.793268</td>\n",
       "      <td>5.316615</td>\n",
       "      <td>0.735387</td>\n",
       "      <td>4.721538</td>\n",
       "      <td>0.780735</td>\n",
       "      <td>9.248016</td>\n",
       "      <td>31.741538</td>\n",
       "      <td>17.719730</td>\n",
       "      <td>4.264615</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>std</th>\n",
       "      <td>3.039516</td>\n",
       "      <td>0.272827</td>\n",
       "      <td>0.697152</td>\n",
       "      <td>2.692407</td>\n",
       "      <td>0.035745</td>\n",
       "      <td>0.103155</td>\n",
       "      <td>0.951845</td>\n",
       "      <td>0.369501</td>\n",
       "      <td>1.149364</td>\n",
       "      <td>0.096781</td>\n",
       "      <td>...</td>\n",
       "      <td>2.298507</td>\n",
       "      <td>0.458869</td>\n",
       "      <td>1.399712</td>\n",
       "      <td>0.195950</td>\n",
       "      <td>1.234845</td>\n",
       "      <td>0.212324</td>\n",
       "      <td>2.712692</td>\n",
       "      <td>9.352855</td>\n",
       "      <td>5.001809</td>\n",
       "      <td>0.990649</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>min</th>\n",
       "      <td>66.460000</td>\n",
       "      <td>0.100000</td>\n",
       "      <td>11.480000</td>\n",
       "      <td>1.320000</td>\n",
       "      <td>0.010000</td>\n",
       "      <td>-0.000000</td>\n",
       "      <td>0.330000</td>\n",
       "      <td>2.580000</td>\n",
       "      <td>2.040000</td>\n",
       "      <td>0.010000</td>\n",
       "      <td>...</td>\n",
       "      <td>3.200000</td>\n",
       "      <td>0.620000</td>\n",
       "      <td>1.880000</td>\n",
       "      <td>0.270000</td>\n",
       "      <td>1.800000</td>\n",
       "      <td>0.290000</td>\n",
       "      <td>6.000000</td>\n",
       "      <td>16.800000</td>\n",
       "      <td>7.700000</td>\n",
       "      <td>2.200000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>25%</th>\n",
       "      <td>72.230000</td>\n",
       "      <td>0.140000</td>\n",
       "      <td>11.760000</td>\n",
       "      <td>1.840000</td>\n",
       "      <td>0.040000</td>\n",
       "      <td>0.050000</td>\n",
       "      <td>0.700000</td>\n",
       "      <td>2.950000</td>\n",
       "      <td>3.920000</td>\n",
       "      <td>0.020000</td>\n",
       "      <td>...</td>\n",
       "      <td>7.600000</td>\n",
       "      <td>1.500000</td>\n",
       "      <td>4.400000</td>\n",
       "      <td>0.610000</td>\n",
       "      <td>3.900000</td>\n",
       "      <td>0.640000</td>\n",
       "      <td>7.000000</td>\n",
       "      <td>26.000000</td>\n",
       "      <td>14.000000</td>\n",
       "      <td>3.800000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>50%</th>\n",
       "      <td>75.060000</td>\n",
       "      <td>0.220000</td>\n",
       "      <td>12.000000</td>\n",
       "      <td>2.540000</td>\n",
       "      <td>0.050000</td>\n",
       "      <td>0.090000</td>\n",
       "      <td>1.030000</td>\n",
       "      <td>3.260000</td>\n",
       "      <td>4.960000</td>\n",
       "      <td>0.030000</td>\n",
       "      <td>...</td>\n",
       "      <td>9.600000</td>\n",
       "      <td>1.890000</td>\n",
       "      <td>5.700000</td>\n",
       "      <td>0.780000</td>\n",
       "      <td>5.000000</td>\n",
       "      <td>0.840000</td>\n",
       "      <td>8.300000</td>\n",
       "      <td>30.000000</td>\n",
       "      <td>17.900000</td>\n",
       "      <td>4.200000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>75%</th>\n",
       "      <td>76.690000</td>\n",
       "      <td>0.390000</td>\n",
       "      <td>12.760000</td>\n",
       "      <td>4.570000</td>\n",
       "      <td>0.070000</td>\n",
       "      <td>0.190000</td>\n",
       "      <td>1.400000</td>\n",
       "      <td>3.480000</td>\n",
       "      <td>5.370000</td>\n",
       "      <td>0.070000</td>\n",
       "      <td>...</td>\n",
       "      <td>10.500000</td>\n",
       "      <td>2.100000</td>\n",
       "      <td>6.200000</td>\n",
       "      <td>0.840000</td>\n",
       "      <td>5.500000</td>\n",
       "      <td>0.920000</td>\n",
       "      <td>11.000000</td>\n",
       "      <td>36.000000</td>\n",
       "      <td>21.000000</td>\n",
       "      <td>4.700000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>max</th>\n",
       "      <td>77.640000</td>\n",
       "      <td>1.120000</td>\n",
       "      <td>14.050000</td>\n",
       "      <td>10.680000</td>\n",
       "      <td>0.160000</td>\n",
       "      <td>0.470000</td>\n",
       "      <td>3.910000</td>\n",
       "      <td>4.180000</td>\n",
       "      <td>6.160000</td>\n",
       "      <td>0.390000</td>\n",
       "      <td>...</td>\n",
       "      <td>14.700000</td>\n",
       "      <td>3.000000</td>\n",
       "      <td>8.900000</td>\n",
       "      <td>1.290000</td>\n",
       "      <td>8.300000</td>\n",
       "      <td>1.200000</td>\n",
       "      <td>16.224530</td>\n",
       "      <td>80.000000</td>\n",
       "      <td>31.000000</td>\n",
       "      <td>7.500000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>8 rows × 38 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "       SiO2 (wt%)  TiO2 (wt%)  Al2O3 (wt%)  Fe2O3 (T) (wt%)  MnO (wt%)  \\\n",
       "count   65.000000   65.000000    65.000000        65.000000  65.000000   \n",
       "mean    74.183231    0.332462    12.286308         3.753692   0.060615   \n",
       "std      3.039516    0.272827     0.697152         2.692407   0.035745   \n",
       "min     66.460000    0.100000    11.480000         1.320000   0.010000   \n",
       "25%     72.230000    0.140000    11.760000         1.840000   0.040000   \n",
       "50%     75.060000    0.220000    12.000000         2.540000   0.050000   \n",
       "75%     76.690000    0.390000    12.760000         4.570000   0.070000   \n",
       "max     77.640000    1.120000    14.050000        10.680000   0.160000   \n",
       "\n",
       "       MgO (wt%)  CaO (wt%)  Na2O (wt%)  K2O (wt%)  P2O5 (wt%)  ...  \\\n",
       "count  65.000000  65.000000   65.000000  65.000000   65.000000  ...   \n",
       "mean    0.125385   1.287692    3.251846   4.600769    0.071015  ...   \n",
       "std     0.103155   0.951845    0.369501   1.149364    0.096781  ...   \n",
       "min    -0.000000   0.330000    2.580000   2.040000    0.010000  ...   \n",
       "25%     0.050000   0.700000    2.950000   3.920000    0.020000  ...   \n",
       "50%     0.090000   1.030000    3.260000   4.960000    0.030000  ...   \n",
       "75%     0.190000   1.400000    3.480000   5.370000    0.070000  ...   \n",
       "max     0.470000   3.910000    4.180000   6.160000    0.390000  ...   \n",
       "\n",
       "        Dy (ppm)   Ho (ppm)   Er (ppm)   Tm (ppm)   Yb (ppm)   Lu (ppm)  \\\n",
       "count  65.000000  65.000000  65.000000  65.000000  65.000000  65.000000   \n",
       "mean    9.144615   1.793268   5.316615   0.735387   4.721538   0.780735   \n",
       "std     2.298507   0.458869   1.399712   0.195950   1.234845   0.212324   \n",
       "min     3.200000   0.620000   1.880000   0.270000   1.800000   0.290000   \n",
       "25%     7.600000   1.500000   4.400000   0.610000   3.900000   0.640000   \n",
       "50%     9.600000   1.890000   5.700000   0.780000   5.000000   0.840000   \n",
       "75%    10.500000   2.100000   6.200000   0.840000   5.500000   0.920000   \n",
       "max    14.700000   3.000000   8.900000   1.290000   8.300000   1.200000   \n",
       "\n",
       "        Hf (ppm)   Pb (ppm)   Th (ppm)    U (ppm)  \n",
       "count  65.000000  65.000000  65.000000  65.000000  \n",
       "mean    9.248016  31.741538  17.719730   4.264615  \n",
       "std     2.712692   9.352855   5.001809   0.990649  \n",
       "min     6.000000  16.800000   7.700000   2.200000  \n",
       "25%     7.000000  26.000000  14.000000   3.800000  \n",
       "50%     8.300000  30.000000  17.900000   4.200000  \n",
       "75%    11.000000  36.000000  21.000000   4.700000  \n",
       "max    16.224530  80.000000  31.000000   7.500000  \n",
       "\n",
       "[8 rows x 38 columns]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "geochem.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How cool is that?!\n",
    "\n",
    "## Working with subplots\n",
    "\n",
    "What we can immediately see is that this dataset contains a lot of geochemical parameters (38!) and if we want to discover how they relate to each other we're going to have to do some more plotting. Ideally, it would be interesting to plot each parameter against every other one but this will get pretty overwhelming so we'll focus on a few interesting ones.\n",
    "\n",
    "To make our plots, we'll use the `plt.subplots()` syntax we've encountered previously, but which we haven't yet used to its full extent. As the name suggests, subplots let you make multiple plots in one. \n",
    "\n",
    "To start with, let's plot the first four major elements in the table ($TiO_2$, $Al_2O_3$, $Fe_2O_3 (T)$, $MnO$) against $Si0_2$.\n",
    "\n",
    "To do this we will need four axes.  We can make four axes using the `plt.subplots` command.  Lets get help for that function.  You can get to the docs for a given function by typing `function?`.  Run the cell below and you should see a pop-up of the docs.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.subplots?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The arguments that we care about at the moment at:\n",
    "- `nrows`, the number of rows of subplots\n",
    "- `ncols`, the number of columns of subplots\n",
    "- `sharex`, whether axes should have the same x-axis.  We will be having $SiO_2$ on all the x-axes, so we will set this to be `True`.\n",
    "\n",
    "Lets make a 2x2 grid of subplots:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "da609a41a1464462a8981a7d5b380e39",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "FigureCanvasNbAgg()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(nrows=2, ncols=2, sharex=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that we īncluded `sharex=True` in our command, which will ensure that all four panels have the same x-axis.\n",
    "\n",
    "This returned two things as usual, `fig`: the `figure` containing all our axes, and `ax`, which is a list of the subplot axes.  We can specify which axes we want to use for each graph by remembering that Python counts from zero and noting that _the axes are indexed by row and then by column_."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "aa856ea696444c2bb958cae8bfb9f420",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "FigureCanvasNbAgg()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(nrows=2, ncols=2, sharex=True)\n",
    "for column in range(2):\n",
    "    for row in range(2):\n",
    "        ax[row][column].set_title(f\"ax[{row}][{column}]\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Right, let's plot actual data using two loops to run through the parameters and axes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "803dd7367b684e418d57022dc5fa08b4",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "FigureCanvasNbAgg()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(nrows=2, ncols=2, sharex=True)\n",
    "\n",
    "elements = [[\"TiO2 (wt%)\", \"Al2O3 (wt%)\"], [\"Fe2O3 (T) (wt%)\", \"MnO (wt%)\"]]\n",
    "for column in range(2):\n",
    "    for row in range(2):\n",
    "        element_name = elements[row][column]\n",
    "        ax[row][column].scatter(geochem[\"SiO2 (wt%)\"], geochem[element_name])\n",
    "        ax[row][column].set_ylabel(element_name)\n",
    "        ax[row][column].set_xlabel(\"$SiO_2$ (wt%)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The y-axis labels are getting in the way a bit but we can easily avoid this as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "e2331212b9174050a054c52df010df86",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "FigureCanvasNbAgg()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(nrows=2, ncols=2, sharex=True)\n",
    "\n",
    "elements = [[\"TiO2 (wt%)\", \"Al2O3 (wt%)\"], [\"Fe2O3 (T) (wt%)\", \"MnO (wt%)\"]]\n",
    "for column in range(2):\n",
    "    for row in range(2):\n",
    "        element_name = elements[row][column]\n",
    "        ax[row][column].scatter(geochem[\"SiO2 (wt%)\"], geochem[element_name])\n",
    "        ax[row][column].set_ylabel(element_name)\n",
    "        ax[row][column].set_xlabel(\"$SiO_2$ (wt%)\")\n",
    "        if column == 1:\n",
    "            ax[row][column].yaxis.tick_right()\n",
    "            ax[row][column].yaxis.set_label_position(\"right\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise:\n",
    "\n",
    "Try writing some code to plot the four parameters of interest against $SiO_2$ in a 1x4 subplot layout with a common y-axis ($SiO_2$)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your answer here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Regression\n",
    "\n",
    "We can infer from the plots above that $Fe_2O_3$ and $TiO_2$ exhibit very similar behaviour with regard to $SiO_2$. Sure enough, we find when we plot one against the other that they have are highly correlated:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "f0c6765a00754b3eb3089b691df38299",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "FigureCanvasNbAgg()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x15293f9e3d50>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "geochem.plot(x=\"Fe2O3 (T) (wt%)\",y=\"TiO2 (wt%)\",kind='scatter')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll use these two parameters to briefly explore linear regression. There are myriad ways to do this in Python but perhaps the simplest to start with is numpy's \"polyfit' algorithm. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.09956925 -0.04129079]\n",
      "[[ 5.62105856e-06 -2.10997243e-05]\n",
      " [-2.10997243e-05  1.19322352e-04]]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "model, covariance = np.polyfit(geochem['Fe2O3 (T) (wt%)'],geochem['TiO2 (wt%)'],1, cov=True)\n",
    "print(model)\n",
    "print(covariance)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What we've done is fit a linear model of the form $y=ax+b$ to the Fe (x) and Ti (y) oxide dat, using least-squares. We could have fit a quadratic model by simply changing the \"1\" at the end of the command to \"2\".\n",
    "\n",
    "We can use the \"poly1d\" function to construct a function from the model parameters that, when re-applied to the Fe data, will return the predicted Ti measurements:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "prediction_function=np.poly1d(model)\n",
    "Ti_predicted=geochem['Fe2O3 (T) (wt%)'].pipe(prediction_function)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's plot the results to see what we've found and whether it makes sense."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "68a3390d698044e2a43294166fb318b9",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "FigureCanvasNbAgg()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/calumch/miniconda3/envs/conda_37/lib/python3.7/site-packages/matplotlib/cbook/__init__.py:1402: FutureWarning: Support for multi-dimensional indexing (e.g. `obj[:, None]`) is deprecated and will be removed in a future version.  Convert to a numpy array before indexing instead.\n",
      "  x[:, None]\n",
      "/home/calumch/miniconda3/envs/conda_37/lib/python3.7/site-packages/matplotlib/axes/_base.py:276: FutureWarning: Support for multi-dimensional indexing (e.g. `obj[:, None]`) is deprecated and will be removed in a future version.  Convert to a numpy array before indexing instead.\n",
      "  x = x[:, np.newaxis]\n",
      "/home/calumch/miniconda3/envs/conda_37/lib/python3.7/site-packages/matplotlib/axes/_base.py:278: FutureWarning: Support for multi-dimensional indexing (e.g. `obj[:, None]`) is deprecated and will be removed in a future version.  Convert to a numpy array before indexing instead.\n",
      "  y = y[:, np.newaxis]\n"
     ]
    }
   ],
   "source": [
    "fig, ax=plt.subplots()\n",
    "ax.plot(geochem['Fe2O3 (T) (wt%)'],geochem['TiO2 (wt%)'],'o',color='red',label='Measurements')\n",
    "ax.plot(geochem['Fe2O3 (T) (wt%)'],Ti_predicted,label='Linear model',linestyle='dotted',linewidth=0.5,color='red')\n",
    "ax.set_xlabel(\"Fe2O3 (T) (wt%)\")\n",
    "ax.set_ylabel(\"TiO2 (wt%)\")\n",
    "ax.set_title(\"An example of linear regression with numpy\")\n",
    "ax.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## More on dates and times and datetimes\n",
    "\n",
    "In the previous notebook, we made our GeoNet query using `str` objects for the start and end-time arguments, but Python has nice native ways of working with dates and times. These do useful things like cope with leap-years, allow you to add seconds (or other time units) to dates and times, and allow you to format dates and times as `str` objects.\n",
    "\n",
    "We should switch from giving `str`s to giving `datetime`s.  `datetime` objects come from Python's native `datetime` library, and include a handy `.strftime` method which is literally string-format-time.  We will use that to make the correctly formatted string for our query.  The query requires something of the format:\n",
    "\n",
    "> year-month-dayThour:minute:second\n",
    "\n",
    "In `datetime` speak the format string for that is:\n",
    "\n",
    "> `%Y-%m-%dT%H:%M:%S`\n",
    "\n",
    "- `%Y` is a four-digit year (e.g.: ..., 2018, 2019, 2020, ...)\n",
    "- `%m` is a two-digit month (e.g.: 01, 02, 03, 04, 05, 06, 07, 08, 09, 10, 11, 12)\n",
    "- `%d` is a two-digit day (zero-padded as months are, e.g.: 01, 02, ... 28, 29, 30, 31)\n",
    "- `T` is just a letter, anything not preceded by a `%` sign is interpreted as a `str`\n",
    "- `%H` is a two-digit hour (as above)\n",
    "- `%M` is a two-digit minute (as above)\n",
    "- `%S` is a two-digit second (as above).\n",
    "\n",
    "Other formatters, for things like day of the week, month as a word, julian-day, milliseconds, etc. can be found in the [offical docs](https://docs.python.org/3/library/datetime.html#strftime-and-strptime-format-codes).\n",
    "\n",
    "Lets see how we could format a `datetime`.  To start we need to make a `datetime` object, we can provide arguments of the year, month, day, hour, minute, second, (millisecond) to make one:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2020-01-10 12:43:10\n"
     ]
    }
   ],
   "source": [
    "# datetime is the module, datetime.datetime is the object itself\n",
    "test_time = datetime.datetime(2020, 1, 10, 12, 43, 10)\n",
    "print(test_time)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now lets format the string the way that we want it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2020-01-10T12:43:10\n"
     ]
    }
   ],
   "source": [
    "print(test_time.strftime(\"%Y-%m-%dT%H:%M:%S\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** format the `datetime` object as \"year/month/day hour:minute:seconds\" and then try some other formats using the information listed in the [offical docs](https://docs.python.org/3/library/datetime.html#strftime-and-strptime-format-codes)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your answer here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we know how it all works, lets put it together into a function with some useful\n",
    "arguments that we can query.\n",
    "\n",
    "## Querying a dataframe using a function\n",
    "\n",
    "We will set some default values for our arguments, so that we do not always have to\n",
    "specify every argument.  These defaults are given in the function definition as:\n",
    "\n",
    "```python\n",
    "def function(argument=value, ...):\n",
    "    contents_of_function\n",
    "```\n",
    "where `argument` is the argument name, and `value` is the default value for that argument:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_geonet_quakes(\n",
    "    min_latitude=-49.0, max_latitude=-40.0,\n",
    "    min_longitude=164.0, max_longitude=182.0,\n",
    "    min_magnitude=1, max_magnitude=9.0,\n",
    "    min_depth=0.0, max_depth=500.0,\n",
    "    start_time=datetime.datetime(1960, 1, 1),\n",
    "    end_time=datetime.datetime(2020, 1, 1),\n",
    "):\n",
    "    \"\"\"\n",
    "    Get a dataframe of the eatrhquakes in the GeoNet catalogue.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    min_latitude\n",
    "        Minimum latitude in degrees for search\n",
    "    max_latitude\n",
    "        Maximum latitude in degrees for search\n",
    "    min_longitude\n",
    "        Minimum longitude in degrees for search\n",
    "    max_longitude\n",
    "        Maximum longitude in degrees for search\n",
    "    min_depth\n",
    "        Minimum depth in km for search\n",
    "    max_depth\n",
    "        Maximum depth in km for search\n",
    "    min_magnitude\n",
    "        Minimum magnitude for search\n",
    "    max_magnitude\n",
    "        Maximum magnitude for search\n",
    "    start_time\n",
    "        Start date and time for search\n",
    "    end_time\n",
    "        End date and time for search\n",
    "        \n",
    "    Returns\n",
    "    -------\n",
    "    pandas.DateFrame of resulting events\n",
    "    \"\"\"\n",
    "    # Convert start_time and end_time to strings\n",
    "    start_time = start_time.strftime(\"%Y-%m-%dT%H:%M:%S\")\n",
    "    end_time = end_time.strftime(\"%Y-%m-%dT%H:%M:%S\")\n",
    "    # Use the more efficient f-string formatting\n",
    "    query_string = (\n",
    "        \"https://quakesearch.geonet.org.nz/csv?bbox=\"\n",
    "        f\"{min_longitude},{min_latitude},{max_longitude},\"\n",
    "        f\"{max_latitude}&minmag={min_magnitude}\"\n",
    "        f\"&maxmag={max_magnitude}&mindepth={min_depth}\"\n",
    "        f\"&maxdepth={max_depth}&startdate={start_time}\"\n",
    "        f\"&enddate={end_time}\")\n",
    "    print(f\"Using query: {query_string}\")\n",
    "    response = requests.get(query_string)\n",
    "    with open(\"data/earthquakes.csv\", \"wb\") as f:\n",
    "        f.write(response.content)\n",
    "    earthquakes = pd.read_csv(\n",
    "        \"data/earthquakes.csv\", \n",
    "        parse_dates=[\"origintime\", \"modificationtime\"],\n",
    "        dtype={\"publicid\": str})\n",
    "    earthquakes = earthquakes.rename(\n",
    "        columns={\" magnitude\": \"magnitude\",\n",
    "                 \" latitude\": \"latitude\",\n",
    "                 \" depth\": \"depth\"})\n",
    "    return earthquakes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets quickly run this function to get the data.  There won't be any output.  Note that I didn't\n",
    "provide these data in the repository because:\n",
    "1. I don't have permission to re-distribute the data and,\n",
    "2. this dataset gets updated frequently!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using query: https://quakesearch.geonet.org.nz/csv?bbox=164.0,-49.0,182.0,-40.0&minmag=3&maxmag=9.0&mindepth=0.0&maxdepth=500.0&startdate=2015-01-01T00:00:00&enddate=2020-01-01T00:00:00\n",
      "Index(['publicid', 'eventtype', 'origintime', 'modificationtime', 'longitude',\n",
      "       'latitude', 'magnitude', 'depth', 'magnitudetype', 'depthtype',\n",
      "       'evaluationmethod', 'evaluationstatus', 'evaluationmode', 'earthmodel',\n",
      "       'usedphasecount', 'usedstationcount', 'magnitudestationcount',\n",
      "       'minimumdistance', 'azimuthalgap', 'originerror',\n",
      "       'magnitudeuncertainty'],\n",
      "      dtype='object')\n"
     ]
    }
   ],
   "source": [
    "earthquakes = get_geonet_quakes(\n",
    "    start_time=datetime.datetime(2015, 1, 1),min_magnitude=3)\n",
    "\n",
    "print(earthquakes.columns)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Do something useful: applying functions to a dataframe\n",
    "\n",
    "We saw in the previous notebook how we could sort and slice a dataframe in various ways.Let's look at\n",
    "an example of doing some actual calculations with a dataframe.  We will use some of what we have\n",
    "learnt to calculate the occurance rate of earthquakes within a region.  In this case\n",
    "we will take the region around the top of South Island, containing the faults that ruptured\n",
    "in the Kaikoura earthquake."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using query: https://quakesearch.geonet.org.nz/csv?bbox=172.37,-43.12,174.95,-41.15&minmag=2&maxmag=9.0&mindepth=0.0&maxdepth=500.0&startdate=2010-01-01T00:00:00&enddate=2020-01-01T00:00:00\n",
      "Downloaded 30854 earthquakes\n"
     ]
    }
   ],
   "source": [
    "kaikoura = get_geonet_quakes(\n",
    "    min_latitude=-43.12, max_latitude=-41.15,\n",
    "    min_longitude=172.37, max_longitude=174.95,\n",
    "    start_time=datetime.datetime(2010, 1, 1),\n",
    "    min_magnitude=2)\n",
    "# Rename those columns\n",
    "kaikoura = kaikoura.rename(\n",
    "    columns={\" magnitude\": \"magnitude\",\n",
    "             \" latitude\": \"latitude\",\n",
    "             \" depth\": \"depth\"})\n",
    "\n",
    "print(f\"Downloaded {len(kaikoura)} earthquakes\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Earthquake rate is the number of earthquakes per unit time.  For a dataset like this we can calculate\n",
    "rate as 1 over the inter-event time.\n",
    "\n",
    "First we need to sort by origin time:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "kaikoura = kaikoura.sort_values(by=[\"origintime\"], ignore_index=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we need to calculate the time between each successive earthquake.  We can do this by\n",
    "taking the `origintime` column away from a one-sample shifted version of the `origintime` column,\n",
    "much like we did for calculating the temperature gradient in the DFDP drillhole. In Pandas we\n",
    "can do this quickly using the `.diff` method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0                          NaT\n",
      "1       0 days 04:49:17.548000\n",
      "2       0 days 05:37:57.876000\n",
      "3       0 days 06:41:39.160000\n",
      "4       0 days 11:29:03.363000\n",
      "                 ...          \n",
      "30849   0 days 04:00:38.425000\n",
      "30850   0 days 03:35:27.277000\n",
      "30851   0 days 10:19:42.972000\n",
      "30852   0 days 10:59:17.003000\n",
      "30853   1 days 00:24:15.080000\n",
      "Name: origintime, Length: 30854, dtype: timedelta64[ns]\n"
     ]
    }
   ],
   "source": [
    "inter_event_time = kaikoura[\"origintime\"].diff()\n",
    "print(inter_event_time)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This column is in `timedelta` format. Lets convert it to seconds:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0              NaN\n",
      "1        17357.548\n",
      "2        20277.876\n",
      "3        24099.160\n",
      "4        41343.363\n",
      "           ...    \n",
      "30849    14438.425\n",
      "30850    12927.277\n",
      "30851    37182.972\n",
      "30852    39557.003\n",
      "30853    87855.080\n",
      "Name: origintime, Length: 30854, dtype: float64\n"
     ]
    }
   ],
   "source": [
    "inter_event_time = inter_event_time.dt.total_seconds()\n",
    "print(inter_event_time)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The rate is simply 1 / `inter_event_time`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "rate = 1 / inter_event_time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0             NaN\n",
      "1        0.000058\n",
      "2        0.000049\n",
      "3        0.000041\n",
      "4        0.000024\n",
      "           ...   \n",
      "30849    0.000069\n",
      "30850    0.000077\n",
      "30851    0.000027\n",
      "30852    0.000025\n",
      "30853    0.000011\n",
      "Name: origintime, Length: 30854, dtype: float64\n"
     ]
    }
   ],
   "source": [
    "print(rate)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now put this column back into the dataframe and call it `\"rate\"`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "kaikoura = kaikoura.merge(rate.rename(\"rate\"), left_index=True, right_index=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can plot rate as a function of time. Remember that the time values we use here are the\n",
    "event origin-times, which do not directly represent the time of the rate calculated, which\n",
    "is an average between events:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "c35c14bc93494aa9bc49b12d915f61ac",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "FigureCanvasNbAgg()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ax = kaikoura.plot(x=\"origintime\", y=\"rate\")\n",
    "ax.set_ylabel(\"Instantaneous Rate (earthquakes per second)\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is lots more that you can do with dataframes.  Pandas is in heavy use in the datascience\n",
    "world, and allows you to quickly explore datasets.  Hopefully you will find it useful at some\n",
    "point in your Geoscience career.\n",
    "\n",
    "As a final thing, you can save your dataframes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "kaikoura.to_csv(\"data/kaikoura.csv\", index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Summary\n",
    "\n",
    "That covers *some* of the basics of `dataframes` in pandas.  You should be able to use them to replace most of what you would have done with spreadsheets to allow you to work in a more programatic and reproducible way. We also demonstrated some of the basics of fitting simple polynomials to data.  Note that numpy's `polyfit` function isn't limited to 1st order polynomials, the sky is the limit!\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Extension: more sophisticated analysis of a dataframe\n",
    "\n",
    "### Extension: computing b-values\n",
    "\n",
    "Earthquake generally follow a Gutenberg-Richter relationship, where the logarithm of the cumulative number of earthquakes above a given magnitude is proportional to the magnitude:\n",
    "\\begin{equation}\n",
    "    \\log_{10}{N} = a - bM\n",
    "\\end{equation}\n",
    "where *M* is magnitude, *N* is the number of events with magnitude >= *M*, and *a* and *b* are constants. This is a nice simple straight-line equation with offset from the origin given by *a* and the gradient by *b*.\n",
    "\n",
    "Some studies (for example, [Nuannin et al., 2005](https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2005GL022679)) have found variations in b-value with time and space, and related this to changes in stress.  Lets see if we can:\n",
    "1. Calculate the b-value for our dataset;\n",
    "2. Do some sliding-window fu to get at b-value variations in time.\n",
    "\n",
    "To kick us off, note that in our analysis we are going to miss one fundamental thing which means that everything we do is wrong.  That thing is catalogue completeness, upon which our b-value calculations depend. To show that completeness, and have a first pass at computing b-values, lets look at a cumulative distribution of earthquake magnitudes.\n",
    "\n",
    "### Plotting cumulative distributions\n",
    "\n",
    "We want an inverse cumulative plot of magnitudes. We can do this with matplotlib's `hist` by setting the `cumulative` argument to `-1`, and the `density` argument set to `True`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "73475e0bddbb4120bc6d2821bb624d4e",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "FigureCanvasNbAgg()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.hist(\n",
    "    kaikoura[\"magnitude\"], bins=len(kaikoura), \n",
    "    histtype=\"step\", density=True, log=True, \n",
    "    cumulative=-1)\n",
    "ax.set_xlabel(\"Magnitude\")\n",
    "ax.set_ylabel(\"Cumulative density\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looks pretty straight on a log-normal plot, as we would expect from the Gutenberg-Richter law.  However, somewhere between M 2 and 3 it stops being straight.  We assume that our catalogue completeness is somewhere in here.  This means that we think that, if we could detect and catalogue all the earthquakes all the way down to the tiny earthquakes, we would continue seeing this log-normal relationship. So, we assume that below our magnitude of completeness ($M_C$) we are missing earthquakes.  This seems reasonable, as earthquakes get smaller they get much harder to detect simply because their amplitudes are greatly reduced.\n",
    "\n",
    "Lets *assume* our catalogue is complete to $M_C=2.5$ and try and fit a straight line to our cumulative-density plot.\n",
    "\n",
    "First we will count how many times each magnitude appears in our dataset, we will use a handy object in Pythons native `collections` library, called `Counter`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[(2.398, 9), (2.508, 9), (2.23, 9), (2.532, 8), (2.549, 8), (2.459, 7), (2.128, 7), (2.293, 7), (2.147, 7), (2.59, 7)]\n"
     ]
    }
   ],
   "source": [
    "from collections import Counter\n",
    "\n",
    "counted_magnitudes = Counter(kaikoura[\"magnitude\"])\n",
    "\n",
    "print(counted_magnitudes.most_common(10))  # Print the most common 10 magnitudes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cool, that gives us a list of the magnitude and the number of occurrences of that magnitude.  What we actually want is magnitudes and the number of occurrences of that magnitude and *any magnitude above that magnitude*.  To do that we will:\n",
    "1. create a unique set of all the magnitudes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "magnitudes = set(kaikoura[\"magnitude\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. Make a sorted `list` from this set:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "magnitudes = sorted(list(magnitudes), reverse=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3. Remove all magnitudes below our completeness using a [list comprehension](https://docs.python.org/3/tutorial/datastructures.html#list-comprehensions) (because they are handy):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "magnitudes = [m for m in magnitudes if m >= 2.5]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "4. Initialise an empty array in which we will put the cumulative density function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "density = np.zeros(len(magnitudes))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "5. Loop through the magnitudes from largest to smallest and add the number of occurrences of that magnitude to the total occurrences of the previous magnitude bin:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "density[0] = counted_magnitudes[magnitudes[0]]\n",
    "for i, magnitude in enumerate(magnitudes[1:]):\n",
    "    density[i + 1] = density[i] + counted_magnitudes[magnitude]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets check that that looks okay by plotting it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "5cf6146075b64e6b928b84220ee61a85",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "FigureCanvasNbAgg()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.semilogy(magnitudes, density)\n",
    "ax.set_ylabel(\"Cumulative density\")\n",
    "ax.set_xlabel(\"Magnitude\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looks good! Now lets try and fit a line to it.  We can use `numpy`'s solvers to do this. Because this is a nice\n",
    "simple equation we will use the [numpy.polyfit](https://docs.scipy.org/doc/numpy/reference/generated/numpy.polyfit.html) function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "a=6.27, b=-0.85\n"
     ]
    }
   ],
   "source": [
    "coefficients, residual, rank, singular_values, rcondition = np.polyfit(\n",
    "    magnitudes, np.log10(density), deg=1, full=True)\n",
    "b, a = coefficients\n",
    "print(f\"a={a:.2f}, b={b:.2f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "b is usually close to 1 (note that the gradient calculated above is negative, which is already taken care of in the Gutenberg-Richter law). \n",
    "\n",
    "Lets estimate the density from our calculated values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "# To make our lives easier we will convert our magnitudes to a numpy array:\n",
    "magnitudes = np.array(magnitudes)\n",
    "estimated_density = 10 ** (a + (magnitudes * b))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Right, lets see if it fits!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "32cad2c506f44f6da0e23def62d25d74",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "FigureCanvasNbAgg()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.semilogy(\n",
    "    magnitudes, density, marker=\"+\", linestyle=\"None\",\n",
    "    label=\"Data\")\n",
    "ax.semilogy(magnitudes, estimated_density, label=\"Model\")\n",
    "ax.set_ylabel(\"Cumulative density\")\n",
    "ax.set_xlabel(\"Magnitude\")\n",
    "ax.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that because we specified `full=True` in our call to `polyfit`, we were returned a range of metrics about how well-fitted our data were.  The easiest one of those to understand is the residual:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[7.36381152]\n"
     ]
    }
   ],
   "source": [
    "print(residual)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is a measure of the misfit between our model and our data.\n",
    "\n",
    "Lets build a simple function to do this with the aim of applying this to distinct time-chunks of our dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calc_b_value(magnitudes, completeness_magnitude=2.5):\n",
    "    \"\"\"\n",
    "    Calculate the b-value for a range of magnitudes.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    magnitudes\n",
    "        List or array of magnitudes\n",
    "    completeness_magnitude\n",
    "        Magnitude of completeness for the dataset\n",
    "        \n",
    "    Returns\n",
    "    -------\n",
    "    b-value\n",
    "    \"\"\"\n",
    "    counted_magnitudes = Counter(magnitudes)\n",
    "    magnitudes = sorted(list(set(magnitudes)), reverse=True)\n",
    "    magnitudes = np.array(magnitudes)\n",
    "    # Remove magnitudes less than completess\n",
    "    magnitudes = magnitudes[magnitudes >= completeness_magnitude]\n",
    "    # Calculate density\n",
    "    density = np.zeros(len(magnitudes))\n",
    "    density[0] = counted_magnitudes[magnitudes[0]]\n",
    "    for i, magnitude in enumerate(magnitudes[1:]):\n",
    "        density[i + 1] = density[i] + counted_magnitudes[magnitude]\n",
    "    coefficients, residual, rank, singular_values, rcondition = np.polyfit(\n",
    "        magnitudes, np.log10(density), deg=1, full=True)\n",
    "    b, a = coefficients\n",
    "    return b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets check that we get the same b-value as we did before:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "b=-0.85\n"
     ]
    }
   ],
   "source": [
    "b = calc_b_value(kaikoura[\"magnitude\"])\n",
    "print(f\"b={b:.2f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Rolling windows with Pandas\n",
    "\n",
    "Pandas has neat ways of doing rolling windows.  We will use this to do two things:\n",
    "1. Calculate the median date for every 2000 earthquakes;\n",
    "2. Calculate the b-value for every 2000 earthquakes.\n",
    "\n",
    "We will then plot these and see if we see any variations.\n",
    "\n",
    "To calculate the median date we will:\n",
    "1. sort the dataframe by `\"origintime\"`\n",
    "2. Extract just the `\"origintime\"` and `\"magnitude\"` columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "window_size = 2000\n",
    "\n",
    "kaikoura = kaikoura.sort_values(by=[\"origintime\"], ignore_index=True)\n",
    "magnitude_times = pd.concat([kaikoura[\"origintime\"], kaikoura[\"magnitude\"]], axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3. Make a new column containing the seconds since the first event - pandas doesn't have a simple way to calculate the median of a range of datetimes, so we will change to working in seconds since a reference time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "seconds_offset = (magnitude_times.origintime - magnitude_times.origintime[0]).dt.total_seconds()\n",
    "magnitude_times = magnitude_times.merge(\n",
    "    seconds_offset.rename(\"seconds_offset\"), left_index=True, right_index=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "4. Compute the rolling median of the seconds_offset column"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "window_median = magnitude_times.seconds_offset.rolling(window_size).median()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "5. Convert this column to timedelta objects"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "window_median = pd.to_timedelta(window_median, unit=\"S\") # Unit is seconds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "6. Add the reference time to these to get back to real-time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "window_median += magnitude_times.origintime[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "7. Put this into the dataframe as a new column"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [],
   "source": [
    "magnitude_times = magnitude_times.merge(\n",
    "    window_median.rename(\"window_median\"), left_index=True, right_index=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Computing moving window b-values\n",
    "\n",
    "Computing the moving b-value is a little simpler to write, but slower to run.  We will use the function we wrote above and pandas `.rolling().apply(func)` chained method to apply our custom `func` to our column:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "b_values = magnitude_times.magnitude.rolling(window_size).apply(calc_b_value)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets quickly convert those from gradients to b-values by multiplying by -1:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [],
   "source": [
    "b_values *= -1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can put those back into the dataframe:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [],
   "source": [
    "magnitude_times = magnitude_times.merge(\n",
    "    b_values.rename(\"b_value\"), left_index=True, right_index=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plotting the results\n",
    "\n",
    "Now lets plot it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "7b8601034f0f46ab9d5e135cab0056ca",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "FigureCanvasNbAgg()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ax = magnitude_times.plot(x=\"window_median\", y=\"b_value\", kind=\"scatter\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### What next?\n",
    "\n",
    "There are some pretty impressive variations there! In particular there are strong variations in 2013 and 2016, right around when the Cook Strait and Kaikoura earthquakes happened. I wonder if there is anything in that...? **before we get ahead of ourselves**, we missed some key things here that mean that this result is not interpretable:\n",
    "1. Not all magnitudes are equal, and we were just using GeoNet's summary magnitude;\n",
    "2. We fixed the magnitude of completeness when in reality completeness depends on a range of factors and is time-varying;\n",
    "3. We haven't taken spatial variations into account - we have looked at quite a large region here.\n",
    "\n",
    "We could get around those factors though and extend our rolling window to compute completeness alongside b-value. Potential student project...?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** Using pandas rolling windows, find the mean earthquake location for every window we used above. You will need to compute the rolling mean for latitude, longitude and depth.  Make three plots to show how latitude, longitude and depth vary with time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your answer here"
   ]
  }
 ],
 "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.9"
  },
  "toc-autonumbering": false,
  "toc-showmarkdowntxt": false,
  "toc-showtags": false
 },
 "nbformat": 4,
 "nbformat_minor": 4
}