Commit 1b8f1456 authored by Maximilian Schanner's avatar Maximilian Schanner
Browse files

Merge branch 'geodetic' into 'master'

Geodetic

See merge request !13
parents a2153a5a 95896eb8
Changelog
=========
2021-08-25 v1.1.0
-----------------
* Add geodetic->geocentric conversion for coordinates
* Add geocentric->geodetic conversion for fields
* Update docs
* Add some license files that were missing
2021-08-11 v1.0.1
-----------------
* Fix the range widget (return closer value if out of bounds; deal with string input)
......
This file complements pymagglobal.
It is licensed under CC0 1.0 (https://creativecommons.org/publicdomain/zero/1.0/)
......@@ -4,18 +4,9 @@
[[_TOC_]]
`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. However, care has to be taken that two header lines of the model files are formatted correctly and the list of spline knot point epochs starts only in line 3. The first header line has to contain start and end epoch of the model as the first two numbers, any further information in that line is ignored. The second header line has to start with three integers, which are the maximum spherical harmonic degree, a dummy that actually is not used, and the number of splines.
By default, `pymagglobal` includes several models. Use
```console
$ pymagglobal --list-models
```
to get a list of these default models or go to [pymagglobal/dat](https://git.gfz-potsdam.de/sec23/korte/pymagglobal/-/tree/master/pymagglobal/dat) for further information. Using
```console
$ pymagglobal ... <path/to/your_model>
```
you can use `pymagglobal` to evaluate your own models, if they come in a similar format. `<path/to/your_model>` specifies the path to your model and is given instead of the name of an included model. You can download additional models [here](ftp://ftp.gfz-potsdam.de/home/mag/arthus/pymagglobal_models/) and use them as above.
`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.
Once installed, `pymagglobal` can be imported and its routines used to access the models from inside your own python code.
For extensive documentation, including a more detailed description of the standard file format, check out the [docs](https://sec23.git-pages.gfz-potsdam.de/korte/pymagglobal).
## License
GNU General Public License, Version 3, 29 June 2007
......@@ -54,62 +45,6 @@ $ pip3 install pymagglobal --extra-index-url https://public:5mz_iyigu-WE3HySBH1J
Since [conda](https://docs.conda.io/) version 4.6, conda and pip get along well. So you can also run `pip3 install ...` from inside your conda environment.
## Documentation
Check out the extended documention [here](https://sec23.git-pages.gfz-potsdam.de/korte/pymagglobal).
pymagglobal comes with a GUI, that can be started from the command line via
```console
$ pymagglobal-gui
```
You can also use pymagglobal directly from the command line, to get various results from the models. For example,
```console
$ pymagglobal dipole gufm1
```
will give a plot of the dipole moment time series for the model gufm1. In general, pymagglobal is called as
```console
$ pymagglobal command --options model
```
where `command` specifies the quantity you want to get from `pymagglobal` and `model` is the respective model. You can use
```console
$ pymagglobal command --options <path/to/your_model>
```
to parse your own model, if it is in a format similar to gufm1. Use
```console
$ pymagglobal --help
```
to get further information. Each command has its own help, so you may also use
```console
$ pymagglobal dipole --help
```
to get information on the options for the dipole time series.
When using `python` you can import the pymagglobal package and access the models directly:
```python
import pymagglobal
```
We can first use `built_in_models`, to access a dictionary of available models:
```python
models = built_in_models()
```
pymagglobal provides a `Model` class. Built-in models can be accessed directly, custom models are set up with a name and a path:
```python
gufm1 = pymagglobal.Model('gufm1')
my_model = pymagglobal.Model('My model', '<path/to/my_model.dat>')
```
The model can be passed to routines in pymagglobal. For example, to get the dipole series from above use
```python
import numpy as np
times = np.linspace(1590, 1990, 201)
gufm1_dipoles = pymagglobal.dipole_series(times, gufm1)
```
Additionally, the object contains several quantities of interest, for example the minimal and maximal time for which the model is valid
```python
>>> gufm1.t_min
1590.0
>>> gufm1.t_max
1990.0
```
## Testing
To test your `pymagglobal` installation, run
......@@ -123,6 +58,9 @@ $ pip install pymagglobal[tests] --extra-index-url https://public:5mz_iyigu-WE3H
[FieldTools]: http://doi.org/10.5880/fidgeo.2019.033
## Documentation
Extensive documentation is available [here]((https://sec23.git-pages.gfz-potsdam.de/korte/pymagglobal).
## Contact
* [Maximilian Schanner](mailto:arthus@gfz-potsdam.de)
Helmholtz Centre Potsdam German Research Centre for Geoscienes GFZ
......
Command line interface
======================
The python package `pymagglobal` also offers a utility command
`pymagglobal` can be accessed directly via the command line. This is documented below.
.. argparse::
:module: pymagglobal.__main__
......
This file complements pymagglobal.
It is licensed under CC0 1.0 (https://creativecommons.org/publicdomain/zero/1.0/)
......@@ -2,23 +2,17 @@ A python interface for global geomagnetic field models
======================================================
.. toctree::
:hidden:
overview
installation
package_documentation
command_line_interface
examples
CHANGELOG
:code:`pymagglobal` serves the purpose of replacing some Fortran scripts, which are used in the geomagnetism community to evaluate global field models.
:ref:`pymagglobal <Overview>` 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.
:code:`pymagglobal` comes with a GUI, that can be started from the command line via
.. code-block:: bash
$ pymagglobal-gui
Installation
------------
.. note:
......@@ -26,16 +20,13 @@ Installation
This should also help if you receive :code:`ImportError: NumPy 1.10+ is required to install cartopy.`
:code:`pymagglobal` is distributed via the PyPI registry of the corresponding repository. It can be installed using
pymagglobal is distributed via the PyPI registry of the corresponding repository. It can be installed using
.. code-block:: bash
$ pip3 install pymagglobal --extra-index-url https://public:5mz_iyigu-WE3HySBH1J@git.gfz-potsdam.de/api/v4/projects/1055/packages/pypi/simple
See also `here <https://git.gfz-potsdam.de/sec23/korte/pymagglobal#installation>`_ and `here <https://git.gfz-potsdam.de/sec23/korte/pymagglobal#testing>`_.
.. include:: ../pymagglobal/dat/README.rst
Additional information can be found in the :ref:`installation section <Installation>`.
License
-------
......@@ -55,11 +46,3 @@ GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see <https://www.gnu.org/licenses/>.
Indices and tables
==================
* :ref:`genindex`
* :ref:`modindex`
* :ref:`search`
* `Repository <https://git.gfz-potsdam.de/sec23/korte/pymagglobal/>`_
.. _Installation:
Installation
============
pymagglobal is distributed via the PyPI registry of the corresponding repository. It can be installed using
.. code-block:: bash
$ pip3 install pymagglobal --extra-index-url https://public:5mz_iyigu-WE3HySBH1J@git.gfz-potsdam.de/api/v4/projects/1055/packages/pypi/simple
Dependencies
------------
pymagglobal depends on
* `numpy <https://numpy.org/>`_
* `scipy <https://www.scipy.org/>`_
* `matplotlib <https://matplotlib.org/>`_
* `cartopy <https://scitools.org.uk/cartopy/>`_
* `PyQt5 <https://www.riverbankcomputing.com/software/pyqt/>`_
Depending on your setup, the latter two may require additional work to install.
The most fluent way to install pymagglobal is to first `install cartopy <https://scitools.org.uk/cartopy/docs/latest/installing.html#installing>`_ and `PyQt5 <https://www.riverbankcomputing.com/static/Docs/PyQt5/installation.html>`__ and then run the pip command given above.
In general you should follow the installation instructions given on the respective package homepages.
Below we illuminate some ways to install the dependencies, that we pursued during troubleshooting.
If you run into problems during installing pymagglobal, do not hesitate to `contact us <mailto:arthus@gfz-potsdam.de>`_.
Please try to install the dependencies beforehand.
Conda
-----
Since `conda <https://docs.conda.io/>`_ version 4.6, conda and pip get along well. So you can also run `pip3 install ...` from inside your conda environment. The dependencies come via
.. code-block:: bash
conda install cartopy pyqt5
Debian
------
First you have to install python and pip:
.. code-block:: bash
apt-get install python3-dev python3-pip
The debian-packages necessary to install cartopy can be installed via
.. code-block:: bash
apt-get install libproj-dev libgeos-dev proj-bin
Afterwards cartopy will be installed alongside pymagglobal, if you run the pip command from above.
This file complements pymagglobal.
It is licensed under CC0 1.0 (https://creativecommons.org/publicdomain/zero/1.0/)
.. _Overview:
Overview
========
:code:`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.
A detailed description of the file format, together with a list of the included models can be found below.
:code:`pymagglobal` comes with a GUI, that can be started from the command line via
.. code-block:: bash
$ pymagglobal-gui
You can also use pymagglobal directly from the command line, to get various results from the models. For example,
.. code-block:: bash
$ pymagglobal dipole gufm1
will give a plot of the dipole moment time series for the model gufm1. In general, pymagglobal is called as
.. code-block:: bash
$ pymagglobal command --options model
where `command` specifies the quantity you want to get from `pymagglobal` and `model` is the respective model. You can use
.. code-block:: bash
$ pymagglobal command --options <path/to/your_model>
to parse your own model, if it is in a format similar to gufm1. Use
.. code-block:: bash
$ pymagglobal --help
to get further information. Each command has its own help, so you may also use
.. code-block:: bash
$ pymagglobal dipole --help
to get information on the options for the dipole time series.
When using `python` you can import the pymagglobal package and access the models directly:
.. code-block:: python
import pymagglobal
pymagglobal provides a `Model` class. Built-in models can be accessed directly, custom models are set up with a name and a path:
.. code-block:: python
gufm1 = pymagglobal.Model('gufm1')
my_model = pymagglobal.Model('My model', '<path/to/my_model.dat>')
The model can be passed to routines in pymagglobal. For example, to get the dipole series from above use
.. code-block:: python
import numpy as np
times = np.linspace(1590, 1990, 201)
gufm1_dipoles = pymagglobal.dipole_series(times, gufm1)
Additionally, the object contains several quantities of interest, for example the minimal and maximal time for which the model is valid
.. code-block:: python
>>> gufm1.t_min
1590.0
>>> gufm1.t_max
1990.0
By default, `pymagglobal` includes several models. A dictionary of the built in models is returned by `built_in_models`. The keys give a list of the model names, while the items point to the paths, where the models are stored. To access a dictionary of available models:
.. code-block:: python
models = built_in_models()
names = models.keys()
paths = models.items()
You can also use
.. code-block:: bash
$ pymagglobal --list-models
to get a list of these default models via the command line interface or check the table below.
.. code-block:: bash
$ pymagglobal ... <path/to/your_model>
you can use `pymagglobal` to evaluate your own models, if they come in a similar format. `<patodel>` specifies the path to your model and is given instead of the name of an included model.
.. include:: ../pymagglobal/dat/README.rst
This file complements pymagglobal.
It is licensed under CC0 1.0 (https://creativecommons.org/publicdomain/zero/1.0/)
......@@ -46,11 +46,4 @@ from pymagglobal.core import local_curve, dipole_series, file2splines, \
Model
from pymagglobal import utils
import warnings
# Monkey-patch the line away from warnings, as it is rather irritating.
formatwarning_orig = warnings.formatwarning
warnings.formatwarning = lambda msg, cat, fname, lineno, line=None: \
formatwarning_orig(msg, cat, fname, lineno, line='')
__version__ = '1.0.1'
__version__ = '1.1.0'
......@@ -27,7 +27,7 @@ import numpy as np
from scipy.interpolate import BSpline
from pymagglobal.utils import nez2dif, REARTH, i2lm_l, i2lm_m, lmax2N, lm2i, \
dsh_basis, grad_d, grad_i, grad_f
dsh_basis, grad_d, grad_i, grad_f, geodetic2geocentric, _rot_mat
__all__ = ['built_in_models', 'Model', 'local_curve', 'coefficients',
......@@ -177,7 +177,7 @@ class Model(object):
def local_curve(times, loc, splines, cov_splines=None, field_type='dif',
check=True):
check=True):
'''Create local curves from a splines object.
Parameters
......@@ -479,7 +479,8 @@ def file2splines(fname):
return c_splines
def field(z_at, splines, cov_splines=None, field_type='nez', check=True):
def field(z_at, splines, cov_splines=None, field_type='nez', check=True,
inp_gd=True, out_gd=True):
'''Evaluate coefficient splines at locations z_at and return the field.
Parameters
......@@ -507,6 +508,10 @@ def field(z_at, splines, cov_splines=None, field_type='nez', check=True):
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
-------
......@@ -522,6 +527,12 @@ def field(z_at, splines, cov_splines=None, field_type='nez', check=True):
if splines._has_covariances:
cov_splines = splines.cov_splines
splines = splines.splines
# if the inputs are given in a geodetic system, transform to geocentric
if inp_gd:
colats = np.copy(z_at[0])
z_at = geodetic2geocentric(z_at)
# count the number of inputs
n_pts = z_at.shape[1]
# count the number of coefficients
......@@ -540,6 +551,14 @@ def field(z_at, splines, cov_splines=None, field_type='nez', check=True):
coeffs,
spatial).T
if out_gd:
angs = z_at[0] - colats
rot_mats = _rot_mat(angs)
field = np.einsum('ij..., j...->i...',
rot_mats,
field)
if cov_splines is None:
if field_type == 'dif':
field = np.array(nez2dif(*field))
......@@ -547,12 +566,18 @@ def field(z_at, splines, cov_splines=None, field_type='nez', check=True):
covs = cov_splines(z_at[3]).T
if field_type == 'dif':
field_cov = np.einsum('ilk, ijl, jlm -> lkm',
spatial,
covs,
spatial)
field_cov = np.einsum('ilk, ijl, jlm -> lkm',
spatial,
covs,
spatial)
if out_gd:
field_cov = np.einsum('ji..., ...jk, kl... -> il...',
rot_mats,
field_cov,
rot_mats)
if field_type == 'dif':
grad_D = grad_d(*field)
cov_D = np.einsum('ij, ijk, ik -> i',
grad_D,
......@@ -574,9 +599,4 @@ def field(z_at, splines, cov_splines=None, field_type='nez', check=True):
return np.array(nez2dif(*field)), \
np.sqrt(np.array([cov_D, cov_I, cov_F]))
field_cov = np.einsum('ilk, ijl, jlk -> lk',
spatial,
covs,
spatial)
return field, np.sqrt(field_cov).T
......@@ -364,6 +364,109 @@ def get_z_at(n, R=REARTH, t=1900., random=False):
return z_at
class WGS84(object):
''' A simple object to store WGS84 parameters. '''
def __init__(self):
self.a = 6378.137
self.f = 1/298.257223563
def _geodetic2geocentric(lat, alt, ellipsoid=WGS84()):
''' From geodetic latitude and altitude calculate the geocentric latitude
and radius.
Parameters
----------
lat : float
The geodetic latitude in degrees.
alt : float
The geodetic altitude, in km over REARTH = 6371.2 km.
ellipsoid : object, optional
An object representing the reference ellipsoid. Should provide
ellipsoid parameters a and f. Default is WGS84.
Returns
-------
clat : float
The geocentric latitude in degrees.
rad : float
The geocentric radius wrt. to the Earth's center in km.
'''
# Geodetic parameters from WGS84
a = ellipsoid.a
f = ellipsoid.f
b = a * (1-f)
# print(a, b)
e2 = 1 - (b/a)**2
# e2 = f * (2 - f)
# The geodetic latitudes in radian
phi = np.deg2rad(lat)
# The altitudes
# The prime vertical radius of curvature
N = a / np.sqrt(1 - e2*np.sin(phi)**2)
# Horizontal length
H = (N + alt) * np.cos(phi)
# Z component
Z = (N * (1-f)**2 + alt) * np.sin(phi)
return np.rad2deg(np.arctan2(Z, H)), np.sqrt(H**2 + Z**2)
def _rot_mat(ang):
''' Given an array of angles ang, calculate 3d-rotation matrices around the
y-axis
'''
ang = np.atleast_1d(ang)
n = len(ang)
mat = np.repeat(np.eye(3), n).reshape(3, 3, n)
c = np.cos(np.deg2rad(ang))
s = np.sin(np.deg2rad(ang))
mat[0, 0, :] = c
mat[2, 2, :] = c
mat[2, 0, :] = s
mat[0, 2, :] = -s
if n == 1:
return mat[:, :, 0]
return mat
def geodetic2geocentric(z, ellipsoid=WGS84()):
''' 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 : array of shape (3, n)
The transformed locations, in geocentric coordinates.
'''
z_gc = np.copy(z)
z_gc[0], z_gc[2] = _geodetic2geocentric(90 - z[0], z[2] - REARTH)
# co-lat from lat
z_gc[0] = 90 - z_gc[0]
return z_gc
def dsh_basis(lmax, z, out, R=REARTH):
'''Write the magnetic field basis functions, evaluated at points z, into
the array out. These are basically the derivatives of the spherical
......
......@@ -21,7 +21,7 @@ import unittest
import numpy as np
from pymagglobal.utils import lm2i, i2lm_l, i2lm_m, lmax2N, REARTH, equi_sph, \
dsh_basis, get_z_at
dsh_basis, get_z_at, _geodetic2geocentric, _rot_mat, WGS84
def get_test_points(lmax=1, n=100):
......@@ -118,6 +118,48 @@ class TestBasisFunction(unittest.TestCase):
self.spatial)
class TestGeodeticGeocentric(unittest.TestCase):
def test_against_old_implementation(self):
''' Test against an old algorithm provided by M. Korte '''
ellipsoid = WGS84()
z = get_z_at(2, random=True)
a = ellipsoid.a
f = ellipsoid.f
a2 = a*a
b2 = (a * (1-f))**2
h = z[2] - REARTH
theta = np.pi/2-np.deg2rad(z[0])
clat = np.cos(theta)
slat = np.sin(theta)
one = a2*clat*clat
two = b2*slat*slat
three = one+two
four = np.sqrt(three)
r_ref = np.sqrt(h*(h+2.*four)+(a2*one+b2*two)/three)
cd = (h+four)/r_ref
sd = (a2-b2)/four*slat*clat/r_ref
sinth = slat*cd-clat*sd
costh = clat*cd+slat*sd
theta = np.arctan2(sinth, costh)
theta_ref = np.rad2deg(theta)
theta_gc, r = _geodetic2geocentric(90 - z[0], h)
self.assertTrue(np.allclose(theta_ref, theta_gc))
self.assertTrue(np.allclose(r_ref, r))
angs = 90 - theta_gc - z[0]
rot_mats = _rot_mat(angs)
self.assertTrue(np.allclose(rot_mats[0, 0], cd))
self.assertTrue(np.allclose(rot_mats[0, 2], -sd))
try:
from orthopoly import spherical_harmonic as sh_ref
......@@ -152,6 +194,7 @@ except ImportError:
try: