|
17 | 17 | "source": [ |
18 | 18 | "Here we demonstrate how to make a simple time-independent main field model by fitting spherical harmonics to ground observatory measurements.\n", |
19 | 19 | "\n", |
20 | | - "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", |
| 20 | + "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](https://chaosmagpy.readthedocs.io/)\n", |
21 | 21 | "\n", |
22 | | - "- [ChaosMagPy](https://chaosmagpy.readthedocs.io/), Clemens Kloss, https://doi.org/10.5281/zenodo.3352398 \n", |
23 | | - " We will use:\n", |
| 22 | + "> Clemens Kloss et al. (2025). ancklo/ChaosMagPy: ChaosMagPy v0.15 (v0.15). Zenodo. https://doi.org/10.5281/zenodo.16091680\n", |
| 23 | + "\n", |
| 24 | + "We will use:\n", |
24 | 25 | " - [design_gauss](https://chaosmagpy.readthedocs.io/en/master/functions/chaosmagpy.model_utils.design_gauss.html) to generate the design matrix for the inversion\n", |
25 | 26 | " - [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", |
26 | 27 | "\n", |
|
39 | 40 | "outputs": [], |
40 | 41 | "source": [ |
41 | 42 | "# Install pooch (not currently in VRE)\n", |
42 | | - "import sys\n", |
43 | | - "!{sys.executable} -m pip install pooch" |
| 43 | + "%pip install pooch" |
44 | 44 | ] |
45 | 45 | }, |
46 | 46 | { |
|
75 | 75 | "id": "638c7a53-8f91-4563-9dd0-66be049c7358", |
76 | 76 | "metadata": {}, |
77 | 77 | "source": [ |
78 | | - "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." |
| 78 | + "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." |
79 | 79 | ] |
80 | 80 | }, |
81 | 81 | { |
|
160 | 160 | "$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", |
161 | 161 | "\n", |
162 | 162 | "$$\n", |
163 | | - "\\begin{align}\\label{eq:spharm_B}\n", |
| 163 | + "\\begin{align}\n", |
164 | 164 | " \\vec{d} =\\ & G \\vec{m} \\\\\n", |
165 | 165 | " B_i =\\ & G \\begin{bmatrix}\n", |
166 | 166 | " g_n^m & h_n^m\n", |
|
197 | 197 | "metadata": {}, |
198 | 198 | "outputs": [], |
199 | 199 | "source": [ |
200 | | - "def build_model(ds=input_data, nmax=10):\n", |
201 | | - " \"\"\"Returns Gauss coefficients for a simple SH model\"\"\"\n", |
202 | | - " theta = ds[\"Colat\"]\n", |
203 | | - " phi = ds[\"Long\"]\n", |
204 | | - " radius = 6371.2 + ds[\"alt(m)\"]/1e3\n", |
205 | | - " B_radius = -ds[\"Z\"]\n", |
206 | | - " B_theta = -ds[\"X\"]\n", |
207 | | - " B_phi = ds[\"Y\"]\n", |
208 | | - " # Build the design matrix\n", |
209 | | - " G_radius, G_theta, G_phi = design_gauss(radius, theta, phi, nmax=nmax, source=\"internal\")\n", |
210 | | - " # Perform the inversion\n", |
211 | | - " G = np.concatenate((G_radius, G_theta, G_phi))\n", |
212 | | - " d = np.concatenate((B_radius, B_theta, B_phi))\n", |
213 | | - " m = np.linalg.inv(G.T @ G) @ (G.T @ d)\n", |
214 | | - " return m\n", |
| 200 | + "# Identify coordinates and measurements\n", |
| 201 | + "theta = input_data[\"Colat\"]\n", |
| 202 | + "phi = input_data[\"Long\"]\n", |
| 203 | + "radius = 6371.2 + input_data[\"alt(m)\"]/1e3\n", |
| 204 | + "B_radius = -input_data[\"Z\"]\n", |
| 205 | + "B_theta = -input_data[\"X\"]\n", |
| 206 | + "B_phi = input_data[\"Y\"]\n", |
215 | 207 | "\n", |
| 208 | + "# Build the design matrix\n", |
| 209 | + "N_max = 6\n", |
| 210 | + "G_radius, G_theta, G_phi = design_gauss(radius, theta, phi, nmax=N_max, source=\"internal\")\n", |
216 | 211 | "\n", |
217 | | - "m = build_model(input_data, nmax=5)\n", |
| 212 | + "# Perform the inversion\n", |
| 213 | + "G = np.concatenate((G_radius, G_theta, G_phi))\n", |
| 214 | + "d = np.concatenate((B_radius, B_theta, B_phi))\n", |
| 215 | + "m = np.linalg.inv(G.T @ G) @ (G.T @ d)\n", |
218 | 216 | "m" |
219 | 217 | ] |
220 | 218 | }, |
|
226 | 224 | "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." |
227 | 225 | ] |
228 | 226 | }, |
| 227 | + { |
| 228 | + "cell_type": "markdown", |
| 229 | + "id": "0d9b50ad-b8dd-408f-8183-fb2892f917c1", |
| 230 | + "metadata": {}, |
| 231 | + "source": [ |
| 232 | + "### Exercise" |
| 233 | + ] |
| 234 | + }, |
| 235 | + { |
| 236 | + "cell_type": "markdown", |
| 237 | + "id": "ce6a3485-7187-4a0d-814b-1147afdb24a0", |
| 238 | + "metadata": {}, |
| 239 | + "source": [ |
| 240 | + "1. What is the shape of the design matrix?\n", |
| 241 | + "1. What is the amplitude of the axial dipole ($g^1_0$) and how does it compare with the [IGRF](https://doi.org/10.5281/zenodo.14012302)?" |
| 242 | + ] |
| 243 | + }, |
229 | 244 | { |
230 | 245 | "cell_type": "markdown", |
231 | 246 | "id": "53a9d939-726c-40f2-a596-dc04409060a2", |
|
239 | 254 | "id": "312a3527-d2b4-4a11-b17d-84f9e9ec73bd", |
240 | 255 | "metadata": {}, |
241 | 256 | "source": [ |
242 | | - "Evaluating a model over the Earth and plotting contours..." |
| 257 | + "Let's evaluate the model over the Earth and plot contours..." |
243 | 258 | ] |
244 | 259 | }, |
245 | 260 | { |
|
249 | 264 | "metadata": {}, |
250 | 265 | "outputs": [], |
251 | 266 | "source": [ |
| 267 | + "def build_model(ds=input_data, nmax=10):\n", |
| 268 | + " \"\"\"Returns Gauss coefficients for a simple SH model\"\"\"\n", |
| 269 | + " theta = ds[\"Colat\"]\n", |
| 270 | + " phi = ds[\"Long\"]\n", |
| 271 | + " radius = 6371.2 + ds[\"alt(m)\"]/1e3\n", |
| 272 | + " B_radius = -ds[\"Z\"]\n", |
| 273 | + " B_theta = -ds[\"X\"]\n", |
| 274 | + " B_phi = ds[\"Y\"]\n", |
| 275 | + " # Build the design matrix\n", |
| 276 | + " G_radius, G_theta, G_phi = design_gauss(radius, theta, phi, nmax=nmax, source=\"internal\")\n", |
| 277 | + " # Perform the inversion\n", |
| 278 | + " G = np.concatenate((G_radius, G_theta, G_phi))\n", |
| 279 | + " d = np.concatenate((B_radius, B_theta, B_phi))\n", |
| 280 | + " m = np.linalg.inv(G.T @ G) @ (G.T @ d)\n", |
| 281 | + " return m\n", |
| 282 | + "\n", |
| 283 | + "\n", |
252 | 284 | "def sample_model_on_grid(m, radius=6371.200):\n", |
253 | 285 | " \"\"\"Evaluate Gauss coefficients over a regular grid over Earth\"\"\"\n", |
254 | 286 | " theta = np.arange(1, 180, 1)\n", |
|
272 | 304 | " )\n", |
273 | 305 | "\n", |
274 | 306 | "\n", |
| 307 | + "m = build_model(input_data, nmax=6)\n", |
275 | 308 | "plot_model_contours(m)" |
276 | 309 | ] |
277 | 310 | }, |
|
280 | 313 | "id": "153e2f1a-87e8-4dee-ac0e-baa5a6a52cf8", |
281 | 314 | "metadata": {}, |
282 | 315 | "source": [ |
283 | | - "The power spectrum is a common way to assess and compare models, showing how much power is concentrated at each spherical harmonic degree:" |
| 316 | + "The power spectrum is a common way to assess and compare models, showing how much power is concentrated at each spherical harmonic degree. We can use [`plot_power_spectrum`](https://chaosmagpy.readthedocs.io/en/master/functions/chaosmagpy.plot_utils.plot_power_spectrum.html#chaosmagpy.plot_utils.plot_power_spectrum) from ChaosMagPy:" |
284 | 317 | ] |
285 | 318 | }, |
286 | 319 | { |
|
347 | 380 | }, |
348 | 381 | { |
349 | 382 | "cell_type": "markdown", |
350 | | - "id": "8b20d365-844e-49d9-ab4a-2823307504a6", |
| 383 | + "id": "964e4fad-ce8e-4667-b482-48699848750b", |
351 | 384 | "metadata": {}, |
352 | 385 | "source": [ |
353 | | - "## Exercises\n", |
354 | | - "\n", |
| 386 | + "## Exercise" |
| 387 | + ] |
| 388 | + }, |
| 389 | + { |
| 390 | + "cell_type": "markdown", |
| 391 | + "id": "0bbf348f-91fb-4ca6-b82c-4b1afc7f061e", |
| 392 | + "metadata": {}, |
| 393 | + "source": [ |
| 394 | + "- Calculate the RMS misfit for each vector component (Hint: we have stored the residuals in the dataset, accessible as `ds[\"res_N\"]`, `ds[\"res_E\"]`, `ds[\"res_C\"]`)\n", |
355 | 395 | "- Experiment with changing the SH degree trunction in the model\n", |
356 | 396 | "- Experiment with using data just over Europe\n", |
357 | 397 | "- Add random noise into the data to observe the effect\n", |
|
0 commit comments