GEOPY-2108: Support 3D vector property group for orientation of rotated gradients

geoapps_utils/utils/transformations.py

@@ -13,13 +13,14 @@
 
 def rotate_xyz(xyz: np.ndarray, center: list, theta: float, phi: float = 0.0):
     """
-    Perform a counterclockwise rotation of scatter points around the x-axis,
-        then z-axis, about a center point.
+    Rotate points counterclockwise around the x then z axes, about a center point.
 
     :param xyz: shape(*, 2) or shape(*, 3) Input coordinates.
     :param center: len(2) or len(3) Coordinates for the center of rotation.
-    :param theta: Angle of rotation around z-axis in degrees.
-    :param phi: Angle of rotation around x-axis in degrees.
+    :param theta: Angle of rotation in counterclockwise degree about the z-axis
+        as viewed from above.
+    :param phi: Angle of rotation in couterclockwise degrees around x-axis
+        as viewed from the east.
     """
     return2d = False
     locs = xyz.copy()
@@ -55,3 +56,80 @@
         # Return 2-dimensional data if the input xyz was 2-dimensional.
         return xyz_out[:, :2]
     return xyz_out
+
+
+def ccw_east_to_cw_north(azimuth: float) -> float:
+    """
+    Convert counterclockwise azimuth from east to clockwise from north
+
+    :param azimuth: Azimuth angle (in xy plane) measured in counterclockwise radians
+        from east.
+
+    :returns: Azimuth angle (in xy plane) measured in clockwise degrees from north.
+    """
+    return (((5 * np.pi) / 2) - azimuth) % (2 * np.pi)
+
+
+def cartesian_to_spherical(points: np.ndarray) -> np.ndarray:
+    """
+    Convert cartesian to spherical coordinates.
+
+    :param points: Array of shape (n, 3) representing x, y, z coordinates of a point
+        in 3D space.
+
+    :returns: Array of shape (n, 2) representing the azimuth and inclination angles
+        in spherical coordinates. The azimuth angle is measured in radians clockwise
+        from north in the range of 0 to 2pi as viewed from above, and inclination
+        angle is measured in positive radians above the horizon and negative below.
+    """
+    inclination = np.arcsin(points[:, 2] / np.linalg.norm(points, axis=1))
+    azimuth = np.sign(points[:, 1]) * (
+        np.arccos(points[:, 0] / np.linalg.norm(points[:, :2], axis=1))
+    )
+    azimuth = ccw_east_to_cw_north(azimuth)
+    return np.column_stack((azimuth, inclination))
+
+
+def spherical_to_direction_and_dip(angles: np.ndarray) -> np.ndarray:
+    """
+    Convert normals in spherical coordinates to dip and direction of the tangent plane.
+
+    Confines the solution to the eastern hemisphere and applies a dip
+    correction for normals originally oriented in the west.
+
+    :param angles: Array of shape (n, 2) representing the azimuth and inclination angles
+        in spherical coordinates. The azimuth angle is measured in radians clockwise
+        from north in the range of 0 to 2pi as viewed from above, and inclination
+        angle is measured in positive radians above the horizon and negative below.
+
+    :returns: Array of shape (n, 2) representing azimuth from 0 to pi radians
+        clockwise from north as viewed from above and dip from -pi to pi in positive
+        radians below the horizon and negative above.
+
+    """
+    azimuth = angles[:, 0]
+    inclination = angles[:, 1]
+    inclination = np.sign(inclination) * ((np.pi / 2) - np.abs(inclination))
+    greater_than_pi = azimuth > np.pi
+    inclination[greater_than_pi] = -1 * (inclination[greater_than_pi])
+    azimuth[greater_than_pi] = azimuth[greater_than_pi] - np.pi
+
+    return np.column_stack((azimuth, inclination))
+
+
+def normal_to_direction_and_dip(points: np.ndarray) -> np.ndarray:
+    """
+    Convert 3D normal vectors to dip and direction within the eastern hemisphere.
+
+    :param points: Array of shape (n, 3) representing x, y, z coordinates of a point
+        in 3D space.
+
+    :returns: Array of shape (n, 2) representing azimuth from 0 to pi radians
+        clockwise from north as viewed from above and dip from -pi to pi in positive
+        radians below the horizon and negative above.
+    """
+
+    spherical_coords = cartesian_to_spherical(points)
+    direction_and_dip = spherical_to_direction_and_dip(spherical_coords)
+
+    return direction_and_dip

tests/transformations_test.py

@@ -12,7 +12,11 @@
 
 import numpy as np
 
-from geoapps_utils.utils.transformations import rotate_xyz
+from geoapps_utils.utils.transformations import (
+    cartesian_to_spherical,
+    rotate_xyz,
+    spherical_to_direction_and_dip,
+)
 
 
 def test_positive_rotation_xyz():
@@ -33,3 +37,107 @@ def test_negative_rotation_xyz():
 
 def test_2d_input():
     assert (rotate_xyz(np.c_[1, 0], [0, 0], 45)).shape[1] == 2, "Error on 2D input."
+
+
+def test_cartesian_to_spherical_first_quadrant_upwards():
+    pole = rotate_xyz(xyz=np.c_[0, 0, 1], center=[0, 0, 0], theta=-60, phi=-30)
+    angles = np.rad2deg(cartesian_to_spherical(pole))
+    assert np.allclose(angles, [[60, 60]])
+
+
+def test_cartesian_to_spherical_second_quadrant_upwards():
+    pole = rotate_xyz(xyz=np.c_[0, 0, 1], center=[0, 0, 0], theta=-150, phi=-30)
+    angles = np.rad2deg(cartesian_to_spherical(pole))
+    assert np.allclose(angles, [[150, 60]])
+
+
+def test_cartesian_to_spherical_third_quadrant_upwards():
+    pole = rotate_xyz(xyz=np.c_[0, 0, 1], center=[0, 0, 0], theta=-240, phi=-30)
+    angles = np.rad2deg(cartesian_to_spherical(pole))
+    assert np.allclose(angles, [[240, 60]])
+
+
+def test_cartesian_to_spherical_fourth_quadrant_upwards():
+    pole = rotate_xyz(xyz=np.c_[0, 0, 1], center=[0, 0, 0], theta=-330, phi=-30)
+    angles = np.rad2deg(cartesian_to_spherical(pole))
+    assert np.allclose(angles, [[330, 60]])
+
+
+def test_cartesian_to_spherical_first_quadrant_downwards():
+    pole = rotate_xyz(xyz=np.c_[0, 0, -1], center=[0, 0, 0], theta=-60, phi=30)
+    angles = np.rad2deg(cartesian_to_spherical(pole))
+    assert np.allclose(angles, [[60, -60]])
+
+
+def test_cartesian_to_spherical_second_quadrant_downwards():
+    pole = rotate_xyz(xyz=np.c_[0, 0, -1], center=[0, 0, 0], theta=-150, phi=30)
+    angles = np.rad2deg(cartesian_to_spherical(pole))
+    assert np.allclose(angles, [[150, -60]])
+
+
+def test_cartesian_to_spherical_third_quadrant_downwards():
+    pole = rotate_xyz(xyz=np.c_[0, 0, -1], center=[0, 0, 0], theta=-240, phi=30)
+    angles = np.rad2deg(cartesian_to_spherical(pole))
+    assert np.allclose(angles, [[240, -60]])
+
+
+def test_cartesian_to_spherical_fourth_quadrant_downwards():
+    pole = rotate_xyz(xyz=np.c_[0, 0, -1], center=[0, 0, 0], theta=-330, phi=30)
+    angles = np.rad2deg(cartesian_to_spherical(pole))
+    assert np.allclose(angles, [[330, -60]])
+
+
+def test_spherical_to_direction_and_dip_first_quadrant_upwards():
+    pole = rotate_xyz(xyz=np.c_[0, 0, 1], center=[0, 0, 0], theta=-60, phi=-30)
+    spherical_coords = cartesian_to_spherical(pole)
+    angles = np.rad2deg(spherical_to_direction_and_dip(spherical_coords))
+    assert np.allclose(angles, [[60, 30]])
+
+
+def test_spherical_to_direction_and_dip_second_quadrant_upwards():
+    pole = rotate_xyz(xyz=np.c_[0, 0, 1], center=[0, 0, 0], theta=-150, phi=-30)
+    spherical_coords = cartesian_to_spherical(pole)
+    angles = np.rad2deg(spherical_to_direction_and_dip(spherical_coords))
+    assert np.allclose(angles, [[150, 30]])
+
+
+def test_spherical_to_direction_and_dip_third_quadrant_upwards():
+    pole = rotate_xyz(xyz=np.c_[0, 0, 1], center=[0, 0, 0], theta=-240, phi=-30)
+    spherical_coords = cartesian_to_spherical(pole)
+    angles = np.rad2deg(spherical_to_direction_and_dip(spherical_coords))
+    assert np.allclose(angles, [[60, -30]])
+
+
+def test_spherical_to_direction_and_dip_fourth_quadrant_upwards():
+    pole = rotate_xyz(xyz=np.c_[0, 0, 1], center=[0, 0, 0], theta=-330, phi=-30)
+    spherical_coords = cartesian_to_spherical(pole)
+    angles = np.rad2deg(spherical_to_direction_and_dip(spherical_coords))
+    assert np.allclose(angles, [[150, -30]])
+
+
+def test_spherical_to_direction_and_dip_first_quadrant_downwards():
+    pole = rotate_xyz(xyz=np.c_[0, 0, -1], center=[0, 0, 0], theta=-60, phi=30)
+    spherical_coords = cartesian_to_spherical(pole)
+    angles = np.rad2deg(spherical_to_direction_and_dip(spherical_coords))
+    assert np.allclose(angles, [[60, -30]])
+
+
+def test_spherical_to_direction_and_dip_second_quadrant_downwards():
+    pole = rotate_xyz(xyz=np.c_[0, 0, -1], center=[0, 0, 0], theta=-150, phi=30)
+    spherical_coords = cartesian_to_spherical(pole)
+    angles = np.rad2deg(spherical_to_direction_and_dip(spherical_coords))
+    assert np.allclose(angles, [[150, -30]])
+
+
+def test_spherical_to_direction_and_dip_third_quadrant_downwards():
+    pole = rotate_xyz(xyz=np.c_[0, 0, -1], center=[0, 0, 0], theta=-240, phi=30)
+    spherical_coords = cartesian_to_spherical(pole)
+    angles = np.rad2deg(spherical_to_direction_and_dip(spherical_coords))
+    assert np.allclose(angles, [[60, 30]])
+
+
+def test_spherical_to_direction_and_dip_fourth_quadrant_downwards():
+    pole = rotate_xyz(xyz=np.c_[0, 0, -1], center=[0, 0, 0], theta=-330, phi=30)
+    spherical_coords = cartesian_to_spherical(pole)
+    angles = np.rad2deg(spherical_to_direction_and_dip(spherical_coords))
+    assert np.allclose(angles, [[150, 30]])