{
"cells": [
{
"cell_type": "markdown",
"id": "15286a2e-6ccf-456b-b41b-79b68657b2e1",
"metadata": {
"tags": []
},
"source": [
"# Making a main field model"
]
},
{
"cell_type": "markdown",
"id": "dcbb6893-825f-4b3b-af2e-47b84dbab5ee",
"metadata": {},
"source": [
"Here we demonstrate how to make a simple time-independent main field model by fitting spherical harmonics to ground observatory measurements.\n",
"\n",
"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.\n",
"\n",
"- [ChaosMagPy](https://chaosmagpy.readthedocs.io/), Clemens Kloss, https://doi.org/10.5281/zenodo.3352398 \n",
" We will use:\n",
" - [design_gauss](https://chaosmagpy.readthedocs.io/en/master/functions/chaosmagpy.model_utils.design_gauss.html) to generate the design matrix for the inversion\n",
" - [synth_values](https://chaosmagpy.readthedocs.io/en/master/functions/chaosmagpy.model_utils.synth_values.html) to run the forward model to get predictions for the magnetic field from our model\n",
"\n",
"[hvPlot](https://hvplot.holoviz.org/) 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."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "030e9e9a-4a31-4fe4-be18-3824b527f905",
"metadata": {
"tags": [
"hide-output"
]
},
"outputs": [],
"source": [
"# Install pooch (not currently in VRE)\n",
"import sys\n",
"!{sys.executable} -m pip install pooch"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "c0fdb8b6-cf3d-4f3d-a64b-d9a10b5dd3d8",
"metadata": {},
"outputs": [],
"source": [
"import pooch\n",
"from chaosmagpy.model_utils import design_gauss, synth_values, power_spectrum\n",
"from chaosmagpy.plot_utils import plot_power_spectrum\n",
"import numpy as np\n",
"import pandas as pd\n",
"import xarray as xr\n",
"import matplotlib.pyplot as plt\n",
"import cartopy.crs as ccrs\n",
"import cartopy.feature as cfeature\n",
"import hvplot.xarray"
]
},
{
"cell_type": "markdown",
"id": "83a29805-1b4e-42b1-9357-aa35548bcc34",
"metadata": {},
"source": [
"## Loading some data"
]
},
{
"cell_type": "markdown",
"id": "638c7a53-8f91-4563-9dd0-66be049c7358",
"metadata": {},
"source": [
"We'll load a pre-prepared dataset containing annual means derived from observatory data. We just select a single year from this dataset. In reality the main field varies on shorter time scales but these annual data points are useful enough for a low resolution model."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "33e38521-db52-4d15-a809-18b848ac7159",
"metadata": {},
"outputs": [],
"source": [
"# Download annual means dataset\n",
"# https://geomag.bgs.ac.uk/data_service/data/annual_means.shtml\n",
"oamjumpsapplied = pooch.retrieve(\n",
" \"https://raw.githubusercontent.com/MagneticEarth/IAGA_SummerSchool2019/master/data/external/oamjumpsapplied.dat\",\n",
" known_hash=\"dee77ded79cb0d86e3366a550830cad9de741322212faf9a74f3d03cb6bdd872\"\n",
")"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "3ccd7f46-c265-4c7c-a653-e915fd27a310",
"metadata": {},
"outputs": [],
"source": [
"def load_data(year=2000.5):\n",
" df = pd.read_csv(oamjumpsapplied, delim_whitespace=True)\n",
" df = df.where(df[\"year\"] == 2000.5).dropna()\n",
" df = df.reset_index(drop=True)\n",
" # Screen out some bad data points\n",
" df = df.where(df[\"Z\"] != 88888).dropna()\n",
" _ds = df.to_xarray()\n",
" _ds[\"Latitude\"] = 90 - _ds[\"Colat\"]\n",
" _ds[\"Longitude\"] = (_ds[\"Long\"] % 360 + 540) % 360 - 180\n",
" _ds = _ds.drop(\"index\")\n",
" return _ds\n",
"\n",
"input_data = load_data(year=2000.5)\n",
"input_data"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "74e9e0ab-06b4-4a6b-aa68-580068d2b094",
"metadata": {},
"outputs": [],
"source": [
"def plot_obs_measurements(ds=input_data):\n",
" return ds.hvplot.scatter(\n",
" x=\"Longitude\", y=\"Latitude\", c=\"Z\",\n",
" coastline=True, projection=ccrs.PlateCarree(), global_extent=True,\n",
" clabel=\"Data: vertical (downwards) magnetic field (nT)\"\n",
" )\n",
"\n",
"\n",
"plot_obs_measurements(input_data)"
]
},
{
"cell_type": "markdown",
"id": "59d32c3d-63b6-4241-b4e5-f74f42b69b90",
"metadata": {
"tags": []
},
"source": [
"## Building a model"
]
},
{
"cell_type": "markdown",
"id": "14b0127c-c4ef-4bdc-8133-547461341d7a",
"metadata": {},
"source": [
"<!-- We assume that measurements are taken within a shell free of electric currents / magnetic sources, and so the magnetic field can be simply expressed as the gradient of its scalar potential, $\\vec{B} = - \\nabla V$. Following from Maxwell's equations, we also have $\\nabla . \\vec{B} = 0$, and so we find that the potential must be a solution to $\\nabla^2 V = 0$ (Laplace's equation). It can be shown that solutions can be expressed as sums of spherical harmonics. -->\n",
"\n",
"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*):\n",
"\n",
"$$\n",
" \\vec{d} = G \\vec{m}\n",
"$$\n",
"\n",
"$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):\n",
"\n",
"$$\n",
"\\begin{align}\\label{eq:spharm_B}\n",
" \\vec{d} =\\ & G \\vec{m} \\\\\n",
" B_i =\\ & G \\begin{bmatrix}\n",
" g_n^m & h_n^m\n",
" \\end{bmatrix} \\\\\n",
" 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 \\\\\n",
" =\\ & -\\partial_i \\begin{bmatrix} \n",
"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}) \\\\\n",
"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})\n",
"\\end{bmatrix}\n",
"\\begin{bmatrix}\n",
"g_n^m & h_n^m\n",
"\\end{bmatrix}\n",
"\\end{align}\n",
"$$\n",
"\n",
"Note that:\n",
"- $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}$\n",
"- $\\partial_i$ is a shorthand which should be replaced by the derivatives in spherical polar coordinates \n",
"$\n",
"\\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)\n",
"$\n",
"- $\\begin{bmatrix} g_n^m & h_n^m \\end{bmatrix}$ are the Gauss coefficients, i.e. the model, $\\vec{m}$\n",
"\n",
"\n",
"A least-squares solution to $\\vec{d} = G \\vec{m}$ can be found with:\n",
"\n",
"$$\\vec{m} = (G^T G)^{-1} G^T \\vec{d}$$"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "f776207a-0392-4d2a-bee7-b6dcdff4b82e",
"metadata": {},
"outputs": [],
"source": [
"def build_model(ds=input_data, nmax=10):\n",
" \"\"\"Returns Gauss coefficients for a simple SH model\"\"\"\n",
" theta = ds[\"Colat\"]\n",
" phi = ds[\"Long\"]\n",
" radius = 6371.2 + ds[\"alt(m)\"]/1e3\n",
" B_radius = -ds[\"Z\"]\n",
" B_theta = -ds[\"X\"]\n",
" B_phi = ds[\"Y\"]\n",
" # Build the design matrix\n",
" G_radius, G_theta, G_phi = design_gauss(radius, theta, phi, nmax=nmax, source=\"internal\")\n",
" # Perform the inversion\n",
" G = np.concatenate((G_radius, G_theta, G_phi))\n",
" d = np.concatenate((B_radius, B_theta, B_phi))\n",
" m = np.linalg.inv(G.T @ G) @ (G.T @ d)\n",
" return m\n",
"\n",
"\n",
"m = build_model(input_data, nmax=5)\n",
"m"
]
},
{
"cell_type": "markdown",
"id": "129456d9-e0a3-4ee9-a3b8-936c5976f8a3",
"metadata": {},
"source": [
"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."
]
},
{
"cell_type": "markdown",
"id": "53a9d939-726c-40f2-a596-dc04409060a2",
"metadata": {},
"source": [
"## How good is our model?"
]
},
{
"cell_type": "markdown",
"id": "312a3527-d2b4-4a11-b17d-84f9e9ec73bd",
"metadata": {},
"source": [
"Evaluating a model over the Earth and plotting contours..."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "90723ab9-cb19-484c-b2fe-73f891e49aa3",
"metadata": {},
"outputs": [],
"source": [
"def sample_model_on_grid(m, radius=6371.200):\n",
" \"\"\"Evaluate Gauss coefficients over a regular grid over Earth\"\"\"\n",
" theta = np.arange(1, 180, 1)\n",
" phi = np.arange(-180, 180, 1)\n",
" theta, phi = np.meshgrid(theta, phi)\n",
" B_r, B_t, B_p = synth_values(m, radius, theta, phi)\n",
" ds_grid = xr.Dataset(\n",
" data_vars={\"B_NEC_model\": ((\"y\", \"x\", \"NEC\"), np.stack((-B_t, B_p, -B_r), axis=2))},\n",
" coords={\"Latitude\": ((\"y\", \"x\"), 90-theta), \"Longitude\": ((\"y\", \"x\"), phi), \"NEC\": np.array([\"N\", \"E\", \"C\"])}\n",
" )\n",
" return ds_grid\n",
"\n",
"\n",
"def plot_model_contours(m):\n",
" ds_grid = sample_model_on_grid(m)\n",
" # Generate contour plot of vertical component\n",
" return ds_grid.sel(NEC=\"C\").hvplot.contourf(\n",
" x=\"Longitude\", y=\"Latitude\", c=\"B_NEC_model\",\n",
" levels=30, coastline=True, projection=ccrs.PlateCarree(), global_extent=True,\n",
" clabel=\"Model: vertical (downwards) magnetic field (nT)\"\n",
" )\n",
"\n",
"\n",
"plot_model_contours(m)"
]
},
{
"cell_type": "markdown",
"id": "153e2f1a-87e8-4dee-ac0e-baa5a6a52cf8",
"metadata": {},
"source": [
"The power spectrum is a common way to assess and compare models, showing how much power is concentrated at each spherical harmonic degree:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "80e05b5c-eec8-4e46-b3d6-c862774b302c",
"metadata": {},
"outputs": [],
"source": [
"plot_power_spectrum(power_spectrum(m));"
]
},
{
"cell_type": "markdown",
"id": "0845b195-28b1-4243-8888-2481145dd5d1",
"metadata": {},
"source": [
"... and we can look at the residuals between the input data and the model:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "509da5ec-0aa6-4a95-8a21-cc0ea9ec1c1e",
"metadata": {},
"outputs": [],
"source": [
"def append_residuals(ds=input_data, m=m):\n",
" # Evaluate model at the data locations\n",
" B_r, B_t, B_p = synth_values(m, 6371.2 + ds[\"alt(m)\"]/1e3, ds[\"Colat\"], ds[\"Longitude\"])\n",
" # Assign data-model residuals to the dataset\n",
" ds[\"res_N\"] = ds[\"X\"] - (-B_t)\n",
" ds[\"res_E\"] = ds[\"Y\"] - (B_p)\n",
" ds[\"res_C\"] = ds[\"Z\"] - (-B_r)\n",
" return ds\n",
"\n",
"ds = append_residuals(input_data, m)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "ea744dc5-7402-4d06-b578-cdaf05a43474",
"metadata": {},
"outputs": [],
"source": [
"plt.hist(ds[\"res_C\"])\n",
"plt.xlabel(\"residual\");"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "3328797d-f563-4d83-aab7-a097ed027615",
"metadata": {},
"outputs": [],
"source": [
"ds.hvplot.scatter(\n",
" x=\"Longitude\", y=\"Latitude\", c=\"res_C\",\n",
" coastline=True, projection=ccrs.PlateCarree(), global_extent=True,\n",
" clabel=\"Data-model residual: vertical (downwards) magnetic field (nT)\"\n",
")"
]
},
{
"cell_type": "markdown",
"id": "8b20d365-844e-49d9-ab4a-2823307504a6",
"metadata": {},
"source": [
"## Exercises\n",
"\n",
"- Experiment with changing the SH degree trunction in the model\n",
"- Experiment with using data just over Europe\n",
"- Add random noise into the data to observe the effect\n",
"- Use viresclient to access the monthly means dataset and produce a model using those instead"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "c7459fce-c5b5-44e9-895f-3d7771929d72",
"metadata": {},
"outputs": [],
"source": [
"def select_europe(ds=input_data):\n",
" \"\"\"Subselect the dataset down to just cover Europe\"\"\"\n",
" colat_minmax = (20, 60)\n",
" lon_minmax = (-10, 40)\n",
" ds = ds.where(\n",
" (ds[\"Colat\"] > colat_minmax[0]) & (ds[\"Colat\"] < colat_minmax[1]) &\n",
" (ds[\"Long\"] > lon_minmax[0]) & (ds[\"Long\"] < lon_minmax[1]),\n",
" ).dropna(dim=\"index\")\n",
" return ds"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"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.10.6"
}
},
"nbformat": 4,
"nbformat_minor": 5
}