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

Add conversion between geodetic and geocentric systems. This changes the default behavior.

parent 87c7c00a
......@@ -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:
from pyfield import CRUSpline as Kernel
class TestAgainstPyfield(unittest.TestCase):
def test_evaluation(self):
......
Markdown is supported
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