GEOPY-3046: Implement borehole EM inversion

geoapps_utils/utils/locations.py

@@ -16,12 +16,14 @@
 import numpy as np
 from geoh5py import Workspace
 from geoh5py.data import Data
-from geoh5py.objects import CellObject, Grid2D, Points
+from geoh5py.objects import CellObject, Curve, Grid2D, Points
 from geoh5py.objects.grid_object import GridObject
 from matplotlib.tri import LinearTriInterpolator, Triangulation
 from scipy.interpolate import NearestNDInterpolator
 from scipy.spatial import cKDTree
 
+from geoapps_utils.utils.transformations import cartesian_to_azimuth_dip
+
 
 _logger = getLogger(__name__)
 
@@ -219,3 +221,40 @@ def left_limits(n_tiles):
         limits = left_limits(n_tiles)
 
     return limits.tolist()
+
+
+def azimuth_dip_from_segments(curve: Curve, reverse=False) -> np.ndarray:
+    """
+    Compute the local orientation of a Curve object at the vertices, in terms of azimuth and dip.
+
+    :param curve: Curve entity to compute azimuth and dip.
+    :param reverse: Reverse the direction of the segments.
+
+    :return: Arrays containing the azimuth angles, positive clockwise from north
+        and dip angles, positive downward from the horizontal plane, in radian.
+    """
+    delta = curve.vertices[curve.cells[:, 1]] - curve.vertices[curve.cells[:, 0]]
+
+    if reverse:
+        delta = -delta
+
+    seg_azm_dip = cartesian_to_azimuth_dip(delta)
+
+    # Average at the node positions
+    azimuth, dip = (
+        np.full((curve.n_vertices, 2), np.nan),
+        np.full((curve.n_vertices, 2), np.nan),
+    )
+    for count, nodes in enumerate(curve.cells.T):
+        azimuth[nodes, count] = seg_azm_dip[:, 0]
+        dip[nodes, count] = seg_azm_dip[:, 1]
+
+    # Deal with (-pi, pi) or (0, 2pi) transition
+    d_azm = np.diff(azimuth, axis=1) % (2 * np.pi)
+    direction = azimuth[:, 0] + d_azm.flatten() / 2.0
+    end_lines = np.where(np.isnan(azimuth).sum(axis=1))
+    direction[end_lines] = np.nansum(azimuth, axis=1)[end_lines]
+
+    dip = np.nansum(dip, axis=1) / (~np.isnan(dip)).sum(axis=1)
+
+    return np.c_[direction, dip]

geoapps_utils/utils/numerical.py

@@ -130,7 +130,7 @@ def weighted_average(
     avg_values = []
     for value in values:
         sub = ~np.isnan(value)
-        rad, ind = cKDTree(xyz_in[sub]).query(xyz_out, n)
+        rad, ind = cKDTree(xyz_in[sub]).query(xyz_out, n, workers=-1)
 
         if n == 1:
             ind = ind[:, np.newaxis]

geoapps_utils/utils/transformations.py

@@ -130,6 +130,24 @@ def cartesian_to_spherical(points: np.ndarray) -> np.ndarray:
     return np.column_stack((magnitude, azimuth, inclination))
 
 
+def cartesian_to_azimuth_dip(vectors: np.ndarray) -> np.ndarray:
+    """
+    Convert vectors in cartesian coordinates to azimuth and dip in radian.
+
+    :param vectors: Array of shape (n, 3) representing x, y, z coordinates of a vector.
+
+    :return: Arrays of azimuth angles measured in radians clockwise from north
+        and dip angles measured in radians positive downward from the horizontal plane.
+    """
+    spherical = cartesian_to_spherical(vectors)
+
+    # Convert angles convention
+    azimuth = ccw_east_to_cw_north(spherical[:, 1])
+    dip = spherical[:, 2] - np.pi / 2
+
+    return np.c_[azimuth, dip]
+
+
 def ccw_east_to_cw_north(azimuth: np.ndarray) -> np.ndarray:
     """
     Convert counterclockwise azimuth from east to clockwise from north.

tests/locations_test.py

@@ -16,10 +16,12 @@
 import numpy as np
 import pytest
 from geoh5py import Workspace
+from geoh5py.groups import PropertyGroup
 from geoh5py.objects import Curve, Grid2D, Points
 from scipy.spatial import Delaunay
 
 from geoapps_utils.utils.locations import (
+    azimuth_dip_from_segments,
     gaussian,
     get_locations,
     get_overlapping_limits,
@@ -194,3 +196,43 @@ def test_overlapping_limits():
     limits = get_overlapping_limits(width * 4, width, overlap=overlap)
 
     assert limits[0][1] == limits[1][0]
+
+
+def test_azimuth_dip_from_segments(tmp_path):
+    with Workspace.create(tmp_path / f"{__name__}.geoh5") as workspace:
+        n = 30
+        theta = np.linspace(0, np.pi * 2, n)
+
+        x, y = np.cos(theta), np.sin(theta)
+        z = np.sin(theta)
+        xyz = np.c_[x.ravel(), y.ravel(), z.ravel()]  # Form a 2D array
+
+        curve = Curve.create(workspace, vertices=xyz)
+
+        azi_dip = np.rad2deg(azimuth_dip_from_segments(curve))
+
+        # For visual check in GA
+        data = curve.add_data(
+            {
+                "azimuth": {"values": azi_dip[:, 0]},
+                "dip": {"values": azi_dip[:, 1]},
+            }
+        )
+        PropertyGroup(
+            parent=curve, property_group_type="Dip direction & dip", properties=data
+        )
+
+        # Check that angles are within expected ranges expect first and last point
+        expected = (450 - np.rad2deg(theta) - 90) % 360
+
+        np.testing.assert_array_almost_equal(
+            azi_dip[1:-1, 0], expected[1:-1], decimal=5
+        )
+
+        # Test reverse flag
+        azi_dip_reversed = np.rad2deg(azimuth_dip_from_segments(curve, reverse=True))
+
+        np.testing.assert_array_almost_equal(
+            (azi_dip_reversed[:, 0] + 180) % 360, azi_dip[:, 0]
+        )
+        np.testing.assert_array_almost_equal(-azi_dip_reversed[:, 1], azi_dip[:, 1])

tests/transformations_test.py

@@ -19,6 +19,7 @@
 from geoapps_utils.utils.transformations import (
     azimuth_to_unit_vector,
     cartesian_normal_to_direction_and_dip,
+    cartesian_to_azimuth_dip,
     cartesian_to_polar,
     cartesian_to_spherical,
     compute_normals,
@@ -227,3 +228,22 @@ def test_cartesian_to_polar():
     np.testing.assert_allclose(
         polar[:, 1], 90 - np.rad2deg(azm) + 180
     )  # All other distances positive
+
+
+@pytest.mark.parametrize(
+    "vector,expected",
+    [
+        ((1, 0, 0), [90, 0]),
+        ((0, 1, 0), [0, 0]),
+        ((0, 0, 1), [90, -90]),
+        ((-1, 0, 0), [270, 0]),
+        ((-1, -1, 0), [225, 0]),
+        ((0.5, 0.5, -0.7071), [45, 45]),
+        ((-0.5, 0.5, -0.7071), [315, 45]),
+    ],
+)
+def test_cartesian_to_azimuth_dip(vector, expected):
+    vector = np.asarray(vector).reshape((1, 3))
+    azm_dip = cartesian_to_azimuth_dip(vector)
+    angles = np.rad2deg(azm_dip).flatten()
+    assert np.allclose(angles, np.r_[expected])