Commit 4758c6e7 authored by Maximilian Schanner's avatar Maximilian Schanner
Browse files

Add an example notebook for custom models.

parent dbf66d82
......@@ -123,10 +123,11 @@ test-debian:
pages:
stage: deploy
script:
- cp examples/*.ipynb docs/
- apt-get update -y -qq
- apt-get install -y -qq python3-dev python3-pip python3-cartopy
- apt-get install -y -qq python3-dev python3-pip python3-cartopy pandoc
- pip install pymagglobal --extra-index-url https://public:5mz_iyigu-WE3HySBH1J@git.gfz-potsdam.de/api/v4/projects/${CI_PROJECT_ID}/packages/pypi/simple
- pip install sphinx sphinx-argparse
- pip install jupyter notebook sphinx nbsphinx sphinx-argparse
- sphinx-build -b html docs public
artifacts:
paths:
......
.. _examples-label:
Example Notebooks
=================
Here we collect relevant examples. The list may grow over time.
.. toctree::
:maxdepth: 2
example_1
......@@ -9,6 +9,7 @@ It can be applied to all cubic-spline based geomagnetic field models stored in t
package_documentation
command_line_interface
examples
Installation
......
%% Cell type:markdown id:681fb4de tags:
# Using pymagglobal with custom models
%% Cell type:markdown id:7c938d6b tags:
This notebook focusses on how to translate a custom model to `pymagglobal`. We focus on the SHA.DIF.14k-model by [Pavón-Carrasco et al.](#pavon), which is publicly available [here](http://pc213fis.fis.ucm.es/sha.dif.14k/).
> <a name="pavon"></a> Pavón-Carrasco, F.J., M. L. Osete, J. M. Torta and A. De Santis (2014)
> A Geomagnetic Field Model for the Holocene based on archaeomagnetic
> and lava flow data. Earth Planet. Sci. Lett., 388, 98 - 109.
To access the model, we download it directly using `urllib` and `numpy`:
%% Cell type:code id:dbcb982b tags:
``` python
import numpy as np
from urllib.request import urlopen
raw_data = urlopen('http://pc213fis.fis.ucm.es/sha.dif.14k/coeff_SHA.DIF.14k.dat')
data = np.genfromtxt(raw_data).T
```
%% Cell type:markdown id:837cd058 tags:
With the additional knowledge that the model is valid in the interval `[-12000, 1800]`, we can set up an object that can be used by `pymagglobal`:
%% Cell type:code id:f4b5705d tags:
``` python
from pymagglobal import Model
from pymagglobal.utils import i2lm_l, i2lm_m, lmax2N
from scipy.interpolate import BSpline
# Inherit from pymagglobals Model-class
class SHADIF14k(Model):
# Override the __init__ method to set up
# the model with our data
def __init__(self):
# set up the name
self.name = 'sha.dif.14k'
# The first column consists of the knots,
# but it contains many duplications
self.knots = np.unique(data[0])
# Set the interval in which the model is valid
self.t_min = -12000
self.t_max = 1800
# The first column contains the degrees, so
# we extract the maxum
self.l_max = int(np.max(data[1]))
# Now we go through the data and access the
# coefficients line by line
# First set up an empty array of the right size
self.coeffs = np.empty((len(self.knots), lmax2N(self.l_max)))
# Then iterate over all coefficients
for it in range(lmax2N(self.l_max)):
# Get degree and order from the index
ell = i2lm_l(it)
m = i2lm_m(it)
# Get the corresponding indices in the data-array
inds = (data[1] == ell) * (data[2] == abs(m))
# Write the data to the coefficient-array
if m < 0:
# scale to nT
self.coeffs[:, it] = data[4][inds]*1000
else:
# scale to nT
self.coeffs[:, it] = data[3][inds]*1000
# Finally initialize the BSpline object
self.splines = BSpline(self.knots,
self.coeffs,
3)
# The model file does not report uncertainties
self.cov_splines = None
```
%% Cell type:markdown id:4a463586 tags:
This model can now be initialized and used like any other model in `pymagglobal`:
%% Cell type:code id:d891c8d5 tags:
``` python
from matplotlib import pyplot as plt
from pymagglobal import dipole_series
plt.rcParams['font.size'] = '14'
shadif14k = SHADIF14k()
times = np.linspace(-1000, 1500, 501)
fig, ax = plt.subplots(1, 1, figsize=(10, 6))
for color, model in zip(['C0', 'C1'], [shadif14k, Model('ARCH10k.1')]):
dip = dipole_series(times, model)
ax.plot(times, dip, color=color, zorder=2, label=model.name)
ax.legend();
ax.set_ylabel(r'Dipole moment [$10^{22}$Am$^2$]');
ax.set_xlabel(r'time [yrs.]');
```
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment