Making a main field model

Making a main field model#

Here we demonstrate how to make a simple time-independent main field model by fitting spherical harmonics to ground observatory measurements.

The principles of spherical harmonic analysis have been shown in the previous pages, so here we will take a shortcut and use utilties from ChaosMagPy.

ChaosMagPy, Clemens Kloss, https://doi.org/10.5281/zenodo.3352398
We will use:
- design_gauss to generate the design matrix for the inversion
- synth_values to run the forward model to get predictions for the magnetic field from our model

hvPlot is used to create some fancy visualisations but this library is more complex to use! so if you are newer to the Python ecosystem, you are better off using Matplotlib for your own work.

# Install pooch (not currently in VRE)
import sys
!{sys.executable} -m pip install pooch

import pooch
from chaosmagpy.model_utils import design_gauss, synth_values, power_spectrum
from chaosmagpy.plot_utils import plot_power_spectrum
import numpy as np
import pandas as pd
import xarray as xr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import hvplot.xarray

Building a model#

Suppose we have data, \(\vec{d}\) (i.e. the set of magnetic measurements), and we want to find a model \(\vec{m}\) (i.e. a set of Gauss coefficients, \(\begin{bmatrix} g_n^m & h_n^m\end{bmatrix}\), to be determined), we can connect these with a matrix \(G\) (the so-called design matrix):

\[ \vec{d} = G \vec{m} \]

\(G\) can be constructed by comparing to the spherical harmonic expansion of the magnetic field \(\vec{B}\) (where we have neglected external sources for simplicity):

\[\begin{split} \begin{align}\label{eq:spharm_B} \vec{d} =\ & G \vec{m} \\ B_i =\ & G \begin{bmatrix} g_n^m & h_n^m \end{bmatrix} \\ B_i =\ & -\partial_i \left( R_E \sum\limits_{n=1}^{N} \sum\limits_{m=0}^{n} \left( g_n^m \cos{m\phi} + h_n^m \sin{m\phi} \right) \left(\frac{R_E}{r}\right)^{n+1} P_n^m (\cos{\theta}) \right) \nonumber \\ =\ & -\partial_i \begin{bmatrix} R_E \sum\limits_{n=1}^{N} \sum\limits_{m=0}^{n} \left( \cos{m\phi} \right) \left(\frac{R_E}{r}\right)^{n+1} P_n^m (\cos{\theta}) \\ R_E \sum\limits_{n=1}^{N} \sum\limits_{m=0}^{n} \left( \sin{m\phi} \right) \left(\frac{R_E}{r}\right)^{n+1} P_n^m (\cos{\theta}) \end{bmatrix} \begin{bmatrix} g_n^m & h_n^m \end{bmatrix} \end{align} \end{split}\]

Note that:

\(B_i\) are the three vector components, \((B_r, B_\theta, B_\phi) = (-B_C, -B_N, B_E)\), i.e. the data, \(\vec{d}\)
\(\partial_i\) is a shorthand which should be replaced by the derivatives in spherical polar coordinates \( \left( \partial_r = \frac{\partial}{\partial r}, \partial_\theta = \frac{1}{r} \frac{\partial}{\partial \theta}, \partial_\phi = \frac{1}{r\sin{\theta}} \frac{\partial}{\partial \phi} \right) \)
\(\begin{bmatrix} g_n^m & h_n^m \end{bmatrix}\) are the Gauss coefficients, i.e. the model, \(\vec{m}\)

A least-squares solution to \(\vec{d} = G \vec{m}\) can be found with:

\[\vec{m} = (G^T G)^{-1} G^T \vec{d}\]

def build_model(ds=input_data, nmax=10):
    """Returns Gauss coefficients for a simple SH model"""
    theta = ds["Colat"]
    phi = ds["Long"]
    radius = 6371.2 + ds["alt(m)"]/1e3
    B_radius = -ds["Z"]
    B_theta = -ds["X"]
    B_phi = ds["Y"]
    # Build the design matrix
    G_radius, G_theta, G_phi = design_gauss(radius, theta, phi, nmax=nmax, source="internal")
    # Perform the inversion
    G = np.concatenate((G_radius, G_theta, G_phi))
    d = np.concatenate((B_radius, B_theta, B_phi))
    m = np.linalg.inv(G.T @ G) @ (G.T @ d)
    return m


m = build_model(input_data, nmax=5)
m

array([-2.96482406e+04, -1.76621579e+03,  5.06863059e+03, -2.47450873e+03,
        3.13168299e+03, -2.64941275e+03,  1.74507015e+03, -3.76171556e+02,
        1.10537201e+03, -2.25018745e+03, -6.30954563e+01,  1.19694272e+03,
        3.54883214e+02,  8.12814264e+02, -5.15153390e+02,  9.99590132e+02,
        9.11464975e+02,  1.54019008e+02,  3.42255212e+02, -3.26457915e+02,
       -3.49774227e+02,  2.34546359e+02,  7.81036147e+01, -3.73915172e+02,
       -2.86846063e+02,  1.90610037e+02, -1.94545936e+02,  3.69767832e+02,
        1.49286822e+02, -1.94868808e+02, -7.07565425e+01, -1.49668101e+02,
       -2.26259266e+01, -5.93299060e+01,  1.62121406e+02])

Since we do not have many data, we would need to be more sophisticated with the inversion procedure (e.g. including regularisation) to get a model of higher resolution. For this reason, the model is truncated at spherical harmonic degree 6. You can explore increasing this limit.

How good is our model?#

Evaluating a model over the Earth and plotting contours…

def sample_model_on_grid(m, radius=6371.200):
    """Evaluate Gauss coefficients over a regular grid over Earth"""
    theta = np.arange(1, 180, 1)
    phi = np.arange(-180, 180, 1)
    theta, phi = np.meshgrid(theta, phi)
    B_r, B_t, B_p = synth_values(m, radius, theta, phi)
    ds_grid = xr.Dataset(
        data_vars={"B_NEC_model": (("y", "x", "NEC"), np.stack((-B_t, B_p, -B_r), axis=2))},
        coords={"Latitude": (("y", "x"), 90-theta), "Longitude": (("y", "x"), phi), "NEC": np.array(["N", "E", "C"])}
    )
    return ds_grid


def plot_model_contours(m):
    ds_grid = sample_model_on_grid(m)
    # Generate contour plot of vertical component
    return ds_grid.sel(NEC="C").hvplot.contourf(
        x="Longitude", y="Latitude", c="B_NEC_model",
        levels=30, coastline=True, projection=ccrs.PlateCarree(), global_extent=True,
        clabel="Model: vertical (downwards) magnetic field (nT)"
    )


plot_model_contours(m)

/opt/conda/lib/python3.9/site-packages/geoviews/operation/projection.py:79: ShapelyDeprecationWarning: Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry.
  polys = [g for g in geom if g.area > 1e-15]

The power spectrum is a common way to assess and compare models, showing how much power is concentrated at each spherical harmonic degree:

plot_power_spectrum(power_spectrum(m));

../_images/72feb276b047fc923057fba71628e54f98bfa9eba32c540b4e7a54939321b2e6.png

… and we can look at the residuals between the input data and the model:

def append_residuals(ds=input_data, m=m):
    # Evaluate model at the data locations
    B_r, B_t, B_p = synth_values(m, 6371.2 + ds["alt(m)"]/1e3, ds["Colat"], ds["Longitude"])
    # Assign data-model residuals to the dataset
    ds["res_N"] = ds["X"] - (-B_t)
    ds["res_E"] = ds["Y"] - (B_p)
    ds["res_C"] = ds["Z"] - (-B_r)
    return ds

ds = append_residuals(input_data, m)

plt.hist(ds["res_C"])
plt.xlabel("residual");

../_images/77db50a57bf16c4530f761692aec39f42fbcb720a53e337702d46d777e82cd79.png

ds.hvplot.scatter(
    x="Longitude", y="Latitude", c="res_C",
    coastline=True, projection=ccrs.PlateCarree(), global_extent=True,
    clabel="Data-model residual: vertical (downwards) magnetic field (nT)"
)

WARNING:param.Parameterized: Use method 'warning' via param namespace

WARNING:param.main: geo option cannot be used with kind='scatter' plot type. Geographic plots are only supported for following plot types: ['bivariate', 'contour', 'contourf', 'dataset', 'hexbin', 'image', 'labels', 'paths', 'points', 'points', 'polygons', 'quadmesh', 'rgb', 'vectorfield']

Exercises#

Experiment with changing the SH degree trunction in the model
Experiment with using data just over Europe
Add random noise into the data to observe the effect
Use viresclient to access the monthly means dataset and produce a model using those instead

def select_europe(ds=input_data):
    """Subselect the dataset down to just cover Europe"""
    colat_minmax = (20, 60)
    lon_minmax = (-10, 40)
    ds = ds.where(
        (ds["Colat"] > colat_minmax[0]) & (ds["Colat"] < colat_minmax[1]) &
        (ds["Long"] > lon_minmax[0]) & (ds["Long"] < lon_minmax[1]),
    ).dropna(dim="index")
    return ds

Making a main field model

Contents

Making a main field model#

Loading some data#

Building a model#

How good is our model?#

Exercises#