Package Documentation

python interface for global geomagnetic field models (pymagglobal)

pymagglobal serves the purpose of replacing some Fortran scripts, which are used in the geomagnetism community to evaluate global field models. It can be applied to all cubic-spline based geomagnetic field models stored in the same file format as gufm1 or the CALSxk model series.

The package exposes several methods that can be applied to models or splines representing a model. To acces the models, use the Model class or its attribute Model.splines. With the models dictionary, several built-in models can be accessed, see also the list of included models.

class pymagglobal.Model(name, fname=None)[source]

Class for representing a global magnetic field model.

Parameters:

  • name (str) – The name of the model. If fname is none, the built-in model will be used.

  • fname (str, optional) – The path leading to a model file in the format of gufm1 and the CALSxk- series.

t_min, t_max[source]

The lower and upper bounds for the model time series. Above and below these the time series will be extrapolated and not robust.

Type:

float

l_max[source]

The maximal spherical harmonics degree included in the model.

Type:

int

knots[source]

The epochs of the spline knots.

Type:

array

coeffs[source]

A large array containing coefficients at the (inner) knots.

Type:

array

splines[source]

A BSpline instance representing the model.

Type:

scipy.interpolate.BSpline

cov_splines[source]

A BSPline instance giving the covariances. Only available for certain models. Default is None.

Type:

scipy.interpolate.BSpline

valid_degrees_orders(degrees, orders)[source]

Check whether the given degrees and orders are valid for the model and calculate the corresponding indices.

Parameters:

  • degrees (int or array-like of ints) – The degrees.

  • orders (int or array-like of ints) – The orders.

Returns:

The standard order indices corresponding to the given degrees and orders.

Return type:

list of ints

valid_epoch(epoch)[source]

Check whether the given epoch is in the range of the model.

Parameters:

epoch (float) – The epoch in years.

Returns:

The epoch in years.

Return type:

float

pymagglobal.coefficients(epoch, splines, cov_splines=None, check=True)[source]

Evaluate splines at an epoch and return the coefficients at the Earth’s surface, together with the appropriate degrees and orders. If coefficients at level different from the Earth’s surface are required, use utils.scaling to get the respective scaling factors.

Parameters:

  • epoch (float or array) – The epoch(s) at which to return the coefficients, in years.

  • splines (scipy.interpolate.BSpline or Model) – An instance of Model or splines specifying the model.

  • cov_splines (scipy.interpolate.BSpline, optional) – A BSpline interpolating the covariance matrices of the coefficients. If a model is passed to splines, this argument will be ignored and set, if the model’s covariances are available. If cov_splines is given, the standard deviations are returned as a fourth return value.

  • check (bool, optional) – If a Model is given, check the input validity.

Returns:

  • list – The list of coefficient degrees

  • list – The list of coefficient orders

  • ndarray – Coefficients of the model at the Earth’s surface for the given epoch, corresponding to the lists of degrees and orders. I.e. if you have

    l, m, coeffs = coefficients(…)

    in your code, coeffs[i] corresponds to degree l[i] and order m[i]. If you want to access a specific degree and order, you can use utils.lm2i like so:

    g_2_minus1 = coeffs[lm2i(2, -1)]

pymagglobal.dipole_series(times, splines, cov_splines=None, check=True)[source]

Create a dipole-moment time series from a splines object.

Parameters:

  • times (array-like) – The times for which to create the dipole series.

  • splines (scipy.interpolate.BSpline or Model) – An instance of Model or splines specifying the model.

  • cov_splines (scipy.interpolate.BSpline, optional) – A BSpline interpolating the covariance matrices of the coefficients. If a model is passed to splines, this argument will be ignored and set, if the model’s covariances are available. If cov_splines is given, the 16- and 84-percentiles are returned along with the series.

  • check (bool, optional) – If a Model is given, check the input validity.

Returns:

An array containing the dipole-moment time series. If cov_splines is given, the 16- and 84-percentiles are returned along with the series. All values are wrt. the Earth’s surface.

Return type:

ndarray

pymagglobal.field(grid, splines, cov_splines=None, field_type='nez', check=True, inp_gd=True, out_gd=True)[source]

Evaluate coefficient splines at locations grid and return the field. This routine is quite flexible and can be used to calculate the field at different radii, e.g. at the core-mantle boundary or in space or with the appropriate input array at a mixture of levels.

Parameters:

  • grid (array-like of shape (4, n) in Fortran order) –

    The locations at which to evaluate the splines. Should be as expected by utils.dsh_basis, i.e.

    • grid[0] contains co-latitudes in degrees.

    • grid[1] contains longitudes in degrees.

    • grid[2] contains radii in km.

    • grid[3] contains dates in years.

    You can use pymagglobal.utils.get_grid to generate an array of input points.

  • splines (scipy.interpolate.BSpline or Model) – An instance of Model or splines specifying the model.

  • cov_splines (scipy.interpolate.BSpline, optional) – A BSpline interpolating the covariance matrices of the coefficients. If a model is passed to splines, this argument will be ignored and set, if the model’s covariances are available. If cov_splines is given, the standard deviations are returned as a second return value.

  • field_type ({'dif', 'nez'}) – The type of output. May be either ‘nez’ (default) for north, east and down component or ‘dif’ for declination, inclination and intensity.

  • check (bool, optional) – If a Model is given, check the input validity.

  • inp_gd (bool, optional) – Whether the inputs are given in a geodetic reference frame.

  • out_gd (bool, optional) – Whether the outputs should be given in a geodetic reference frame.

Returns:

An array containing the field components at the input locations. If cov_splines is given, the standard deviations are returned as a second return value of identical shape.

Return type:

array of shape (3, n)

pymagglobal.file2splines(fname)[source]

Read a coefficient file in the gufm1/CAlSxk format and return a scipy.interpolate.BSpline object representing the model.

pymagglobal.local_curve(times, loc, splines, cov_splines=None, field_type='dif', check=True)[source]

Create local curves from a splines object. A local curve consists of a time series of magnetic field values at a specific location at the Earth’s surface.

Parameters:

  • times (array-like) – The times for which to create the local curve.

  • loc (tuple) – lat, lon tuple of the location at which to create the local curve.

  • splines (scipy.interpolate.BSpline or Model) – An instance of Model or splines specifying the model.

  • field_type ({'dif', 'nez'}) – The type of the local curves. May be either ‘dif’ (default) for declination, inclination and intensity or ‘nez’ for north, east, down.

  • cov_splines (scipy.interpolate.BSpline, optional) – A BSpline interpolating the covariance matrices of the coefficients. If a model is passed to splines, this argument will be ignored and set, if the model’s covariances are available.

  • check (bool, optional) – If a Model is given, check the input validity.

Returns:

  • ndarray – The first component local curve. Either declination or north, depending on the field_type kwarg. If cov_splines is given, a tuple containing the component and the standard deviation is returned.

  • ndarray – The second component local curve. Either inclination or east, depending on the field_type kwarg. If cov_splines is given, a tuple containing the component and the standard deviation is returned.

  • ndarray – The third component local curve. Either intensity or down, depending on the field_type kwarg. If cov_splines is given, a tuple containing the component and the standard deviation is returned.

pymagglobal.power_spectrum(epoch, splines, check=True, r=6371.2)[source]

Evaluate splines at an epoch and calculate the spatial power spectrum.

Parameters:

  • epoch (float or array) – The epoch(s) at which to return the coefficients, in years.

  • splines (scipy.interpolate.BSpline or Model) – An instance of Model or splines specifying the model.

  • check (bool, optional) – If a Model is given, check the input validity.

  • r (float, optional) – The radius at which to calculate the power spectrum.

Returns:

  • list – The list of coefficient degrees

  • array – The power spectrum over the given coefficients, with one power per degree.

pymagglobal.secular_variation(epoch, splines, cov_splines=None, check=True)[source]

Evaluate splines at an epoch and return the secular variation coefficients at the Earth’s surface, together with the appropriate degrees and orders.

Parameters:

  • epoch (float or array) – The epoch(s) at which to return the coefficients, in years.

  • splines (scipy.interpolate.BSpline or Model) – An instance of Model or splines specifying the model.

  • check (bool, optional) – If a Model is given, check the input validity.

Returns:

  • list – The list of coefficient degrees

  • list – The list of coefficient orders

  • ndarray – Secular variation coefficients of the model at the Earth’s surface for the given epoch, corresponding to the lists of degrees and orders.

Utilities

Utility functions for pymagglobal.

class pymagglobal.utils.WGS84[source]

A simple object to store WGS84 parameters.

pymagglobal.utils.dsh_basis(lmax, z, out=None, R=6371.2)[source]

Write the magnetic field basis functions, evaluated at points z, into the array out. These are basically the derivatives of the spherical harmonics in spherical coordinates with some scaling factors applied. Note that it is assumed that the coefficients that are multiplied to this basis are given at the radius R. This implementation is based on a specific recursion of the Legendre polynomials, to guarantee sane behavior at the poles. See [Du] for further details.

Parameters:

  • lmax (int) – The maximal spherical harmonics degree.

  • z (array of shape (3, n)) –

    The points at which to evaluate the basis functions, given as

    • z[0, :] contains the co-latitude in degrees.

    • z[1, :] contains the longitude in degrees.

    • z[2, :] contains the radius in kilometers.

  • out (array of shape (lmax2N(lmax), 3*n), optional) – The output array in which the basis functions are to be stored. This is included for historic reasons. Originally the function was to be replacable by a C-based function from the related pyfield library at a later stage [FieldTools].

  • R (float, optional) – The reference radius for the coefficients.

Returns:

out –

Array in which the basis functions are stored, as follows:
  • out[i, 0::3] contains the North component of the basis

    corresponding to degree l and order m, where i = lm2i(l, m).

  • out[i, 1::3] contains the East component of the basis.

  • out[i, 2::3] contains the Down component of the basis.

Return type:

array of shape (lmax2N(lmax), 3*n)

References

[Du]

J. Du, “Non-singular spherical harmonic expressions of geomagnetic vector and gradient tensor fields in the local north-oriented reference frame.”, Geosci. Model Dev., vol. 8, pages 1979-1990, 2014.

[FieldTools]

H. Matuschek and S. Mauerberger, “FieldTools - A toolbox for manipulating vector fields on the sphere”, GFZ Data Services, 2019. DOI: 10.5880/fidgeo.2019.033

pymagglobal.utils.equi_sph(n, twopi=True)[source]
Regularly place points on the surface of the unit sphere [Deserno]. the

number of output points may differ from n.

Parameters:

  • n (int) – The approximate number of points

  • twopi (bool, optional) – Whether to include 2*pi as a duplicate of 0 for the longitutinal angle (useful for plotting purposes).

Returns:

polar and azimutal angles in radians of m points, that are equally spaced on the sphere. m is approximately n.

Return type:

array of shape (2, m)

References

[Deserno]

M. Deserno, “How to generate equidistributed points on the surface of a sphere”, Max-Planck-Institut für Polymerforschung, 2004.

pymagglobal.utils.geodetic2geocentric(z, ellipsoid=<pymagglobal.utils.WGS84 object>)[source]

Assuming locations z are given in geodetic coordinates, transform them to geocentric ones.

Details are scattered in this article: https://en.wikipedia.org/wiki/Geographic_coordinate_conversion

Parameters:

  • z (array of shape (3, n)) –

    The points at which to evaluate the basis functions, given as

    • z[0, :] contains the geodetic co-latitude in degrees.

    • z[1, :] contains the geodetic longitude in degrees.

    • z[2, :] contains the geodetic altitude in kilometers, wrt. the

      Earth’s center.

  • ellipsoid (object, optional) – An object representing the reference ellipsoid. Should provide ellipsoid parameters a and f. Default is WGS84.

Returns:

z_gc – The transformed locations, in geocentric coordinates.

Return type:

array of shape (3, n)

pymagglobal.utils.get_grid(n, R=6371.2, t=1900.0, random=False, twopi=True)[source]
Get input points on the sphere. This is a convenience routine to allow

easy construction of field maps and synthetic data from the models.

Parameters:

  • n (int) – The approximate number of points. Due to the points being equally spaced on the sphere, the actual number may be slightly higher.

  • R (float, optional) – The radius of the sphere. By default this is the Earth’s radius.

  • t (float, optional) – The epoch at which the points are generated.

  • random (bool, optional) – If true, exactly n random points are returned. This is useful for generating synthetic data.

  • twopi (bool, optional) – Whether to include 2*pi as a duplicate of 0 for the longitutinal angle (useful for plotting purposes). Only used if random is False.

Returns:

grid –

  • grid[0] contains co-latitudes in degrees.

  • grid[1] contains longitudes in degrees.

  • grid[2] contains radii in km.

  • grid[3] contains dates in years.

n’ is approximately n.

Return type:

array of shape (4, n’)

pymagglobal.utils.grad_d(n, e, z)[source]

Calculate the gradient of D.

pymagglobal.utils.grad_i(n, e, z)[source]

Calculate the gradient of I.

pymagglobal.utils.i2lm_l(i)[source]

Returns the degree l of a Gauss coefficient at index i, in the “standard order”. The “standard order” is as follows

i

l

m

0

1

0

1

1

1

2

1

-1

3

2

0

4

2

1

5

2

-1

6

2

2

7

2

-2

…
Parameters:

i (int) – The index of a Gauss coefficient.

Returns:

The corresponding degree.

Return type:

int

Examples

>>> i2lm_l(0)
1
>>> i2lm_l(21)
4
pymagglobal.utils.i2lm_m(i)[source]

Returns the order m of a Gauss coefficient at index i, in the “standard order”. The “standard order” is as follows

i

l

m

0

1

0

1

1

1

2

1

-1

3

2

0

4

2

1

5

2

-1

6

2

2

7

2

-2

…
Parameters:

i (int) – The index of a Gauss coefficient.

Returns:

The corresponding order.

Return type:

int

Examples

>>> i2lm_m(0)
0
>>> i2lm_m(21)
-3
pymagglobal.utils.lm2i(ell, m)[source]

Returns the “standard order” index corresponding to a Gauss coefficient of degree l and order m. The “standard order” is as follows

i

l

m

0

1

0

1

1

1

2

1

-1

3

2

0

4

2

1

5

2

-1

6

2

2

7

2

-2

…
Parameters:

  • ell (int) – The degree of a Gauss coefficient.

  • m (int) – The order of a Gauss coefficient.

Returns:

The corresponding index.

Return type:

int

Examples

>>> lm2i(1, 0)
0
>>> lm2i(4, -1)
17
pymagglobal.utils.lmax2N(lmax)[source]

Returns the number of Gauss coefficients up to and including the degree lmax.

Parameters:

lmax (int) – The maximal spherical harmonic degree.

Returns:

The number of Gauss coefficients.

Return type:

int

pymagglobal.utils.lt2yr(times)[source]

Translate kilo-years before present into years CE.

Parameters:

times (float or int) – Kilo-years before present (backward counting from 1/1/1950).

Returns:

Years CE.

Return type:

float or int

Examples

>>> lt2yr(0)
1950
>>> lt2yr(2.)
-50.0
pymagglobal.utils.nez2dif(n, e, z)[source]

Transform the magnetic field components north, east and vertically downward to declination, inclination and intensity.

Parameters:

  • n – Pointing towards north.

  • e – Pointing towards east.

  • z – Pointing radially downwards.

Returns:

  • ndarray – Declination.

  • ndarray – Inclination.

  • ndarray – Intensity.

Example

>>> nez2dif(1, 0, 0)
(0.0, 0.0, 1.0)
pymagglobal.utils.scaling(r_from, r_to, lmax)[source]

Calculate the scaling matrix for Gauss-coefficients

pymagglobal.utils.yr2lt(times)[source]

Translate times given in years CE into kilo-years before present often used for longterm models.

Parameters:

times (float or int) – Years CE.

Returns:

Kilo years before present (backward counting from 1/1/1950).

Return type:

float

Examples

>>> yr2lt(1950)
0.0
>>> yr2lt(-50)
2.0