Repository: fatiando/harmonica
Pull request #402 (2023-04-22 01:36:30)

Convert magnetic vector to inclination and declination
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.


santisoler: > Looks awesome! Just some minor suggestions for the docstring. Also, should we make a note somewhere, that our coordinate is upward up, so we have a negative magnetic vector (magnetic_u) with downward magnetization? Sometimes the upward is down, which might be a little confusing for people, haha.

Thanks for the review @LL-Geo! You're right that we should probably be more explicit about directions. For Harmonica I'm considering that both for the magnetic and magnetization vectors, the vertical component points **upwards**. That's why we are using `magnetic_u` and not `magnetic_z`.

> Sometimes the upward is down, which might be a little confusing for people, haha.

Sounds like an [Escher's painting](https://upload.wikimedia.org/wikipedia/en/a/a3/Escher%27s_Relativity.jpg) haha. Sometimes people have coordinate systems with vertical pointing down. We try to overcome that confusion by explicitly putting the word **upward**.

The tricky thing in these functions is that the inclination angle is defined as positive when the vector points downwards. So, the equations for transforming them are not exactly the ones in Blakely (1995): in there they use z->down.

LL-Geo: 
> Sounds like an [Escher's painting](https://upload.wikimedia.org/wikipedia/en/a/a3/Escher%27s_Relativity.jpg) haha. Sometimes people have coordinate systems with vertical pointing down. We try to overcome that confusion by explicitly putting the word **upward**.
> 
> The tricky thing in these functions is that the inclination angle is defined as positive when the vector points downwards. So, the equations for transforming them are not exactly the ones in Blakely (1995): in there they use z->down.

I like that painting! hahaha. That's the world in my dream🤣

Yeh... 'Traditional' coordinates usually point vertically down. And always confuse people. One day I suddenly realised that the positive density generates negative gravity in SimPEG due to the definition of "positive up and down". hahaha...


santisoler: > > Sounds like an [Escher's painting](https://upload.wikimedia.org/wikipedia/en/a/a3/Escher%27s_Relativity.jpg) haha. Sometimes people have coordinate systems with vertical pointing down. We try to overcome that confusion by explicitly putting the word **upward**.
> > The tricky thing in these functions is that the inclination angle is defined as positive when the vector points downwards. So, the equations for transforming them are not exactly the ones in Blakely (1995): in there they use z->down.
> 
> I like that painting! hahaha. That's the world in my dreamrofl
> 
> Yeh... 'Traditional' coordinates usually point vertically down. And always confuse people. One day I suddenly realised that the positive density generates negative gravity in SimPEG due to the definition of "positive up and down". hahaha...

Indeed! That's why we compute the `g_z` in Harmonica, but explicitly saying that it's the **downward** component, although our coordinate system has the vertical pointing upwards.
santisoler: @aguspesce I think this is ready to be merged after I took the liberty to make a few changes. Would you like to review it and share your thoughts?

Although it would be nice to have some example for these two functions, I think we can avoid including them here and just use these functions in other example that forward models some magnetic field. For example, the ones in #369 

Thanks a lot for the time you put into this!
aguspesce: @santisoler and @LL-Geo thanks you for your collaboration on this and finished this code for me. 

Patch:
diff --git a/doc/api/index.rst b/doc/api/index.rst
index cf87e02..cca7df4 100644
--- a/doc/api/index.rst
+++ b/doc/api/index.rst
@@ -125,4 +125,5 @@ Utilities
 .. autosummary::
    :toctree: generated/
 
-    test
+    magnetic_vec_to_angles
+    magnetic_angles_to_vec
diff --git a/harmonica/__init__.py b/harmonica/__init__.py
index 62f6901..97de1f8 100644
--- a/harmonica/__init__.py
+++ b/harmonica/__init__.py
@@ -29,6 +29,7 @@ from ._transformations import (
     reduction_to_pole,
     upward_continuation,
 )
+from ._utils import magnetic_angles_to_vec, magnetic_vec_to_angles
 from ._version import __version__
 
 
diff --git a/harmonica/_utils.py b/harmonica/_utils.py
new file mode 100644
index 0000000..f5107be
--- /dev/null
+++ b/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
diff --git a/pyproject.toml b/pyproject.toml
index 6a5b0e2..edf1519 100644
--- a/pyproject.toml
+++ b/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"
 ]


Test patch:
diff --git a/harmonica/tests/test_utils.py b/harmonica/tests/test_utils.py
new file mode 100644
index 0000000..8c08082
--- /dev/null
+++ b/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)
+        )