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]>

doc/api/index.rst

@@ -125,4 +125,5 @@ Utilities
 .. autosummary::
  :toctree: generated/
 
- test
+ magnetic_vec_to_angles
+ magnetic_angles_to_vec

harmonica/__init__.py

@@ -29,6 +29,7 @@
  reduction_to_pole,
  upward_continuation,
 )
+from ._utils import magnetic_angles_to_vec, magnetic_vec_to_angles
 from ._version import __version__
 
 

harmonica/_utils.py

@@ -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

harmonica/tests/test_utils.py

@@ -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)
+ )

pyproject.toml

@@ -10,6 +10,7 @@ local_scheme = "node-and-date"
 write_to = "harmonica/_version_generated.py"
 
 [tool.pytest.ini_options]
+doctest_optionflags = "NUMBER"
 markers = [
  "use_numba: mark test functions that call Numba jitted functions"
 ]