Fix vector to angles for small intensities (#565)

When the intensity is very small (<1e-10), the vector to angles
calculation fails because of the masking that was used. It rounds these
numbers to zero and returns always inc=90, dec=0. To handle this, use
the arctangent instead of arcsine in the equations. This was, we can use
arctan2 which is able to handle 0 in the denominator. Add tests to make
sure this doesn't happen again. These small values appear in magnetic
microscopy and we ran into this issue in Magali.

Author: leouieda

harmonica/_utils.py

@@ -68,7 +68,6 @@ def magnetic_vec_to_angles(magnetic_e, magnetic_n, magnetic_u, degrees=True):
         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
@@ -119,7 +118,7 @@ def magnetic_vec_to_angles(magnetic_e, magnetic_n, magnetic_u, degrees=True):
 
     .. math::
 
-        D = \arcsin \frac{B_e}{\sqrt{B_e^2 + B_n^2}}
+        D = \arctan \frac{B_e}{B_n}
 
     Examples
     --------
@@ -130,23 +129,14 @@ def magnetic_vec_to_angles(magnetic_e, magnetic_n, magnetic_u, degrees=True):
     # 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)
+    horizontal_component = np.sqrt(magnetic_e**2 + magnetic_n**2)
+    # Calculate the inclination and declination
+    inclination = np.arctan2(-magnetic_u, horizontal_component)
+    declination = np.arctan2(magnetic_e, magnetic_n)
     # 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

@@ -17,9 +17,22 @@
     [0, 0, -1],  # Over -z axis
     [1, 0, 0],  # Over east (y) axis
     [0, 1, 0],  # Over north (x) axis
+    [0, 0, -1e-13],  # Over -z axis
+    [1e-13, 0, 0],  # Over east (y) axis
+    [0, 1e-13, 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]]
+ANGLES = [
+    [1, 45, 45],
+    [1, -45, 45.0],
+    [1, 45, -45],
+    [1, 90, 0],
+    [1, 0, 90],
+    [1, 0, 0],
+    [1e-13, 90, 0],
+    [1e-13, 0, 90],
+    [1e-13, 0, 0],
+]
 
 
 @pytest.mark.parametrize("angles, vector", [(a, v) for a, v in zip(ANGLES, VECTORS)])