{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1><center>IAGA Summer School 2019</center></h1>\n",
    "\n",
    "<h1><center>Visualising geomagnetic data: time series, Bartels plots and contour plots</center></h1>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import notebook dependencies\n",
    "\n",
    "import os\n",
    "import sys\n",
    "import datetime as dt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "from matplotlib import cm\n",
    "\n",
    "sys.path.append('..')\n",
    "from src import mag_lib as mag\n",
    "\n",
    "# Set the data paths\n",
    "ESK_2003_PATH = os.path.abspath('../data/external/ESK_2003/')\n",
    "ESK_HOURLY_PATH = os.path.abspath('../data/external/ESK_hourly/')\n",
    "K_INDS_PATH = os.path.abspath('../data/external/k_inds/')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1. Plotting 1-minute and 1-hour observatory means"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load the 2003 1-minute data for ESK and resample to hourly means"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "obs = \"ESK\"\n",
    "# 1-minute data\n",
    "obs_min = mag.load_year(observatory='esk', year=2003, path=ESK_2003_PATH)\n",
    "# 1-hour data\n",
    "obs_hourly = obs_min.resample('1h').mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Print the first few lines of each data set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "obs_min.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "obs_hourly.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### >> USER INPUT HERE: Choose start and end times of the time series plots in the format '2003-mm-dd hh'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "start_date = '2003-01-01 00'\n",
    "end_date = '2003-01-31 23'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the 1-minute data for the selected period"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = obs_min.loc[start_date:end_date, 'ESKX':'ESKZ'].plot(subplots=True, figsize=(10,7),\n",
    "                                                          title='1-minute means: ESK 2003', legend=False)\n",
    "plt.xlabel('Date')\n",
    "ax[0].set_ylabel('X (nT)')\n",
    "ax[1].set_ylabel('Y (nT)')\n",
    "ax[2].set_ylabel('Z (nT)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the hourly mean data for the same period"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = obs_hourly.loc[start_date:end_date, 'ESKX':'ESKZ'].plot(subplots=True, figsize=(10,7),\n",
    "                                                             title='1-hour means: ESK 2003', legend=False)\n",
    "plt.xlabel('Date')\n",
    "ax[0].set_ylabel('X (nT)')\n",
    "ax[1].set_ylabel('Y (nT)')\n",
    "ax[2].set_ylabel('Z (nT)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare hourly means with the annual means"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Resample to annual means\n",
    "obs_annual = obs_min.resample('1Y').mean()\n",
    "dates = obs_hourly.loc[start_date:end_date].reset_index()['DATE_TIME']\n",
    "annual = pd.DataFrame(index=dates, data=[[obs_annual.loc['2003-12-31','ESKX'], obs_annual.loc['2003-12-31','ESKY'],\n",
    "                                       obs_annual.loc['2003-12-31','ESKZ']]], columns=['ESKX', 'ESKY', 'ESKZ'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, (ax1, ax2, ax3) = plt.subplots(nrows=3, sharex=True, figsize=(10,9))\n",
    "# X component\n",
    "ax1.plot(dates, obs_hourly.loc[start_date:end_date,'ESKX'], 'k')\n",
    "ax1.plot(annual['ESKX'], 'k')\n",
    "ax1.fill_between(dates, obs_hourly.loc[start_date:end_date,'ESKX'], annual['ESKX'].values,\n",
    "                 where=obs_hourly.loc[start_date:end_date,'ESKX'] < annual['ESKX'], facecolor='lightblue',\n",
    "                 interpolate=True)\n",
    "ax1.fill_between(dates, obs_hourly.loc[start_date:end_date,'ESKX'], annual['ESKX'],\n",
    "                 where=obs_hourly.loc[start_date:end_date,'ESKX'] >= annual['ESKX'], facecolor='pink',\n",
    "                 interpolate=True)\n",
    "ax1.set_ylabel('X (nT)')\n",
    "# Y component\n",
    "ax2.plot(dates, obs_hourly.loc[start_date:end_date,'ESKY'], 'k')\n",
    "ax2.plot(annual['ESKY'], 'k')\n",
    "ax2.fill_between(dates, obs_hourly.loc[start_date:end_date,'ESKY'], annual['ESKY'].values,\n",
    "                where=obs_hourly.loc[start_date:end_date,'ESKY'] < annual['ESKY'], facecolor='lightblue',\n",
    "                 interpolate=True)\n",
    "ax2.fill_between(dates, obs_hourly.loc[start_date:end_date,'ESKY'], annual['ESKY'],\n",
    "                 where=obs_hourly.loc[start_date:end_date,'ESKY'] >= annual['ESKY'], facecolor='pink',\n",
    "                 interpolate=True)\n",
    "ax2.set_ylabel('Y (nT)')\n",
    "# Z component\n",
    "ax3.plot(dates, obs_hourly.loc[start_date:end_date,'ESKZ'], 'k')\n",
    "ax3.plot(annual['ESKZ'], 'k')\n",
    "ax3.fill_between(dates, obs_hourly.loc[start_date:end_date,'ESKZ'], annual['ESKZ'].values,\n",
    "                where=obs_hourly.loc[start_date:end_date,'ESKZ'] < annual['ESKZ'], facecolor='lightblue',\n",
    "                 interpolate=True)\n",
    "ax3.fill_between(dates, obs_hourly.loc[start_date:end_date,'ESKZ'], annual['ESKZ'],\n",
    "                where=obs_hourly.loc[start_date:end_date,'ESKZ'] >= annual['ESKZ'], facecolor='pink',\n",
    "                 interpolate=True)\n",
    "ax3.set_ylabel('Z (nT)')\n",
    "plt.xticks(rotation = 45)\n",
    "plt.xlabel('Date', )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. Daily, seasonal and solar variations in declination\n",
    "\n",
    "Here we plot the daily variation in declination at Eskdalemuir (with the secular variation removed) over a selectable number of years (files from 1911 to 2017 are available in this repository). The seasonal variation in the amplitude of the daily variation is well illustrated by the plot and, if you extend the plot over a complete solar cycle, it is clear that there is a solar cycle modulation as well. Try 1986-1997 as an example (see also https://en.wikipedia.org/wiki/List_of_solar_cycles for a list of suitable date ranges.)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### >> USER INPUT HERE: Choose the start and end years of the contour plot.\n",
    "\n",
    "#### Files from 1911 to 2017 are avaliable in this repository. We recommend plotting at least one 11 year solar cycle ( e.g. 1986 to 1997) in order to see all the variations. The plotting may be a little slow if you plot the entire data set - be patient!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "start_year = 1986\n",
    "end_year = 1997"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read declination data from IAGA-2002 format files, or calculate it from the given components if necessary\n",
    "hourly = mag.read_obs_hmv_declination(obscode='esk', year_st=start_year, year_fn=end_year, folder=ESK_HOURLY_PATH)\n",
    "hourly.reset_index(inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate the time as minutes since midnight and date as Julian date\n",
    "# This is needed for plotting purposes as the axes must be numbers, not datetime objects or strings\n",
    "hourly['mins_count'] = hourly['DATE_TIME'].map(lambda x: x.time().hour*60 + x.time().minute)\n",
    "hourly['julian_date'] = hourly['DATE_TIME'].map(lambda x: int(x.to_julian_date()))\n",
    "\n",
    "# Calculate the monthly mean for each hourly interval of the day\n",
    "monthly = hourly.groupby([hourly['DATE_TIME'].dt.year, hourly['DATE_TIME'].dt.month, hourly['DATE_TIME'].dt.hour]).mean()\n",
    "\n",
    "# Calculate the monthly mean over all hourly intervals\n",
    "monthly_all = hourly.groupby([hourly['DATE_TIME'].dt.year, hourly['DATE_TIME'].dt.month]).mean()\n",
    "\n",
    "# Calculate the daily declination variations as the monthly average of the hourly intervals minus the total monthly mean\n",
    "variations = monthly['D'].values - monthly_all['D'].values.repeat(24)\n",
    "\n",
    "# New variables for the dates and times\n",
    "times = monthly.mins_count.values\n",
    "days = monthly.julian_date.values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convert the Julian dates back to datetimes and strings for axis labels \n",
    "epoch = pd.to_datetime(0, unit='s').to_julian_date()\n",
    "dates = pd.to_datetime(days[0:-1:24]-epoch, unit='D').strftime('%Y')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot declination variations using trisurf. If you used '%matplotlib notebook' above, the plot will be interactive - you\n",
    "can pan and zoom for different perspectives of the 3D surface.\n",
    "\n",
    "The spacing of the tick labels on the date axis can be controlled by changing *n* below, e.g. n=5 gives one tick per 5 years."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(10,10))\n",
    "ax = fig.add_subplot(1, 1, 1, projection='3d')\n",
    "ax.plot_trisurf(days, times, variations, cmap=cm.jet, vmin=-0.15, vmax=0.15, antialiased=True)\n",
    "\n",
    "# Set x tick labels to 1 per n year\n",
    "n=1\n",
    "ax.set_xticks(days[0::24*12*n])\n",
    "ax.set_xticklabels(dates[0::12*n],rotation=0)\n",
    "ax.set_xlabel('Year', rotation=60, labelpad=8)\n",
    "ax.set_xlim([days[0], days[-1]])\n",
    "\n",
    "# Set y labels to hours\n",
    "ax.set_yticks(np.arange(0, 1560, 120))\n",
    "ax.set_yticklabels(np.arange(0, 26, 2))\n",
    "ax.set_ylabel('Hour (UT)')\n",
    "\n",
    "# z axis setup\n",
    "ax.set_zlabel('D variations (degrees)', labelpad=8)\n",
    "ax.set_zlim([-0.2, 0.2])\n",
    "ax.tick_params(axis='z', pad=5)\n",
    "\n",
    "# Set initial viewpoint elevation and azimuth\n",
    "ax.view_init(elev=60, azim=210)\n",
    "\n",
    "# Set the background colour to white\n",
    "ax.w_xaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))\n",
    "ax.w_yaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))\n",
    "ax.w_zaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. Plotting observatory data by the Bartels rotation number\n",
    "\n",
    "Solar rotations are numbered by the Bartels solar rotation number, with Day 1 of rotation 1 chosen as 8th February 1832. In this section, hourly mean values for 2003 at Eskdalemuir Observatory are plotted ordered by Bartels rotation number (the X, Y and Z geomagnetic components may be selected below.) The plot shows a number of the features of geomagnetic field behaviour:\n",
    "\n",
    "1.\tThe annual mean is plotted as a horizontal line in each row dividing the plots into sections ‘above’ and ‘below’ the mean. The proportions of the plots above and below changes as the year progresses because of the slow core field changes with time the secular variation.\n",
    "\n",
    "2.\tThe daily variation in each element is clear. Note the differences between winter and summer, which we also saw in the 3D contour plot above.\n",
    "\n",
    "3.\tAlthough substantially attenuated by taking hourly means, times of magnetic disturbances are obvious.\n",
    "\n",
    "4.\tThe rows are plotted 27 days long because the equatorial rotation period, as seen from Earth, is approximately 27 days. As a result, if a region on the Sun responsible for a disturbance on one day survives a full rotation, it may cause a further disturbance 27 days later and this will line up vertically in the plots.\n",
    "\n",
    "\n",
    "This plot reproduces an example found in the Eskdalemuir monthly bulletin of December 2003, which can be found at http://geomag.bgs.ac.uk/data_service/data/bulletins/esk/2003/esk_dec03.pdf (pages 22-24)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#These functions are used to produce the Bartels plots\n",
    "def bartels_rotation(datetime, return_indexes=True):\n",
    "    \"\"\"\n",
    "    \n",
    "    Args:\n",
    "        date (datetime/DatetimeIndex)\n",
    "        \n",
    "    Returns:\n",
    "        tuple (rotation number, day within rotation, hour within day)\n",
    "        \n",
    "    \"\"\"\n",
    "    if type(datetime) is pd.DatetimeIndex:\n",
    "        date, hour = datetime.date, datetime.hour\n",
    "    elif type(datetime) is dt.datetime:\n",
    "        date, hour = datetime.date(), datetime.hour\n",
    "    # Number of days since Bartels 0\n",
    "    ndays = pd.to_timedelta(date - dt.date(1832, 2, 8)).days\n",
    "    bartels_rotation = ndays//27 + 1\n",
    "    day = ndays%27 + 1\n",
    "    if return_indexes:\n",
    "        return bartels_rotation, day, hour\n",
    "    else:\n",
    "        return tuple([item.values for item in (bartels_rotation, day, hour)])\n",
    "\n",
    "\n",
    "def create_bartels_index(df, nlevels=3):\n",
    "    \"\"\"Replaces a DatetimeIndex with a Bartels MultiIndex\n",
    "    \n",
    "    If nlevels is 3, the MultiIndex has levels:\n",
    "        0: (bartrot) Bartel rotation number\n",
    "        1: (bartrotday) Day within the current Bartel rotation\n",
    "        2: (hourofday) Hour within the current day\n",
    "    else if nlevels is 2:\n",
    "        0: (bartrot) Bartel rotation number\n",
    "        1: (bartrotday) Fractional day\n",
    "    \n",
    "    Args:\n",
    "        df (DataFrame)\n",
    "    \n",
    "    Returns:\n",
    "        DataFrame\n",
    "\n",
    "    \"\"\"\n",
    "    bartrot, bartrotday, hourofday = bartels_rotation(df.index)\n",
    "    if nlevels == 3:\n",
    "        df.index = pd.MultiIndex.from_arrays(\n",
    "            (bartrot, bartrotday, hourofday), names=(\"bartrot\", \"bartrotday\", \"hourofday\")\n",
    "        )\n",
    "    elif nlevels == 2:\n",
    "        df.index = pd.MultiIndex.from_arrays(\n",
    "            (bartrot, bartrotday + hourofday/24), names=(\"bartrot\", \"bartrotday\")\n",
    "        )\n",
    "    else:\n",
    "        raise NotImplementedError()\n",
    "    return df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate annual means and append them to the hourly dataframe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "obs = 'ESK'\n",
    "annual_means = (\n",
    "    obs_min[[f\"{obs}{i}\" for i in \"XYZF\"]].resample(\"1Y\").mean()\n",
    "    .rename(columns={f\"{obs}{i}\":f\"{obs}{i}_mean\" for i in \"XYZF\"})\n",
    "    .reindex(index=obs_hourly.index, method=\"nearest\")\n",
    ")\n",
    "obs_hourly = obs_hourly.join(annual_means)\n",
    "obs_hourly.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Change the dataframe index to a Bartels-rotation-based MultiIndex and preserve the time index as a column instead"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "obs_hourly[\"time\"] = obs_hourly.index\n",
    "obs_hourly = create_bartels_index(obs_hourly, nlevels=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Identify the calendar month starts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "month_starts = obs_hourly[\"time\"].where(\n",
    "    (obs_hourly[\"time\"].dt.day == 1) & (obs_hourly[\"time\"].dt.hour == 0)\n",
    ").dropna()\n",
    "month_starts"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### >> USER INPUT HERE: Choose the component to plot\n",
    "**var** options: <br>\n",
    "X <br>\n",
    "Y <br>\n",
    "Z <br>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "var = \"Y\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Loop through each bartrot and plot it on its own axis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "bartrots = range(obs_hourly.index[0][0], obs_hourly.index[-1][0] + 1)\n",
    "\n",
    "fig, axes = plt.subplots(\n",
    "    nrows=len(bartrots), ncols=1, figsize=(10, 10),\n",
    "    gridspec_kw = {'hspace':0},\n",
    "    sharex=True, sharey=True\n",
    ")\n",
    "for bartrot, ax in zip(bartrots, axes):\n",
    "    x = obs_hourly.loc[bartrot].index\n",
    "    y0 = obs_hourly.loc[bartrot][f\"{obs}{var}_mean\"]\n",
    "    y1 = obs_hourly.loc[bartrot][f\"{obs}{var}\"]\n",
    "    ax.plot(x, y0, color=\"black\", linewidth=0.4)\n",
    "    ax.plot(x, y1, color=\"black\", linewidth=0.8, clip_on=False)\n",
    "    ax.fill_between(\n",
    "        x, y0, y1, where=y1 < y0,\n",
    "        facecolor='lightblue', interpolate=True, zorder=9#, clip_on=False currently causes a bug\n",
    "    )\n",
    "    ax.fill_between(\n",
    "        x, y0, y1, where=y1 >= y0,\n",
    "        facecolor='pink', interpolate=True, zorder=10#, clip_on=False\n",
    "    )\n",
    "    ax_r = ax.twinx()\n",
    "    ax_r.set_ylabel(bartrot, fontsize=10)\n",
    "    ax_r.set_yticks([])\n",
    "    # some magic which enables lines from one axis to show on top of other axes\n",
    "    ax.patch.set_visible(False)\n",
    "    \n",
    "    # Add text identifying the start of each month\n",
    "    if bartrot in month_starts.keys():\n",
    "        bartrotday = float(month_starts.loc[bartrot].index.values)\n",
    "        month = month_starts.loc[bartrot].dt.strftime(\"%b\").values[0]\n",
    "        ax.text(bartrotday, y0.iloc[0] - 85, month, verticalalignment=\"top\")\n",
    "\n",
    "ax.set_ylim((y0.iloc[0] - 75, y0.iloc[0] + 75))\n",
    "ax.set_xlabel(\"Day within Bartels rotation\");\n",
    "axes[0].set_title(f\"{obs}-{var} (nT, left axis), by Bartels rotation number (right axis)\",\n",
    "                  fontsize=15);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4. Geomagnetic jerks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Geomagnetic jerks are rapid changes in the trend of secular variation (SV), traditionally thought of as a 'V' (or inverted 'V') shape punctuating several years or decades of roughly linear SV. They are of internal origin, possibly linked to hydromagnetic wave motions in Earth's fluid outer core, and occur at irregular intervals on average about once per decade. Their magnitudes vary according to location, with some events observed globally (e.g. 1969) and others confined to regional scales (e.g. 2003). In this section, we plot SV calculated as first differences of observatory annual means to see various geomagnetic jerks. \n",
    "\n",
    "The code below plots the $X$, $Y$ and $Z$ components of SV for a user-specified observatory; we have provided a file containing the observatory annual means held by BGS, with baseline corrections already applied. These are the data available at http://www.geomag.bgs.ac.uk/data_service/data/annual_means.shtml - that webpage contains a list of observatories for which BGS holds annual mean values, which you may find helpful when selecting an observatory below. To get started, try looking at the $Y$ component of a European observatory around 1969 (e.g. Eskdalemuir \\[ESK\\], Niemegk \\[NGK\\] or Chambon-la-Foret \\[CLF\\])."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### >> USER INPUT HERE: Choose the magnetic observatory (three digit IAGA code)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "obs_name = 'ESK'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Read in the observatory annual means"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "obs_annual = mag.read_obs_ann_mean(obscode=obs_name, filename='../data/external/oamjumpsapplied.dat')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "obs_annual.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate the SV"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "obs_sv = obs_annual[['X', 'Y', 'Z']].astype('float').diff()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "obs_sv.insert(0, 'year', obs_annual['year'].astype('float') - 0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "obs_sv.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot all three components of SV"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "ax = obs_sv.plot(x='year', subplots=True, figsize=(10,7), title='Secular variation at %s' % obs_name,\n",
    "                 legend=False, marker='x')\n",
    "plt.xlabel('Year')\n",
    "ax[0].set_ylabel('dX/dt (nT/yr)')\n",
    "ax[1].set_ylabel('dY/dt (nT/yr)')\n",
    "ax[2].set_ylabel('dZ/dt (nT/yr)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "Plot the secular acceleration (SA, the second time derivative of the gemagnetic field. This is the first time derivative of the SV) at your chosen observatory. What do jerks look like in the SA? **Hint:** This is a very similar calculation to that used above to obtain the SV from field values..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Geomagnetic pulsations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Geomagnetic pulsations are ultra-low frequency (ULF) plasma waves in Earth's magnetosphere. They cover the range from approximately 1 mHz to 1Hz. These magnetic field line resonances are of external origin and are considered a useful tool for remotely diagnosing activity occurring in the plasmasphere, for example, density variations in plasma that can disrupt satellite communications. For interest, we have included a directory called `pulsations` in this repository. The directory contains 1-sec resolution data at Hartland (HAD) for two days on which pulsation activity was observed: 14-12-2018 and 11-06-2019, along with figures showing the pulsations on these days. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Challenge\n",
    "Can you load in the one second magnetic data and filter them to obtain similar plots of pulsation activity?"
   ]
  },
  {
   "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}