Commit b24e515f authored by Maximilian Schanner's avatar Maximilian Schanner
Browse files

Add generation for fixed coeffs that are no model.

parent a6c5bdeb
...@@ -19,7 +19,7 @@ ...@@ -19,7 +19,7 @@
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": null, "execution_count": null,
"id": "5380320c", "id": "f52a366b",
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
...@@ -30,7 +30,7 @@ ...@@ -30,7 +30,7 @@
}, },
{ {
"cell_type": "markdown", "cell_type": "markdown",
"id": "98302e49", "id": "13988168",
"metadata": {}, "metadata": {},
"source": [ "source": [
"Next we have to generate a data distribution. We use `pymagglobal`s `get_z_at` routine, to generate `n_at` random locations, that are uniformly distributed on the sphere. Times are drawn uniformly as well." "Next we have to generate a data distribution. We use `pymagglobal`s `get_z_at` routine, to generate `n_at` random locations, that are uniformly distributed on the sphere. Times are drawn uniformly as well."
...@@ -39,7 +39,7 @@ ...@@ -39,7 +39,7 @@
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": null, "execution_count": null,
"id": "68d9a6d8", "id": "6f57b72a",
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
...@@ -64,7 +64,7 @@ ...@@ -64,7 +64,7 @@
}, },
{ {
"cell_type": "markdown", "cell_type": "markdown",
"id": "cbf80490", "id": "806081dc",
"metadata": {}, "metadata": {},
"source": [ "source": [
"Finally we can evaluate the model at the inputs, to get synthetic ''observations'' of the field" "Finally we can evaluate the model at the inputs, to get synthetic ''observations'' of the field"
...@@ -73,7 +73,7 @@ ...@@ -73,7 +73,7 @@
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": null, "execution_count": null,
"id": "ea10fa42", "id": "5842b198",
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
...@@ -83,7 +83,7 @@ ...@@ -83,7 +83,7 @@
}, },
{ {
"cell_type": "markdown", "cell_type": "markdown",
"id": "0c28838f", "id": "14502d90",
"metadata": {}, "metadata": {},
"source": [ "source": [
"Let's have a look at the records:" "Let's have a look at the records:"
...@@ -92,7 +92,7 @@ ...@@ -92,7 +92,7 @@
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": null, "execution_count": null,
"id": "119a535b", "id": "67cd6ed7",
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
...@@ -108,7 +108,7 @@ ...@@ -108,7 +108,7 @@
}, },
{ {
"cell_type": "markdown", "cell_type": "markdown",
"id": "99e9b776", "id": "e96b4ec7",
"metadata": {}, "metadata": {},
"source": [ "source": [
"Declination, inclination and intensity records are also easily generated, by using the `field` kwargs:" "Declination, inclination and intensity records are also easily generated, by using the `field` kwargs:"
...@@ -117,7 +117,7 @@ ...@@ -117,7 +117,7 @@
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": null, "execution_count": null,
"id": "8bda1baf", "id": "57b38352",
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
...@@ -131,6 +131,116 @@ ...@@ -131,6 +131,116 @@
" axs[it].set_global()\n", " axs[it].set_global()\n",
" axs[it].coastlines();" " axs[it].coastlines();"
] ]
},
{
"cell_type": "markdown",
"id": "a16058b9",
"metadata": {},
"source": [
"## Data from an arbitray set of coefficients"
]
},
{
"cell_type": "markdown",
"id": "003f37e3",
"metadata": {},
"source": [
"Sometimes it is helpful to generate data not from a model that's included in pymagglobal, but from an arbitrary set of coefficients. Provided they are in the ''standard order'', i.e. $g_1^0, g_1^1, g_1^{-1}, g_2^0, g_2^1, g_2^{-1}, ...$, this is also straight-forward using utility routines from pymagglobal:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "fdbdfd1c",
"metadata": {},
"outputs": [],
"source": [
"from pymagglobal.core import coefficients\n",
"\n",
"epoch = -3000\n",
"# First generate coefficients at epoch\n",
"_, _, coeffs = coefficients(epoch, myModel)\n",
"\n",
"# Generate new input locations at the given epoch\n",
"z_at = get_z_at(n_at, random=True, t=epoch)"
]
},
{
"cell_type": "markdown",
"id": "e81a742d",
"metadata": {},
"source": [
"We use the pymagglobal routine `dsh_basis`, which evaluates the derivatives of the spherical harmonics up to a given degree at given inputs:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "5924dde3",
"metadata": {},
"outputs": [],
"source": [
"from pymagglobal.utils import dsh_basis\n",
"\n",
"# dsh_basis was designed to work with a C backend,\n",
"# so the return array has to be allocated beforehand\n",
"# and is filled during the function call\n",
"base = np.empty((len(coeffs), 3*z_at.shape[1]))\n",
"dsh_basis(myModel.l_max, z_at, base)"
]
},
{
"cell_type": "markdown",
"id": "0d875959",
"metadata": {},
"source": [
"Field values are then generated via a dot product:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a0a43e60",
"metadata": {},
"outputs": [],
"source": [
"obs = coeffs @ base"
]
},
{
"cell_type": "markdown",
"id": "1c804ae2",
"metadata": {},
"source": [
"The outputs are now given by obs, every third value is N, E, Z. To have a more convenient form, we reshape the output:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a6eb5bc3",
"metadata": {},
"outputs": [],
"source": [
"obs = obs.reshape(n_at, 3).T"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "d3ffc91a",
"metadata": {},
"outputs": [],
"source": [
"fig, axs = plt.subplots(1, 3, subplot_kw={'projection': ccrs.Mollweide()},\n",
" figsize=(15, 4))\n",
"titles = [r'$B_N$', r'$B_E$', r'$B_Z$']\n",
"for it in range(3):\n",
" axs[it].set_title(titles[it])\n",
" axs[it].scatter(z_at[1], 90-z_at[0], c=obs[it], transform=ccrs.PlateCarree())\n",
" axs[it].set_global()\n",
" axs[it].coastlines();"
]
} }
], ],
"metadata": { "metadata": {
......
%% Cell type:markdown id:681fb4de tags: %% Cell type:markdown id:681fb4de tags:
# Generate synthetic data # Generate synthetic data
%% Cell type:markdown id:7c938d6b tags: %% Cell type:markdown id:7c938d6b tags:
Here we briefly show how `pymagglobal` can be used to generate synthetic data. We first set up the model we want to use to generate the synthetic data: Here we briefly show how `pymagglobal` can be used to generate synthetic data. We first set up the model we want to use to generate the synthetic data:
%% Cell type:code id:5380320c tags: %% Cell type:code id:f52a366b tags:
``` python ``` python
from pymagglobal import Model from pymagglobal import Model
myModel = Model('CALS10k.2') myModel = Model('CALS10k.2')
``` ```
%% Cell type:markdown id:98302e49 tags: %% Cell type:markdown id:13988168 tags:
Next we have to generate a data distribution. We use `pymagglobal`s `get_z_at` routine, to generate `n_at` random locations, that are uniformly distributed on the sphere. Times are drawn uniformly as well. Next we have to generate a data distribution. We use `pymagglobal`s `get_z_at` routine, to generate `n_at` random locations, that are uniformly distributed on the sphere. Times are drawn uniformly as well.
%% Cell type:code id:68d9a6d8 tags: %% Cell type:code id:6f57b72a tags:
``` python ``` python
from pymagglobal.utils import get_z_at from pymagglobal.utils import get_z_at
import numpy as np import numpy as np
from matplotlib import pyplot as plt from matplotlib import pyplot as plt
from cartopy import crs as ccrs from cartopy import crs as ccrs
# the number of artificial records # the number of artificial records
n_at = 400 n_at = 400
z_at = get_z_at(n_at, random=True) z_at = get_z_at(n_at, random=True)
# set the times to something arbitrary # set the times to something arbitrary
z_at[3] = np.random.uniform(-1000, 1900, size=n_at) z_at[3] = np.random.uniform(-1000, 1900, size=n_at)
# plot the records, convert co-lat to lat # plot the records, convert co-lat to lat
fig, ax = plt.subplots(1, 1, subplot_kw={'projection': ccrs.Mollweide()}) fig, ax = plt.subplots(1, 1, subplot_kw={'projection': ccrs.Mollweide()})
ax.scatter(z_at[1], 90-z_at[0], transform=ccrs.PlateCarree()) ax.scatter(z_at[1], 90-z_at[0], transform=ccrs.PlateCarree())
ax.set_global() ax.set_global()
ax.coastlines(); ax.coastlines();
``` ```
%% Cell type:markdown id:cbf80490 tags: %% Cell type:markdown id:806081dc tags:
Finally we can evaluate the model at the inputs, to get synthetic ''observations'' of the field Finally we can evaluate the model at the inputs, to get synthetic ''observations'' of the field
%% Cell type:code id:ea10fa42 tags: %% Cell type:code id:5842b198 tags:
``` python ``` python
from pymagglobal.core import field from pymagglobal.core import field
obs = field(z_at, myModel) obs = field(z_at, myModel)
``` ```
%% Cell type:markdown id:0c28838f tags: %% Cell type:markdown id:14502d90 tags:
Let's have a look at the records: Let's have a look at the records:
%% Cell type:code id:119a535b tags: %% Cell type:code id:67cd6ed7 tags:
``` python ``` python
fig, axs = plt.subplots(1, 3, subplot_kw={'projection': ccrs.Mollweide()}, fig, axs = plt.subplots(1, 3, subplot_kw={'projection': ccrs.Mollweide()},
figsize=(15, 4)) figsize=(15, 4))
titles = [r'$B_N$', r'$B_E$', r'$B_Z$'] titles = [r'$B_N$', r'$B_E$', r'$B_Z$']
for it in range(3): for it in range(3):
axs[it].set_title(titles[it]) axs[it].set_title(titles[it])
axs[it].scatter(z_at[1], 90-z_at[0], c=obs[it], transform=ccrs.PlateCarree()) axs[it].scatter(z_at[1], 90-z_at[0], c=obs[it], transform=ccrs.PlateCarree())
axs[it].set_global() axs[it].set_global()
axs[it].coastlines(); axs[it].coastlines();
``` ```
%% Cell type:markdown id:99e9b776 tags: %% Cell type:markdown id:e96b4ec7 tags:
Declination, inclination and intensity records are also easily generated, by using the `field` kwargs: Declination, inclination and intensity records are also easily generated, by using the `field` kwargs:
%% Cell type:code id:8bda1baf tags: %% Cell type:code id:57b38352 tags:
``` python ``` python
obs_dif = field(z_at, myModel, field_type='dif') obs_dif = field(z_at, myModel, field_type='dif')
fig, axs = plt.subplots(1, 3, subplot_kw={'projection': ccrs.Mollweide()}, fig, axs = plt.subplots(1, 3, subplot_kw={'projection': ccrs.Mollweide()},
figsize=(15, 4)) figsize=(15, 4))
titles = [r'$D$', r'$I$', r'$F$'] titles = [r'$D$', r'$I$', r'$F$']
for it in range(3): for it in range(3):
axs[it].set_title(titles[it]) axs[it].set_title(titles[it])
axs[it].scatter(z_at[1], 90-z_at[0], c=obs_dif[it], transform=ccrs.PlateCarree()) axs[it].scatter(z_at[1], 90-z_at[0], c=obs_dif[it], transform=ccrs.PlateCarree())
axs[it].set_global() axs[it].set_global()
axs[it].coastlines(); axs[it].coastlines();
``` ```
%% Cell type:markdown id:a16058b9 tags:
## Data from an arbitray set of coefficients
%% Cell type:markdown id:003f37e3 tags:
Sometimes it is helpful to generate data not from a model that's included in pymagglobal, but from an arbitrary set of coefficients. Provided they are in the ''standard order'', i.e. $g_1^0, g_1^1, g_1^{-1}, g_2^0, g_2^1, g_2^{-1}, ...$, this is also straight-forward using utility routines from pymagglobal:
%% Cell type:code id:fdbdfd1c tags:
``` python
from pymagglobal.core import coefficients
epoch = -3000
# First generate coefficients at epoch
_, _, coeffs = coefficients(epoch, myModel)
# Generate new input locations at the given epoch
z_at = get_z_at(n_at, random=True, t=epoch)
```
%% Cell type:markdown id:e81a742d tags:
We use the pymagglobal routine `dsh_basis`, which evaluates the derivatives of the spherical harmonics up to a given degree at given inputs:
%% Cell type:code id:5924dde3 tags:
``` python
from pymagglobal.utils import dsh_basis
# dsh_basis was designed to work with a C backend,
# so the return array has to be allocated beforehand
# and is filled during the function call
base = np.empty((len(coeffs), 3*z_at.shape[1]))
dsh_basis(myModel.l_max, z_at, base)
```
%% Cell type:markdown id:0d875959 tags:
Field values are then generated via a dot product:
%% Cell type:code id:a0a43e60 tags:
``` python
obs = coeffs @ base
```
%% Cell type:markdown id:1c804ae2 tags:
The outputs are now given by obs, every third value is N, E, Z. To have a more convenient form, we reshape the output:
%% Cell type:code id:a6eb5bc3 tags:
``` python
obs = obs.reshape(n_at, 3).T
```
%% Cell type:code id:d3ffc91a tags:
``` python
fig, axs = plt.subplots(1, 3, subplot_kw={'projection': ccrs.Mollweide()},
figsize=(15, 4))
titles = [r'$B_N$', r'$B_E$', r'$B_Z$']
for it in range(3):
axs[it].set_title(titles[it])
axs[it].scatter(z_at[1], 90-z_at[0], c=obs[it], transform=ccrs.PlateCarree())
axs[it].set_global()
axs[it].coastlines();
```
......
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