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.

Hide code cell content 
## Uncomment here to autoreload local modules as you edit them
# %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

Python implementation: CPython
Python version       : 3.9.7
IPython version      : 8.0.1

numpy      : 1.21.5
pandas     : 1.4.1
xarray     : 2023.3.0
matplotlib : 3.5.1
ipywidgets : 7.6.5
viresclient: 0.11.1

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. Check http://intermagnet.org/new_data_download.html for other ways to access these data.

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.

Hide code cell source 
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 IAGA_code Longitude X Y Z
Timestamp
2003-01-01 00:00:00 6.363967e+06 55.119672 ESK 356.8 17196.514832 -1473.2 46252.152020
2003-01-01 00:01:00 6.363967e+06 55.119672 ESK 356.8 17196.014835 -1473.4 46252.150446
2003-01-01 00:02:00 6.363967e+06 55.119672 ESK 356.8 17196.514832 -1473.4 46252.152020
2003-01-01 00:03:00 6.363967e+06 55.119672 ESK 356.8 17196.715461 -1473.8 46251.952650
2003-01-01 00:04:00 6.363967e+06 55.119672 ESK 356.8 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.

Hide code cell source 
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

Hide code cell content 
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.

Hide code cell source 
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.619355 -466.096774 46683.338710 -1.529885 19.5 2020.958333 0.000651
20 6363961.0 55.136 356.8 17450.877419 -461.903226 46683.909677 -1.516192 20.5 2020.958333 0.014344
21 6363961.0 55.136 356.8 17450.264516 -459.548387 46684.096774 -1.508518 21.5 2020.958333 0.022018
22 6363961.0 55.136 356.8 17450.845161 -459.032258 46683.454839 -1.506776 22.5 2020.958333 0.023760
23 6363961.0 55.136 356.8 17452.651613 -459.387097 46682.177419 -1.507789 23.5 2020.958333 0.022747

31272 rows × 10 columns

Hide code cell source 
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).

Hide code cell source 
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

Hide code cell content 
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]

/tmp/ipykernel_89/3898201863.py:14: FutureWarning: Dropping invalid columns in DataFrameGroupBy.mean is deprecated. In a future version, a TypeError will be raised. Before calling .mean, select only columns which should be valid for the function.
  df = df.resample("1y").mean()

Fetching CLF [3/3]

/tmp/ipykernel_89/3898201863.py:14: FutureWarning: Dropping invalid columns in DataFrameGroupBy.mean is deprecated. In a future version, a TypeError will be raised. Before calling .mean, select only columns which should be valid for the function.
  df = df.resample("1y").mean()
Hide code cell content 
# 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.909176 -494.591985 46661.667919 11.399655 64.205389 34.416070 -3.381144 -1.446209 0.476892

110 rows × 12 columns

annual_means.keys()

dict_keys(['ESK', 'NGK', 'CLF'])
Hide code cell source 
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.