IAGA Summer School 2019

Visualising geomagnetic data: time series, Bartels plots and contour plots

In [ ]:
%matplotlib inline
In [ ]:
# Import notebook dependencies

import os
import sys
import datetime as dt
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm

sys.path.append('..')
from src import mag_lib as mag

# Set the data paths
ESK_2003_PATH = os.path.abspath('../data/external/ESK_2003/')
ESK_HOURLY_PATH = os.path.abspath('../data/external/ESK_hourly/')
K_INDS_PATH = os.path.abspath('../data/external/k_inds/')

1. Plotting 1-minute and 1-hour observatory means

Load the 2003 1-minute data for ESK and resample to hourly means

In [ ]:
obs = "ESK"
# 1-minute data
obs_min = mag.load_year(observatory='esk', year=2003, path=ESK_2003_PATH)
# 1-hour data
obs_hourly = obs_min.resample('1h').mean()

Print the first few lines of each data set

In [ ]:
obs_min.head()
In [ ]:
obs_hourly.head()

>> USER INPUT HERE: Choose start and end times of the time series plots in the format '2003-mm-dd hh'

In [ ]:
start_date = '2003-01-01 00'
end_date = '2003-01-31 23'

Plot the 1-minute data for the selected period

In [ ]:
ax = obs_min.loc[start_date:end_date, 'ESKX':'ESKZ'].plot(subplots=True, figsize=(10,7),
                                                          title='1-minute means: ESK 2003', legend=False)
plt.xlabel('Date')
ax[0].set_ylabel('X (nT)')
ax[1].set_ylabel('Y (nT)')
ax[2].set_ylabel('Z (nT)')

Plot the hourly mean data for the same period

In [ ]:
ax = obs_hourly.loc[start_date:end_date, 'ESKX':'ESKZ'].plot(subplots=True, figsize=(10,7),
                                                             title='1-hour means: ESK 2003', legend=False)
plt.xlabel('Date')
ax[0].set_ylabel('X (nT)')
ax[1].set_ylabel('Y (nT)')
ax[2].set_ylabel('Z (nT)')

Compare hourly means with the annual means

In [ ]:
# Resample to annual means
obs_annual = obs_min.resample('1Y').mean()
dates = obs_hourly.loc[start_date:end_date].reset_index()['DATE_TIME']
annual = pd.DataFrame(index=dates, data=[[obs_annual.loc['2003-12-31','ESKX'], obs_annual.loc['2003-12-31','ESKY'],
                                       obs_annual.loc['2003-12-31','ESKZ']]], columns=['ESKX', 'ESKY', 'ESKZ'])
In [ ]:
fig, (ax1, ax2, ax3) = plt.subplots(nrows=3, sharex=True, figsize=(10,9))
# X component
ax1.plot(dates, obs_hourly.loc[start_date:end_date,'ESKX'], 'k')
ax1.plot(annual['ESKX'], 'k')
ax1.fill_between(dates, obs_hourly.loc[start_date:end_date,'ESKX'], annual['ESKX'].values,
                 where=obs_hourly.loc[start_date:end_date,'ESKX'] < annual['ESKX'], facecolor='lightblue',
                 interpolate=True)
ax1.fill_between(dates, obs_hourly.loc[start_date:end_date,'ESKX'], annual['ESKX'],
                 where=obs_hourly.loc[start_date:end_date,'ESKX'] >= annual['ESKX'], facecolor='pink',
                 interpolate=True)
ax1.set_ylabel('X (nT)')
# Y component
ax2.plot(dates, obs_hourly.loc[start_date:end_date,'ESKY'], 'k')
ax2.plot(annual['ESKY'], 'k')
ax2.fill_between(dates, obs_hourly.loc[start_date:end_date,'ESKY'], annual['ESKY'].values,
                where=obs_hourly.loc[start_date:end_date,'ESKY'] < annual['ESKY'], facecolor='lightblue',
                 interpolate=True)
ax2.fill_between(dates, obs_hourly.loc[start_date:end_date,'ESKY'], annual['ESKY'],
                 where=obs_hourly.loc[start_date:end_date,'ESKY'] >= annual['ESKY'], facecolor='pink',
                 interpolate=True)
ax2.set_ylabel('Y (nT)')
# Z component
ax3.plot(dates, obs_hourly.loc[start_date:end_date,'ESKZ'], 'k')
ax3.plot(annual['ESKZ'], 'k')
ax3.fill_between(dates, obs_hourly.loc[start_date:end_date,'ESKZ'], annual['ESKZ'].values,
                where=obs_hourly.loc[start_date:end_date,'ESKZ'] < annual['ESKZ'], facecolor='lightblue',
                 interpolate=True)
ax3.fill_between(dates, obs_hourly.loc[start_date:end_date,'ESKZ'], annual['ESKZ'],
                where=obs_hourly.loc[start_date:end_date,'ESKZ'] >= annual['ESKZ'], facecolor='pink',
                 interpolate=True)
ax3.set_ylabel('Z (nT)')
plt.xticks(rotation = 45)
plt.xlabel('Date', )

2. Daily, seasonal and solar variations in declination

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.)

>> USER INPUT HERE: Choose the start and end years of the contour plot.

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!

In [ ]:
start_year = 1986
end_year = 1997
In [ ]:
# Read declination data from IAGA-2002 format files, or calculate it from the given components if necessary
hourly = mag.read_obs_hmv_declination(obscode='esk', year_st=start_year, year_fn=end_year, folder=ESK_HOURLY_PATH)
hourly.reset_index(inplace=True)
In [ ]:
# Calculate the time as minutes since midnight and date as Julian date
# This is needed for plotting purposes as the axes must be numbers, not datetime objects or strings
hourly['mins_count'] = hourly['DATE_TIME'].map(lambda x: x.time().hour*60 + x.time().minute)
hourly['julian_date'] = hourly['DATE_TIME'].map(lambda x: int(x.to_julian_date()))

# Calculate the monthly mean for each hourly interval of the day
monthly = hourly.groupby([hourly['DATE_TIME'].dt.year, hourly['DATE_TIME'].dt.month, hourly['DATE_TIME'].dt.hour]).mean()

# Calculate the monthly mean over all hourly intervals
monthly_all = hourly.groupby([hourly['DATE_TIME'].dt.year, hourly['DATE_TIME'].dt.month]).mean()

# Calculate the daily declination variations as the monthly average of the hourly intervals minus the total monthly mean
variations = monthly['D'].values - monthly_all['D'].values.repeat(24)

# New variables for the dates and times
times = monthly.mins_count.values
days = monthly.julian_date.values
In [ ]:
# Convert the Julian dates back to datetimes and strings for axis labels 
epoch = pd.to_datetime(0, unit='s').to_julian_date()
dates = pd.to_datetime(days[0:-1:24]-epoch, unit='D').strftime('%Y')

Plot declination variations using trisurf. If you used '%matplotlib notebook' above, the plot will be interactive - you can pan and zoom for different perspectives of the 3D surface.

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.

In [ ]:
fig = plt.figure(figsize=(10,10))
ax = fig.add_subplot(1, 1, 1, projection='3d')
ax.plot_trisurf(days, times, variations, cmap=cm.jet, vmin=-0.15, vmax=0.15, antialiased=True)

# Set x tick labels to 1 per n year
n=1
ax.set_xticks(days[0::24*12*n])
ax.set_xticklabels(dates[0::12*n],rotation=0)
ax.set_xlabel('Year', rotation=60, labelpad=8)
ax.set_xlim([days[0], days[-1]])

# Set y labels to hours
ax.set_yticks(np.arange(0, 1560, 120))
ax.set_yticklabels(np.arange(0, 26, 2))
ax.set_ylabel('Hour (UT)')

# z axis setup
ax.set_zlabel('D variations (degrees)', labelpad=8)
ax.set_zlim([-0.2, 0.2])
ax.tick_params(axis='z', pad=5)

# Set initial viewpoint elevation and azimuth
ax.view_init(elev=60, azim=210)

# Set the background colour to white
ax.w_xaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))
ax.w_yaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))
ax.w_zaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))

3. Plotting observatory data by the Bartels rotation number

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:

  1. The 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.

  2. The daily variation in each element is clear. Note the differences between winter and summer, which we also saw in the 3D contour plot above.

  3. Although substantially attenuated by taking hourly means, times of magnetic disturbances are obvious.

  4. The 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.

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).

In [ ]:
#These functions are used to produce the Bartels plots
def bartels_rotation(datetime, return_indexes=True):
    """
    
    Args:
        date (datetime/DatetimeIndex)
        
    Returns:
        tuple (rotation number, day within rotation, hour within day)
        
    """
    if type(datetime) is pd.DatetimeIndex:
        date, hour = datetime.date, datetime.hour
    elif type(datetime) is dt.datetime:
        date, hour = datetime.date(), datetime.hour
    # Number of days since Bartels 0
    ndays = pd.to_timedelta(date - dt.date(1832, 2, 8)).days
    bartels_rotation = ndays//27 + 1
    day = ndays%27 + 1
    if return_indexes:
        return bartels_rotation, day, hour
    else:
        return tuple([item.values for item in (bartels_rotation, day, hour)])


def create_bartels_index(df, nlevels=3):
    """Replaces a DatetimeIndex with a Bartels MultiIndex
    
    If nlevels is 3, the MultiIndex has levels:
        0: (bartrot) Bartel rotation number
        1: (bartrotday) Day within the current Bartel rotation
        2: (hourofday) Hour within the current day
    else if nlevels is 2:
        0: (bartrot) Bartel rotation number
        1: (bartrotday) Fractional day
    
    Args:
        df (DataFrame)
    
    Returns:
        DataFrame

    """
    bartrot, bartrotday, hourofday = bartels_rotation(df.index)
    if nlevels == 3:
        df.index = pd.MultiIndex.from_arrays(
            (bartrot, bartrotday, hourofday), names=("bartrot", "bartrotday", "hourofday")
        )
    elif nlevels == 2:
        df.index = pd.MultiIndex.from_arrays(
            (bartrot, bartrotday + hourofday/24), names=("bartrot", "bartrotday")
        )
    else:
        raise NotImplementedError()
    return df

Calculate annual means and append them to the hourly dataframe

In [ ]:
obs = 'ESK'
annual_means = (
    obs_min[[f"{obs}{i}" for i in "XYZF"]].resample("1Y").mean()
    .rename(columns={f"{obs}{i}":f"{obs}{i}_mean" for i in "XYZF"})
    .reindex(index=obs_hourly.index, method="nearest")
)
obs_hourly = obs_hourly.join(annual_means)
obs_hourly.head()

Change the dataframe index to a Bartels-rotation-based MultiIndex and preserve the time index as a column instead

In [ ]:
obs_hourly["time"] = obs_hourly.index
obs_hourly = create_bartels_index(obs_hourly, nlevels=2)

Identify the calendar month starts

In [ ]:
month_starts = obs_hourly["time"].where(
    (obs_hourly["time"].dt.day == 1) & (obs_hourly["time"].dt.hour == 0)
).dropna()
month_starts

>> USER INPUT HERE: Choose the component to plot

var options:
X
Y
Z

In [ ]:
var = "Y"

Loop through each bartrot and plot it on its own axis

In [ ]:
bartrots = range(obs_hourly.index[0][0], obs_hourly.index[-1][0] + 1)

fig, axes = plt.subplots(
    nrows=len(bartrots), ncols=1, figsize=(10, 10),
    gridspec_kw = {'hspace':0},
    sharex=True, sharey=True
)
for bartrot, ax in zip(bartrots, axes):
    x = obs_hourly.loc[bartrot].index
    y0 = obs_hourly.loc[bartrot][f"{obs}{var}_mean"]
    y1 = obs_hourly.loc[bartrot][f"{obs}{var}"]
    ax.plot(x, y0, color="black", linewidth=0.4)
    ax.plot(x, y1, color="black", linewidth=0.8, clip_on=False)
    ax.fill_between(
        x, y0, y1, where=y1 < y0,
        facecolor='lightblue', interpolate=True, zorder=9#, clip_on=False currently causes a bug
    )
    ax.fill_between(
        x, y0, y1, where=y1 >= y0,
        facecolor='pink', interpolate=True, zorder=10#, clip_on=False
    )
    ax_r = ax.twinx()
    ax_r.set_ylabel(bartrot, fontsize=10)
    ax_r.set_yticks([])
    # some magic which enables lines from one axis to show on top of other axes
    ax.patch.set_visible(False)
    
    # Add text identifying the start of each month
    if bartrot in month_starts.keys():
        bartrotday = float(month_starts.loc[bartrot].index.values)
        month = month_starts.loc[bartrot].dt.strftime("%b").values[0]
        ax.text(bartrotday, y0.iloc[0] - 85, month, verticalalignment="top")

ax.set_ylim((y0.iloc[0] - 75, y0.iloc[0] + 75))
ax.set_xlabel("Day within Bartels rotation");
axes[0].set_title(f"{obs}-{var} (nT, left axis), by Bartels rotation number (right axis)",
                  fontsize=15);

4. Geomagnetic jerks

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.

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]).

>> USER INPUT HERE: Choose the magnetic observatory (three digit IAGA code)

In [ ]:
obs_name = 'ESK'

Read in the observatory annual means

In [ ]:
obs_annual = mag.read_obs_ann_mean(obscode=obs_name, filename='../data/external/oamjumpsapplied.dat')
In [ ]:
obs_annual.head()

Calculate the SV

In [ ]:
obs_sv = obs_annual[['X', 'Y', 'Z']].astype('float').diff()
In [ ]:
obs_sv.insert(0, 'year', obs_annual['year'].astype('float') - 0.5)
In [ ]:
obs_sv.head()

Plot all three components of SV

In [ ]:
ax = obs_sv.plot(x='year', subplots=True, figsize=(10,7), title='Secular variation at %s' % obs_name,
                 legend=False, marker='x')
plt.xlabel('Year')
ax[0].set_ylabel('dX/dt (nT/yr)')
ax[1].set_ylabel('dY/dt (nT/yr)')
ax[2].set_ylabel('dZ/dt (nT/yr)')

Exercise

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...

Geomagnetic pulsations

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.

Challenge

Can you load in the one second magnetic data and filter them to obtain similar plots of pulsation activity?

In [ ]: