example_1.ipynb 5.06 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "681fb4de",
   "metadata": {},
   "source": [
    "# Using pymagglobal with custom models"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7c938d6b",
   "metadata": {},
   "source": [
    "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/).\n",
    "\n",
    "> <a name=\"pavon\"></a> Pavón-Carrasco, F.J., M. L. Osete, J. M. Torta and A. De Santis (2014)  \n",
    "> A Geomagnetic Field Model for the Holocene based on archaeomagnetic  \n",
    "> and lava flow data. Earth Planet. Sci. Lett., 388, 98 - 109.\n",
    "\n",
    "To access the model, we download it directly using `urllib` and `numpy`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dbcb982b",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from urllib.request import urlopen\n",
    "\n",
    "raw_data = urlopen('http://pc213fis.fis.ucm.es/sha.dif.14k/coeff_SHA.DIF.14k.dat')\n",
    "data = np.genfromtxt(raw_data).T"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "837cd058",
   "metadata": {},
   "source": [
    "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",
   "execution_count": null,
   "id": "f4b5705d",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pymagglobal import Model\n",
    "from pymagglobal.utils import i2lm_l, i2lm_m, lmax2N\n",
    "\n",
    "from scipy.interpolate import BSpline\n",
    "\n",
    "\n",
    "# Inherit from pymagglobals Model-class\n",
    "class SHADIF14k(Model):\n",
    "    # Override the __init__ method to set up\n",
    "    # the model with our data\n",
    "    def __init__(self):\n",
    "        # set up the name\n",
    "        self.name = 'sha.dif.14k'\n",
    "        # The first column consists of the knots,\n",
    "        # but it contains many duplications\n",
    "        self.knots = np.unique(data[0])\n",
    "        # Set the interval in which the model is valid\n",
    "        self.t_min = -12000\n",
    "        self.t_max = 1800\n",
    "        # The first column contains the degrees, so\n",
    "        # we extract the maxum\n",
    "        self.l_max = int(np.max(data[1]))\n",
    "        # Now we go through the data and access the\n",
    "        # coefficients line by line\n",
    "        \n",
    "        # First set up an empty array of the right size\n",
    "        self.coeffs = np.empty((len(self.knots), lmax2N(self.l_max)))\n",
    "        # Then iterate over all coefficients\n",
    "        for it in range(lmax2N(self.l_max)):\n",
    "            # Get degree and order from the index\n",
    "            ell = i2lm_l(it)\n",
    "            m = i2lm_m(it)\n",
    "            # Get the corresponding indices in the data-array\n",
    "            inds = (data[1] == ell) * (data[2] == abs(m))\n",
    "            # Write the data to the coefficient-array\n",
    "            if m < 0:\n",
    "                # scale to nT\n",
    "                self.coeffs[:, it] = data[4][inds]*1000\n",
    "            else:\n",
    "                # scale to nT\n",
    "                self.coeffs[:, it] = data[3][inds]*1000\n",
    "\n",
    "        # Finally initialize the BSpline object\n",
    "        self.splines = BSpline(self.knots,\n",
    "                               self.coeffs,\n",
    "                               3)\n",
    "        # The model file does not report uncertainties\n",
    "        self.cov_splines = None"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4a463586",
   "metadata": {},
   "source": [
    "This model can now be initialized and used like any other model in `pymagglobal`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d891c8d5",
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "from pymagglobal import dipole_series\n",
    "\n",
    "plt.rcParams['font.size'] = '14'\n",
    "\n",
    "shadif14k = SHADIF14k()\n",
    "\n",
    "times = np.linspace(-1000, 1500, 501)\n",
    "\n",
    "fig, ax = plt.subplots(1, 1, figsize=(10, 6))\n",
    "for color, model in zip(['C0', 'C1'], [shadif14k, Model('ARCH10k.1')]):\n",
    "    dip = dipole_series(times, model)\n",
    "    ax.plot(times, dip, color=color, zorder=2, label=model.name)\n",
    "\n",
    "ax.legend();\n",
    "ax.set_ylabel(r'Dipole moment [$10^{22}$Am$^2$]');\n",
    "ax.set_xlabel(r'time [yrs.]');"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.9.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}