Visualising geomagnetic data

Tip

If you are viewing this on the web, you will not be able to use the interactive components - you will need to be in an active JupyterLab. For this you can use the Swarm Virtual Research Environment service (after first creating a free account). Just hover over the rocket button at the top of this page and select JupyterHub to get started.

# Install latest development version of viresclient
!pip install --upgrade git+https://github.com/ESA-VirES/VirES-Python-Client.git@staging

# %load_ext autoreload
# %autoreload 2

import datetime as dt
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import ipywidgets as widgets
from viresclient import SwarmRequest

## Use to make interactive matplotlib plots
##  - doesn't work smoothly
##  - investigating using plotly/bokeh instead
# %matplotlib widget

%load_ext watermark
%watermark -i -v -p numpy,pandas,xarray,matplotlib,ipywidgets,viresclient
Collecting git+https://github.com/ESA-VirES/VirES-Python-Client.git@staging
  Cloning https://github.com/ESA-VirES/VirES-Python-Client.git (to revision staging) to /tmp/pip-req-build-cnl_214i
  Running command git clone -q https://github.com/ESA-VirES/VirES-Python-Client.git /tmp/pip-req-build-cnl_214i
  Running command git checkout -b staging --track origin/staging
  Switched to a new branch 'staging'
  Branch 'staging' set up to track remote branch 'staging' from 'origin'.
Requirement already satisfied: Jinja2<3.0.0,>=2.10 in /opt/conda/lib/python3.8/site-packages (from viresclient==0.8.1a1) (2.11.3)
Requirement already satisfied: pandas<2.0.0,>=0.18 in /opt/conda/lib/python3.8/site-packages (from viresclient==0.8.1a1) (1.2.3)
Requirement already satisfied: cdflib<=0.3.20,>=0.3.9 in /opt/conda/lib/python3.8/site-packages (from viresclient==0.8.1a1) (0.3.20)
Requirement already satisfied: tables<4.0.0,>=3.4.4 in /opt/conda/lib/python3.8/site-packages (from viresclient==0.8.1a1) (3.6.1)
Requirement already satisfied: tqdm<5.0.0,>=4.23.0 in /opt/conda/lib/python3.8/site-packages (from viresclient==0.8.1a1) (4.60.0)
Requirement already satisfied: xarray<1.0.0,>=0.11.0 in /opt/conda/lib/python3.8/site-packages (from viresclient==0.8.1a1) (0.17.0)
Requirement already satisfied: requests<3.0.0,>=2.0.0 in /opt/conda/lib/python3.8/site-packages (from viresclient==0.8.1a1) (2.25.1)
Requirement already satisfied: attrs>=19.2.0 in /opt/conda/lib/python3.8/site-packages (from cdflib<=0.3.20,>=0.3.9->viresclient==0.8.1a1) (20.3.0)
Requirement already satisfied: numpy in /opt/conda/lib/python3.8/site-packages (from cdflib<=0.3.20,>=0.3.9->viresclient==0.8.1a1) (1.20.2)
Requirement already satisfied: MarkupSafe>=0.23 in /opt/conda/lib/python3.8/site-packages (from Jinja2<3.0.0,>=2.10->viresclient==0.8.1a1) (1.1.1)
Requirement already satisfied: python-dateutil>=2.7.3 in /opt/conda/lib/python3.8/site-packages (from pandas<2.0.0,>=0.18->viresclient==0.8.1a1) (2.8.1)
Requirement already satisfied: pytz>=2017.3 in /opt/conda/lib/python3.8/site-packages (from pandas<2.0.0,>=0.18->viresclient==0.8.1a1) (2021.1)
Requirement already satisfied: six>=1.5 in /opt/conda/lib/python3.8/site-packages (from python-dateutil>=2.7.3->pandas<2.0.0,>=0.18->viresclient==0.8.1a1) (1.15.0)
Requirement already satisfied: certifi>=2017.4.17 in /opt/conda/lib/python3.8/site-packages (from requests<3.0.0,>=2.0.0->viresclient==0.8.1a1) (2020.12.5)
Requirement already satisfied: urllib3<1.27,>=1.21.1 in /opt/conda/lib/python3.8/site-packages (from requests<3.0.0,>=2.0.0->viresclient==0.8.1a1) (1.26.4)
Requirement already satisfied: idna<3,>=2.5 in /opt/conda/lib/python3.8/site-packages (from requests<3.0.0,>=2.0.0->viresclient==0.8.1a1) (2.10)
Requirement already satisfied: chardet<5,>=3.0.2 in /opt/conda/lib/python3.8/site-packages (from requests<3.0.0,>=2.0.0->viresclient==0.8.1a1) (4.0.0)
Requirement already satisfied: numexpr>=2.6.2 in /opt/conda/lib/python3.8/site-packages (from tables<4.0.0,>=3.4.4->viresclient==0.8.1a1) (2.7.3)
Requirement already satisfied: setuptools>=40.4 in /opt/conda/lib/python3.8/site-packages (from xarray<1.0.0,>=0.11.0->viresclient==0.8.1a1) (49.6.0.post20210108)
Building wheels for collected packages: viresclient
  Building wheel for viresclient (setup.py) ... ?25l-
 \
 done
?25h  Created wheel for viresclient: filename=viresclient-0.8.1a1-py3-none-any.whl size=72462 sha256=db65d783534c403376bcc3faad68fc2fc14fa79797bbfe30abf88820048ebac5
  Stored in directory: /tmp/pip-ephem-wheel-cache-ryuxguz8/wheels/fd/03/72/a6e54843a458edd4535b58da3c752ba6444624da09cfa67205
Successfully built viresclient
Installing collected packages: viresclient
  Attempting uninstall: viresclient
    Found existing installation: viresclient 0.8.0a2
    Uninstalling viresclient-0.8.0a2:
      Successfully uninstalled viresclient-0.8.0a2
Successfully installed viresclient-0.8.1a1
Python implementation: CPython
Python version       : 3.8.8
IPython version      : 7.22.0

numpy      : 1.20.2
pandas     : 1.2.3
xarray     : 0.17.0
matplotlib : 3.4.1
ipywidgets : 7.6.3
viresclient: 0.8.1a1

Fetching & loading INTERMAGNET data

INTERMAGNET is the International Real-time Magnetic Observatory Network, a global network of observatories continually monitoring Earth’s magnetic field.

Here we use VirES to access some ground observatory data (see here for more). The following code helps you download the data for a given observatory for a given year, and load them as a Pandas Dataframe.

We will select data from 2003 measured at the Eskdalemuir observatory (ESK) in Scotland, which has a latitude of 55.314° N. If you wish you can select data from a different year or observatory (check this map of INTERMAGNET observatories to find the three-letter code that specifies the observatory).

NB: Data is provided in the NEC (North, East, Centre) geocentrically-defined frame, in contrast to INTERMAGNET data from other sources where the geodetic frame is used. This results in a small rotation in the North (X) and Centre (Z) vectors and a slightly different “latitude” for the observatory location.

def fetch_obs_data_for_years(
    observatory="ESK", year_start=2003, year_end=2003, cadence="M",
    use_xarray=False, **kwargs
):
    """Fetch given years of observatory data from VirES, at minute (M) or hour (H) cadence
    
    Args:
        observatory (str): 3-letter IAGA code
        year_start (int): Chosen year to start
        year_end (int): Year to end with (inclusive)
        cadence (str): "M" for minute, or "H" for hour
        use_xarray (bool): Return xarray.Dataset instead
        
    Returns:
        DataFrame
    """
    # Fetch data from VirES
    request = SwarmRequest()
    request.set_collection(f"SW_OPER_AUX_OBS{cadence}2_:{observatory}", verbose=False)
    request.set_products(measurements=["B_NEC", "IAGA_code"])
    data = request.get_between(
        dt.datetime(year_start, 1, 1),
        dt.datetime(year_end+1, 1, 1),
        **kwargs
    )
    if use_xarray:
        ds = data.as_xarray().drop("Spacecraft")
        return ds
    else:
        # Load data in Pandas Dataframe with X, Y, Z columns
        df = data.as_dataframe(expand=True).drop(columns="Spacecraft")
        df = df.rename(columns={f"B_NEC_{i}": j for i, j in zip("NEC", "XYZ")})
        return df

obs_minute = fetch_obs_data_for_years(
    observatory="ESK",
    year_start=2003,
    year_end=2003,
    cadence="M",
)
obs_minute.head()
Radius Latitude Longitude IAGA_code X Y Z
Timestamp
2003-01-01 00:00:00 6.363967e+06 55.119672 356.8 ESK 17196.514832 -1473.2 46252.152020
2003-01-01 00:01:00 6.363967e+06 55.119672 356.8 ESK 17196.014835 -1473.4 46252.150446
2003-01-01 00:02:00 6.363967e+06 55.119672 356.8 ESK 17196.514832 -1473.4 46252.152020
2003-01-01 00:03:00 6.363967e+06 55.119672 356.8 ESK 17196.715461 -1473.8 46251.952650
2003-01-01 00:04:00 6.363967e+06 55.119672 356.8 ESK 17197.115459 -1473.7 46251.953909

Plotting 1-minute data and their hourly means

This next code cell creates an interactive element that lets you view the data loaded above.

def plot_subset_timeseries(start_date, end_date, hourly_mean=False, show_annual_mean=False, df=obs_minute):
    """Configurably plot a subset of the data
    
    Args:
        start_date (datetime)
        end_date (datetime)
        hourly_mean (bool): Evaluate and plot the hourly means instead
        show_annual_mean (bool): Show offset from annual mean
        df (DataFrame): Assumed to be of form output from fetch_obs_data_for_year
        
    Returns:
        Figure
    """
    # Evaluate annual means to use later,
    #   reindexed with the ending index points for each year
    annual_mean = df.resample("1y").mean()
    annual_mean.index = [df.loc[str(year)].index[-1] for year in df.index.year.unique()]
    # Subset dataframe to selection
    df = df.loc[start_date:end_date]
    # Cut the annual mean (and reindex) to match df
    #  so we can use it directly in ax.fill_between
    annual_mean = annual_mean.loc[str(start_date.year):str(end_date.year)]
    annual_mean = annual_mean.reindex(index=df.index, method="backfill")
    observatory = df["IAGA_code"][0]
    title = f"Minute data from {observatory}"
    if hourly_mean:
        df = df.resample("1h").mean()
        annual_mean = annual_mean.reindex(index=df.index)
        title += ": averaged over each hour"
    if show_annual_mean:
        title += "\nshowing offset from annual mean"
    fig, axes = plt.subplots(nrows=3, figsize=(10, 7), sharex=True)
    for i, cpt in enumerate("XYZ"):
        if show_annual_mean:
            axes[i].fill_between(df.index, df[cpt], annual_mean[cpt])
        else:
            axes[i].plot(df[cpt])
        axes[i].set_ylabel(f"{cpt} (nT)")
        axes[i].grid()
    fig.suptitle(title)
    axes[2].set_xlabel("Date")
    fig.tight_layout()
    return fig

def make_widgets_minute_data():
    """Use ipywidgets to interact with plot_subset_timeseries"""
    mini, maxi = obs_minute.index.min().date(), obs_minute.index.max().date()
    start_date = widgets.DatePicker(
        value=mini,
        description='Start Date',
    )
    end_date = widgets.DatePicker(
        value=dt.datetime(mini.year, mini.month+1, mini.day),
        description='End Date',
    )
    hourly_mean = widgets.Checkbox(
        value=True,
        description='Hourly mean',
    )
    annual_mean = widgets.Checkbox(
        value=True,
        description='Show annual mean',
    )
    return widgets.VBox(
        [widgets.Label(f"Select dates within range: {mini}, {maxi}"),
         widgets.HBox([start_date, end_date, widgets.VBox([hourly_mean, annual_mean])]),
         widgets.interactive_output(
             plot_subset_timeseries,
             {'start_date': start_date, 'end_date': end_date, 'hourly_mean': hourly_mean, 'show_annual_mean': annual_mean}
         )]
    )

make_widgets_minute_data()

Note

What signals can you see in this data?

  • Daily oscillation: due to the rotation of the Earth driving ionospheric change through the day/night - this is the Sq variation (“solar quiet-day” variation)

  • Shift in baseline over the year: due to the change in the main magnetic field from the core - this is the secular variation (SV)

  • More random variations due to geomagnetic activity

Daily, seasonal and solar variations in declination

Let’s now fetch the hourly dataset - these data are specially processed to improve data quality, over the straightforward hourly means calculated above from the minute data. For more information, see Macmillan, S., Olsen, N. Observatory data and the Swarm mission. Earth Planet Sp 65, 15 (2013). https://doi.org/10.5047/eps.2013.07.011

obs_hourly = fetch_obs_data_for_years(
    observatory="ESK",
    year_start=1900,
    year_end=2020,
    cadence="H",
    # additional kwargs for vireslient.SwarmRequest.get_between()
    asynchronous=False,  # Make synchronous requests (faster)
                         #  - only works for smaller data chunks
                         #  - implicitly disables "Processing" progress bar
    show_progress=False,        # Disable intermediate progress bars
#     leave_intermediate_progress_bars=False,  # Clean up bars as we go
#     show_progress_chunks=False  # Disable "Processing chunks" progress bar
)

We will evaluate the declination angle, D, the horizontal deviation of the field from geographic North (what are the geomagnetic components?)

Next we summarise the data further by aggregating measurements over each month, evaluating the mean values over hourly intervals. For example, the mean declination at 10am across all days in January, the mean at 11am, and so on, repeated for each time of day and for each month. We then evaluate the offset of these declinations from the mean over the whole of each month - this is stored in D_variation in the resulting dataframe.

def monthly_means(df=obs_hourly):
    """Return MultIndex DataFrame of monthly means over each hourly interval"""
    # Append hour of day, and approx fractional year, to use for plotting
    df["t_hour"] = df.index.map(lambda x: x.hour + x.minute/60)
    epoch = pd.to_datetime(0, unit='s').to_julian_date()
    df["t_year"] = df.index.map(lambda x: x.year + (x.month-.5)/12)
    # Calculate the monthly mean for each hourly interval of the day
    monthly = df.groupby([df.index.year, df.index.month, df.index.hour]).mean()
    monthly.index.names = ["Year", "Month", "Hour"]
    # Calculate the monthly mean over all hourly intervals
    monthly_all = df.groupby([df.index.year, df.index.month]).mean()
    # Calculate the daily declination variations:
    #   the monthly average of the hourly intervals minus the total monthly mean
    monthly["D_variation"] = monthly['D'].values - monthly_all['D'].values.repeat(24)
    return monthly

obs_hourly["D"] = np.rad2deg(np.arctan2(obs_hourly["Y"], obs_hourly["X"]))
obs_monthly = monthly_means(obs_hourly)
obs_monthly
Radius Latitude Longitude X Y Z D t_hour t_year D_variation
Year Month Hour
1911 1 0 6363961.0 55.136 356.8 15864.441935 -5270.635484 45410.229032 -18.378015 0.5 1911.041667 0.035155
1 6363961.0 55.136 356.8 15861.590323 -5273.319355 45405.193548 -18.389796 1.5 1911.041667 0.023375
2 6363961.0 55.136 356.8 15859.645161 -5276.290323 45404.948387 -18.401581 2.5 1911.041667 0.011589
3 6363961.0 55.136 356.8 15860.135484 -5282.454839 45403.725806 -18.421103 3.5 1911.041667 -0.007933
4 6363961.0 55.136 356.8 15861.967742 -5279.909677 45406.032258 -18.410849 4.5 1911.041667 0.002322
... ... ... ... ... ... ... ... ... ... ... ... ...
2020 12 19 6363961.0 55.136 356.8 17451.764516 -465.945161 46683.003226 -1.529375 19.5 2020.958333 0.000556
20 6363961.0 55.136 356.8 17451.058065 -461.709677 46683.667742 -1.515541 20.5 2020.958333 0.014390
21 6363961.0 55.136 356.8 17450.429032 -459.390323 46683.912903 -1.507985 21.5 2020.958333 0.021946
22 6363961.0 55.136 356.8 17451.012903 -458.819355 46683.193548 -1.506063 22.5 2020.958333 0.023869
23 6363961.0 55.136 356.8 17452.741935 -459.187097 46681.954839 -1.507126 23.5 2020.958333 0.022806

31272 rows × 10 columns

def plot_dec_variation(df=obs_monthly, year_range=(1986, 1997)):
    """Make surface plot of declination variation against UT and Year"""
    year_start, year_end = year_range
    # Make subselection to plot, slicing along years
    df = df.loc[slice(year_start, year_end), :]
    fig, ax = plt.subplots(1, 1, figsize=(10, 10), subplot_kw={"projection": "3d"})
    ax.plot_trisurf(
        df["t_year"], df["t_hour"], df["D_variation"],
        cmap=plt.cm.jet, vmin=-0.15, vmax=0.15, antialiased=True
    )
    ax.set_xlabel("Year")
    ax.set_ylabel("Hour (UT)")
    ax.set_zlabel("D variations (degrees)")
    ax.view_init(elev=60, azim=210)
    return fig, ax

def make_widgets_dec_variation():
    year_min = obs_monthly.index[0][0]
    year_max = obs_monthly.index[-1][0]
    year_range = widgets.IntRangeSlider(
        value=(1986, 2008),
        min=year_min, max=year_max,
        description="Year range:"
    )
    return widgets.VBox([
        year_range,
        widgets.interactive_output(
            plot_dec_variation,
            {"year_range": year_range}
        )
    ])

make_widgets_dec_variation()

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

def bartels_rotation(datetime):
    """Return the Bartels rotation numbers and day and hour within each rotation

    Args:
        date (datetime/DatetimeIndex)

    Returns:
        tuple (rotation number, day within rotation, hour within day)

    """
    if isinstance(datetime, pd.DatetimeIndex):
        date, hour = datetime.date, datetime.hour
    elif isinstance(datetime, 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
    return bartels_rotation, day, hour

def reform_dataframe_for_Bartels_plot(df=obs_hourly):
    """Reforms an observatory dataframe ready for the Bartels plot

    1. Appends annual means
    2. Replaces the DatetimeIndex with a MultiIndex of rotation number and day

    Args:
        df (DataFrame)

    Returns:
        DataFrame

    """
    df = df.copy()

    # Evaluate the annual means
    annual_means = df[["X", "Y", "Z"]].resample("1y").mean()
    annual_means = annual_means.rename(columns={"X": "X_mean", "Y": "Y_mean", "Z": "Z_mean"})
    # Identify the year-end index times
    year_ends = [df.loc[str(year)].index[-1] for year in df.index.year.unique()]
    # Remove unused years from annual means
    missing_years = set(annual_means.index.year) - set([i.year for i in year_ends])
    annual_means = annual_means.drop(
        [pd.to_datetime(f"{year}-12-31") for year in missing_years]
    )
    # Match annual_means index to input index
    annual_means.index = year_ends
    # Infill the means into the input dataframe
    annual_means = annual_means.reindex(index=df.index, method="backfill")
    df = df.join(annual_means)
    
    # Preserve the original index as a column
    df["time"] = df.index
    # Replace the index with a MultiIndex of the Bartels rotation number
    # and the fractional day within each rotation
    bartrot, bartrotday, hourofday = bartels_rotation(df.index)
    df.index = pd.MultiIndex.from_arrays(
        (bartrot, bartrotday + hourofday/24), names=("bartrot", "bartrotday")
    )
    return df

def plot_by_Bartels_rotation(df=obs_hourly, year=2003, var="Y"):
    """Plot the hourly data organised by Bartels rotation"""
    # Trim dataframe to chosen year and reform it
    df = reform_dataframe_for_Bartels_plot(obs_hourly.loc[str(year)])

    bartrots = range(df.index[0][0], df.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 = df.loc[bartrot].index
        y0 = df.loc[bartrot][f"{var}_mean"]
        y1 = df.loc[bartrot][f"{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
        month_starts = df["time"].where(
            (df["time"].dt.day == 1) & (df["time"].dt.hour == 0)
        ).dropna()
        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"{df['IAGA_code'].iloc[0]}: {year}: {var} (nT, left axis), by Bartels rotation number (right axis)",
        fontsize=15
    )
    return fig, axes

def make_widgets_Bartels_plot():
    year_min = obs_hourly.index.year.min()
    year_max = obs_hourly.index.year.max()
    year = widgets.IntSlider(value=2003, min=year_min, max=year_max, description="Year:")
    component = widgets.Dropdown(options="XYZ", value="Y", description="Component:")
    return widgets.VBox([
        widgets.HBox([year, component]),
        widgets.interactive_output(
            plot_by_Bartels_rotation,
            {"year": year, "var": component}
        )
    ])

make_widgets_Bartels_plot()

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.

We will fetch hourly data again, but this time for three observatories: Eskdalemuir [ESK], Niemegk [NGK] and Chambon-la-Foret [CLF]. We evaluate the annual mean and the year-to-year differences (i.e. the annual secular variation, in nT/year). Annual means are also available directly from http://www.geomag.bgs.ac.uk/data_service/data/annual_means.shtml

def fetch_annual_means(observatories):
    """Collect annual means from given observatories"""
    annual_means = {}
    for i, obs in enumerate(observatories):
        print(f"Fetching {obs} [{i+1}/{len(observatories)}]")
        df = fetch_obs_data_for_years(
            observatory=obs,
            year_start=1900,
            year_end=2020,
            cadence="H",
            asynchronous=False,
            show_progress=False,
        )
        df = df.resample("1y").mean()
        df.index = df.index.year
        df.index.name = "Year"
        annual_means[obs] = df
    return annual_means

observatories = ["ESK", "NGK", "CLF"]
annual_means = fetch_annual_means(observatories)
Fetching ESK [1/3]
Fetching NGK [2/3]
Fetching CLF [3/3]
# Append annual differences
for obs in annual_means.keys():
    # First differences = first derivative = secular variation
    diffs = annual_means[obs][["X", "Y", "Z"]].diff()
    # Second differences = seconda derivative = secular acceleration
    diffs2 = diffs.diff()
    diffs = diffs.rename(columns={"X": "dX", "Y": "dY", "Z": "dZ"})
    diffs2 = diffs2.rename(columns={"X": "ddX", "Y": "ddY", "Z": "ddZ"})
    annual_means[obs] = annual_means[obs].join(diffs)
    annual_means[obs] = annual_means[obs].join(diffs2)

# Now each dataframe from each observatory looks like this:
annual_means["ESK"]
Radius Latitude Longitude X Y Z dX dY dZ ddX ddY ddZ
Year
1911 6363961.0 55.136 356.8 15863.714996 -5264.613667 45396.045733 NaN NaN NaN NaN NaN NaN
1912 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
1913 6363961.0 55.136 356.8 15867.208015 -5166.959745 45326.189360 NaN NaN NaN NaN NaN NaN
1914 6363961.0 55.136 356.8 15864.780980 -5125.410973 45239.835457 -2.427035 41.548772 -86.353903 NaN NaN NaN
1915 6363961.0 55.136 356.8 15862.133080 -5076.216160 45224.384173 -2.647900 49.194814 -15.451284 -0.220866 7.646042 70.902619
... ... ... ... ... ... ... ... ... ... ... ... ...
2016 6363961.0 55.136 356.8 17391.096334 -751.640710 46534.882001 17.924063 58.215682 21.308725 12.765204 0.585545 -9.712919
2017 6363961.0 55.136 356.8 17405.208105 -687.978196 46563.534612 14.111771 63.662514 28.652611 -3.812292 5.446832 7.343885
2018 6363961.0 55.136 356.8 17421.728721 -624.448973 46593.312671 16.520616 63.529224 29.778059 2.408846 -0.133290 1.125449
2019 6363961.0 55.136 356.8 17436.509521 -558.797374 46627.251849 14.780799 65.651598 33.939178 -1.739817 2.122374 4.161119
2020 6363961.0 55.136 356.8 17447.969843 -494.473782 46661.702789 11.460322 64.323593 34.450940 -3.320477 -1.328006 0.511762

110 rows × 12 columns

annual_means.keys()
dict_keys(['ESK', 'NGK', 'CLF'])
def plot_secular_variation(df_set=annual_means, observatory="ESK"):
    # Select dataframe for chosen observatory
    df = df_set[observatory]
    fig, axes = plt.subplots(nrows=6, figsize=(10,7), sharex=True)
    colors = plt.rcParams["axes.prop_cycle"].by_key()["color"]
    labels = {
        "dX": "dX/dt (nT/yr)", "dY": "dY/dt (nT/yr)", "dZ": "dZ/dt (nT/yr)",
        "ddX": "d$^2$X/dt$^2$", "ddY": "d$^2$Y/dt$^2$", "ddZ": "d$^2$Z/dt$^2$"
    }
    for (i, (var1, var2)) in enumerate(zip(["dX", "dY", "dZ"], ["ddX", "ddY", "ddZ"])):
        axes[2*i].plot(df[var1], color=colors[i], marker="x")
        axes[2*i].set_ylabel(labels[var1])
        axes[2*i+1].plot(df[var2], color=colors[i], alpha=0.5, marker="x")
        axes[2*i+1].set_ylabel(labels[var2])
        axes[2*i].grid()
        axes[2*i+1].grid()
    axes[-1].set_xlabel("Year")
    fig.suptitle(f"Secular variation at {observatory}")

def make_widgets_secular_variation():
    obs_choice = widgets.Dropdown(options=annual_means.keys())
    return widgets.VBox([
        obs_choice,
        widgets.interactive_output(
            plot_secular_variation,
            {"observatory": obs_choice}
        )
    ])

make_widgets_secular_variation()

A sudden change in the field is quite clear in the secular variation (e.g. dY/dt for ESK around 1970), but the change is not obvious in the secular acceleration in the way that we have calculated it (one would expect a sudden jump from negative to positive acceleration there)! Some more careful processing is needed in order to make the change in the secular acceleration clear.