-
Notifications
You must be signed in to change notification settings - Fork 69
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Convert magnetic vector to inclination and declination (#402)
Create functions to convert the three components of magnetic vectors to the intensity, inclination a declination, and vice-versa. Add new functions to the API reference and add tests for them. --------- Co-authored-by: Lu Li <[email protected]> Co-authored-by: Santiago Soler <[email protected]>
- Loading branch information
1 parent
17cd09e
commit e885239
Showing
5 changed files
with
270 additions
and
1 deletion.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -125,4 +125,5 @@ Utilities | |
.. autosummary:: | ||
:toctree: generated/ | ||
|
||
test | ||
magnetic_vec_to_angles | ||
magnetic_angles_to_vec |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,150 @@ | ||
# Copyright (c) 2018 The Harmonica Developers. | ||
# Distributed under the terms of the BSD 3-Clause License. | ||
# SPDX-License-Identifier: BSD-3-Clause | ||
# | ||
# This code is part of the Fatiando a Terra project (https://www.fatiando.org) | ||
# | ||
import numpy as np | ||
|
||
|
||
def magnetic_angles_to_vec(intensity, inclination, declination): | ||
""" | ||
Convert magnetic field angles to magnetic field vector | ||
Convert intensity, inclination and declination angles of the magnetic field | ||
to a 3-component magnetic vector. | ||
.. note:: | ||
Inclination is measured positive downward from the horizontal plane, | ||
and declination is measured with respect to north, where positive | ||
angles indicate east. | ||
Parameters | ||
---------- | ||
intensity: float or array | ||
Intensity (norm) of the magnetic vector in A/m. | ||
inclination : float or array | ||
Inclination angle of the magnetic vector in degree. | ||
It must be in ``degrees``. | ||
declination : float or array | ||
Declination angle of the magnetic vector. | ||
It must be in ``degrees``. | ||
Returns | ||
------- | ||
magnetic_e : float or array | ||
Easting component of the magnetic vector. | ||
magnetic_n : float or array | ||
Northing component of the magnetic vector. | ||
magnetic_u : float or array | ||
Upward component of the magnetic vector. | ||
Examples | ||
-------- | ||
>>> mag_e, mag_n, mag_u = magnetic_angles_to_vec(3.0, 45.0, 45.0) | ||
>>> print(mag_e, mag_n, mag_u) | ||
1.5 1.5 -2.121 | ||
""" | ||
# Transform to radians | ||
inc_rad = np.radians(inclination) | ||
dec_rad = np.radians(declination) | ||
# Calculate the 3 components | ||
magnetic_e = intensity * np.cos(inc_rad) * np.sin(dec_rad) | ||
magnetic_n = intensity * np.cos(inc_rad) * np.cos(dec_rad) | ||
magnetic_u = -intensity * np.sin(inc_rad) | ||
return magnetic_e, magnetic_n, magnetic_u | ||
|
||
|
||
def magnetic_vec_to_angles(magnetic_e, magnetic_n, magnetic_u, degrees=True): | ||
r""" | ||
Convert magnetic field vector to magnetic field angles | ||
Convert the 3-component magnetic vector to intensity, and inclination and | ||
declination angles. | ||
.. note:: | ||
Inclination is measured positive downward from the horizontal plane and | ||
declination is measured with respect to North and it is positive east. | ||
Parameters | ||
---------- | ||
magnetic_e : float or array | ||
Easting component of the magnetic vector. | ||
magnetic_n : float or array | ||
Northing component of the magnetic vector. | ||
magnetic_u : float or array | ||
Upward component of the magnetic vector. | ||
degrees : bool (optional) | ||
If True, the angles are returned in degrees. | ||
If False, the angles are returned in radians. | ||
Default True. | ||
Returns | ||
------- | ||
intensity: float or array | ||
Intensity of the magnetic vector. | ||
inclination : float or array | ||
Inclination angle of the magnetic vector. | ||
If ``degrees`` is True, then the angle is returned in degree, else it's | ||
returned in radians. | ||
declination : float or array | ||
Declination angle of the magnetic vector. | ||
If ``degrees`` is True, then the angle is returned in degrees, else | ||
it's returned in radians. | ||
Notes | ||
----- | ||
The intensity of the magnetic vector is calculated as: | ||
.. math:: | ||
T = \sqrt{B_e^2 + B_n^2 + B_u^2} | ||
where :math:`B_e`, :math:`B_n`, :math:`B_u` are the easting, northing and | ||
upward components of the magnetic vector, respectively. | ||
The inclination angle is defined as the angle between the magnetic field | ||
vector and the horizontal plane: | ||
.. math:: | ||
I = \arctan \frac{- B_u}{\sqrt{B_e^2 + B_n^2}} | ||
And the declination angle is defined as the azimuth of the projection of | ||
the magnetic field vector onto the horizontal plane (starting from the | ||
northing direction, positive to the east and negative to the west): | ||
.. math:: | ||
D = \arcsin \frac{B_e}{\sqrt{B_e^2 + B_n^2}} | ||
Examples | ||
-------- | ||
>>> intensity, inc, dec = magnetic_vec_to_angles(1.5, 1.5, -2.12132) | ||
>>> print(intensity, inc, dec) | ||
3.0 45.0 45.0 | ||
""" | ||
# Compute the intensity as a norm | ||
intensity = np.sqrt(magnetic_e**2 + magnetic_n**2 + magnetic_u**2) | ||
# Compute the horizontal component of the magnetic vector | ||
horizontal_component = np.atleast_1d(np.sqrt(magnetic_e**2 + magnetic_n**2)) | ||
# Mask the values equal to zero in the horizontal component | ||
horizontal_component = np.ma.masked_values(horizontal_component, 0.0) | ||
# Calculate the inclination and declination using the mask | ||
inclination = np.arctan(-magnetic_u / horizontal_component) | ||
declination = np.arcsin(magnetic_e / horizontal_component) | ||
# Fill the masked values | ||
inclination = inclination.filled(-np.sign(magnetic_u) * np.pi / 2) | ||
declination = declination.filled(0) | ||
# Convert to degree if needed | ||
if degrees: | ||
inclination = np.degrees(inclination) | ||
declination = np.degrees(declination) | ||
# Cast to floats if all components are floats | ||
if intensity.ndim == 0: | ||
(inclination,) = inclination | ||
(declination,) = declination | ||
return intensity, inclination, declination |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,116 @@ | ||
# Copyright (c) 2018 The Harmonica Developers. | ||
# Distributed under the terms of the BSD 3-Clause License. | ||
# SPDX-License-Identifier: BSD-3-Clause | ||
# | ||
# This code is part of the Fatiando a Terra project (https://www.fatiando.org) | ||
# | ||
import numpy as np | ||
import numpy.testing as npt | ||
import pytest | ||
|
||
from .. import magnetic_angles_to_vec, magnetic_vec_to_angles | ||
|
||
VECTORS = [ | ||
[0.5, 0.5, -0.70710678], | ||
[0.5, 0.5, 0.70710678], | ||
[-0.5, 0.5, -0.70710678], | ||
[0, 0, -1], # Over -z axis | ||
[1, 0, 0], # Over east (y) axis | ||
[0, 1, 0], # Over north (x) axis | ||
] | ||
|
||
ANGLES = [[1, 45, 45], [1, -45, 45.0], [1, 45, -45], [1, 90, 0], [1, 0, 90], [1, 0, 0]] | ||
|
||
|
||
@pytest.mark.parametrize("angles, vector", [(a, v) for a, v in zip(ANGLES, VECTORS)]) | ||
def test_magnetic_ang_to_vec_float(angles, vector): | ||
""" | ||
Check if the function returns the expected values for a given intensity | ||
inclination and declination as float | ||
""" | ||
intensity, inclination, declination = angles | ||
magnetic_e, magnetic_n, magnetic_u = vector | ||
npt.assert_almost_equal( | ||
magnetic_angles_to_vec(intensity, inclination, declination), | ||
(magnetic_e, magnetic_n, magnetic_u), | ||
) | ||
|
||
|
||
@pytest.mark.parametrize("degrees", (False, True), ids=("radians", "degrees")) | ||
@pytest.mark.parametrize("angles, vector", [(a, v) for a, v in zip(ANGLES, VECTORS)]) | ||
def test_magnetic_vec_to_angles_float(angles, vector, degrees): | ||
""" | ||
Check if the function returns the expected values for a given magnetic | ||
vector as float | ||
""" | ||
intensity, inclination, declination = angles | ||
magnetic_e, magnetic_n, magnetic_u = vector | ||
if not degrees: | ||
inclination, declination = np.radians(inclination), np.radians(declination) | ||
npt.assert_allclose( | ||
magnetic_vec_to_angles(magnetic_e, magnetic_n, magnetic_u, degrees=degrees), | ||
(intensity, inclination, declination), | ||
) | ||
|
||
|
||
@pytest.fixture(name="arrays", params=["single-element", "multi-element"]) | ||
def angles_vectors_as_arrays(request): | ||
""" | ||
Generate magnetic angles and vectors as arrays | ||
""" | ||
if request.param == "single-element": | ||
# Generate arrays with a single element | ||
intensity, inclination, declination = tuple(np.atleast_1d(i) for i in ANGLES[0]) | ||
magnetic_e, magnetic_n, magnetic_u = tuple(np.atleast_1d(i) for i in VECTORS[0]) | ||
else: | ||
intensity, inclination, declination = np.vstack(ANGLES).T | ||
magnetic_e, magnetic_n, magnetic_u = np.vstack(VECTORS).T | ||
return (intensity, inclination, declination), (magnetic_e, magnetic_n, magnetic_u) | ||
|
||
|
||
def test_magnetic_ang_to_vec_array(arrays): | ||
""" | ||
Check if the function returns the expected values for a given intensity, | ||
inclination and declination a array | ||
""" | ||
intensity, inclination, declination = arrays[0] | ||
magnetic_e, magnetic_n, magnetic_u = arrays[1] | ||
npt.assert_almost_equal( | ||
magnetic_angles_to_vec(intensity, inclination, declination), | ||
(magnetic_e, magnetic_n, magnetic_u), | ||
) | ||
|
||
|
||
@pytest.mark.parametrize("degrees", (False, True), ids=("radians", "degrees")) | ||
def test_magnetic_vec_to_angles_array(arrays, degrees): | ||
""" | ||
Check if the function returns the expected values for the given magnetic | ||
vector as arrays | ||
""" | ||
intensity, inclination, declination = arrays[0] | ||
magnetic_e, magnetic_n, magnetic_u = arrays[1] | ||
if not degrees: | ||
inclination, declination = np.radians(inclination), np.radians(declination) | ||
npt.assert_allclose( | ||
magnetic_vec_to_angles(magnetic_e, magnetic_n, magnetic_u, degrees=degrees), | ||
(intensity, inclination, declination), | ||
) | ||
|
||
|
||
@pytest.mark.parametrize("start_with", ("angles", "vectors")) | ||
def test_identity(arrays, start_with): | ||
""" | ||
Check if applying both conversions return the original set of vectors | ||
""" | ||
if start_with == "angles": | ||
intensity, inclination, declination = arrays[0] | ||
vector = magnetic_angles_to_vec(intensity, inclination, declination) | ||
npt.assert_almost_equal( | ||
magnetic_vec_to_angles(*vector), (intensity, inclination, declination) | ||
) | ||
else: | ||
magnetic_e, magnetic_n, magnetic_u = arrays[1] | ||
angles = magnetic_vec_to_angles(magnetic_e, magnetic_n, magnetic_u) | ||
npt.assert_almost_equal( | ||
magnetic_angles_to_vec(*angles), (magnetic_e, magnetic_n, magnetic_u) | ||
) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters