diff --git a/geoapps_utils/utils/locations.py b/geoapps_utils/utils/locations.py index 573530e..8f9b6f6 100644 --- a/geoapps_utils/utils/locations.py +++ b/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] diff --git a/geoapps_utils/utils/numerical.py b/geoapps_utils/utils/numerical.py index 2c351fb..50b6732 100644 --- a/geoapps_utils/utils/numerical.py +++ b/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] diff --git a/geoapps_utils/utils/transformations.py b/geoapps_utils/utils/transformations.py index ce45b04..c4aa3fb 100644 --- a/geoapps_utils/utils/transformations.py +++ b/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. diff --git a/tests/locations_test.py b/tests/locations_test.py index eb06284..2518fd5 100644 --- a/tests/locations_test.py +++ b/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]) diff --git a/tests/transformations_test.py b/tests/transformations_test.py index 3621330..cc1163b 100644 --- a/tests/transformations_test.py +++ b/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])