From b351df08ca0920b137b8bc515becda42bca70edc Mon Sep 17 00:00:00 2001 From: domfournier Date: Thu, 20 Aug 2026 10:18:48 -0700 Subject: [PATCH 1/7] Add new function to compute azimuth dip from curve segments --- geoapps_utils/utils/locations.py | 38 +++++++++++++++++++++++++++++++- 1 file changed, 37 insertions(+), 1 deletion(-) diff --git a/geoapps_utils/utils/locations.py b/geoapps_utils/utils/locations.py index 573530e..3e1d602 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_spherical + _logger = getLogger(__name__) @@ -219,3 +221,37 @@ def left_limits(n_tiles): limits = left_limits(n_tiles) return limits.tolist() + + +def azimuth_dip_from_segments(curve: Curve) -> tuple[np.ndarray, np.ndarray]: + """ + Compute the local orientation of a Curve object at the vertices, in terms of azimuth and dip. + + Azimuth angles are positive counter-clockwise, referenced from East. + Dip angles are positive downward, referenced from the horizontal plane. + + :param curve: Curve entity to compute azimuth and dip. + + :return: Arrays containing the azimuth and dip angles, in degrees. + """ + delta = curve.vertices[curve.cells[:, 1]] - curve.vertices[curve.cells[:, 0]] + + spherical = cartesian_to_spherical(delta) + + # Convert angles convention + seg_azimuth = (450 - np.rad2deg(spherical[:, 1])) % 360 + seg_dip = np.rad2deg(spherical[:, 2]) - 90 + + # 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_azimuth + dip[nodes, count] = seg_dip + + azimuth = np.nansum(azimuth, axis=1) / (~np.isnan(azimuth)).sum(axis=1) + dip = np.nansum(dip, axis=1) / (~np.isnan(dip)).sum(axis=1) + + return azimuth, dip From 7e615848e9dc7cb2a11b93c1a113ab8ded7f9f84 Mon Sep 17 00:00:00 2001 From: domfournier Date: Thu, 20 Aug 2026 10:20:02 -0700 Subject: [PATCH 2/7] Add unit test --- geoapps_utils/utils/numerical.py | 2 +- tests/locations_test.py | 32 ++++++++++++++++++++++++++++++++ 2 files changed, 33 insertions(+), 1 deletion(-) 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/tests/locations_test.py b/tests/locations_test.py index eb06284..dfc5c7e 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,33 @@ 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 = azimuth_dip_from_segments(curve) + + # For visual check in GA + data = curve.add_data( + { + "azimuth": {"values": azi}, + "dip": {"values": dip}, + } + ) + 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[1:-1], expected[1:-1], decimal=5) From 819fe5c9f7da3843b5659dc302225a224e44c850 Mon Sep 17 00:00:00 2001 From: domfournier Date: Thu, 20 Aug 2026 14:17:19 -0700 Subject: [PATCH 3/7] Add function to convert cartesian to azimuth dip --- geoapps_utils/utils/locations.py | 20 +++++++------------- geoapps_utils/utils/transformations.py | 18 ++++++++++++++++++ tests/locations_test.py | 10 ++++++---- 3 files changed, 31 insertions(+), 17 deletions(-) diff --git a/geoapps_utils/utils/locations.py b/geoapps_utils/utils/locations.py index 3e1d602..d60faff 100644 --- a/geoapps_utils/utils/locations.py +++ b/geoapps_utils/utils/locations.py @@ -22,7 +22,7 @@ from scipy.interpolate import NearestNDInterpolator from scipy.spatial import cKDTree -from geoapps_utils.utils.transformations import cartesian_to_spherical +from geoapps_utils.utils.transformations import cartesian_to_azimuth_dip _logger = getLogger(__name__) @@ -227,20 +227,14 @@ def azimuth_dip_from_segments(curve: Curve) -> tuple[np.ndarray, np.ndarray]: """ Compute the local orientation of a Curve object at the vertices, in terms of azimuth and dip. - Azimuth angles are positive counter-clockwise, referenced from East. - Dip angles are positive downward, referenced from the horizontal plane. - :param curve: Curve entity to compute azimuth and dip. - :return: Arrays containing the azimuth and dip angles, in degrees. + :return: Arrays containing the azimuth anglesm, positive counter-clockwise from East + and dip angles, positive downward from the horizontal plane, in radian. """ delta = curve.vertices[curve.cells[:, 1]] - curve.vertices[curve.cells[:, 0]] - spherical = cartesian_to_spherical(delta) - - # Convert angles convention - seg_azimuth = (450 - np.rad2deg(spherical[:, 1])) % 360 - seg_dip = np.rad2deg(spherical[:, 2]) - 90 + seg_azm_dip = cartesian_to_azimuth_dip(delta) # Average at the node positions azimuth, dip = ( @@ -248,10 +242,10 @@ def azimuth_dip_from_segments(curve: Curve) -> tuple[np.ndarray, np.ndarray]: np.full((curve.n_vertices, 2), np.nan), ) for count, nodes in enumerate(curve.cells.T): - azimuth[nodes, count] = seg_azimuth - dip[nodes, count] = seg_dip + azimuth[nodes, count] = seg_azm_dip[:, 0] + dip[nodes, count] = seg_azm_dip[:, 1] azimuth = np.nansum(azimuth, axis=1) / (~np.isnan(azimuth)).sum(axis=1) dip = np.nansum(dip, axis=1) / (~np.isnan(dip)).sum(axis=1) - return azimuth, dip + return np.c_[azimuth, dip] diff --git a/geoapps_utils/utils/transformations.py b/geoapps_utils/utils/transformations.py index ce45b04..0d93eaa 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 are measured in radians counterclockwise from north. + and dip angles are 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 dfc5c7e..526d039 100644 --- a/tests/locations_test.py +++ b/tests/locations_test.py @@ -209,13 +209,13 @@ def test_azimuth_dip_from_segments(tmp_path): curve = Curve.create(workspace, vertices=xyz) - azi, dip = azimuth_dip_from_segments(curve) + azi_dip = np.rad2deg(azimuth_dip_from_segments(curve)) # For visual check in GA data = curve.add_data( { - "azimuth": {"values": azi}, - "dip": {"values": dip}, + "azimuth": {"values": azi_dip[:, 0]}, + "dip": {"values": azi_dip[:, 1]}, } ) PropertyGroup( @@ -225,4 +225,6 @@ def test_azimuth_dip_from_segments(tmp_path): # 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[1:-1], expected[1:-1], decimal=5) + np.testing.assert_array_almost_equal( + azi_dip[1:-1, 0], expected[1:-1], decimal=5 + ) From a542007862c5d259e17e07f5c3d2bb55e3981fed Mon Sep 17 00:00:00 2001 From: domfournier Date: Thu, 20 Aug 2026 14:41:42 -0700 Subject: [PATCH 4/7] Add unit test --- tests/transformations_test.py | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) 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]) From 50bdf7d14f8c33c4aee3f6f07652a45066703b5f Mon Sep 17 00:00:00 2001 From: domfournier Date: Wed, 26 Aug 2026 13:39:58 -0700 Subject: [PATCH 5/7] Add option to get the flipped direction on segments --- geoapps_utils/utils/locations.py | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/geoapps_utils/utils/locations.py b/geoapps_utils/utils/locations.py index d60faff..ff77337 100644 --- a/geoapps_utils/utils/locations.py +++ b/geoapps_utils/utils/locations.py @@ -223,17 +223,23 @@ def left_limits(n_tiles): return limits.tolist() -def azimuth_dip_from_segments(curve: Curve) -> tuple[np.ndarray, np.ndarray]: +def azimuth_dip_from_segments( + curve: Curve, reverse=False +) -> tuple[np.ndarray, 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 anglesm, positive counter-clockwise from East + :return: Arrays containing the azimuth angles, positive counter-clockwise from East 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 From 92c85f97226605171890b43a5fb6dc1a3bfa03ae Mon Sep 17 00:00:00 2001 From: domfournier Date: Mon, 31 Aug 2026 11:52:40 -0700 Subject: [PATCH 6/7] Handle transition between -pi to pi --- geoapps_utils/utils/locations.py | 9 +++++++-- tests/locations_test.py | 8 ++++++++ 2 files changed, 15 insertions(+), 2 deletions(-) diff --git a/geoapps_utils/utils/locations.py b/geoapps_utils/utils/locations.py index ff77337..71f9431 100644 --- a/geoapps_utils/utils/locations.py +++ b/geoapps_utils/utils/locations.py @@ -251,7 +251,12 @@ def azimuth_dip_from_segments( azimuth[nodes, count] = seg_azm_dip[:, 0] dip[nodes, count] = seg_azm_dip[:, 1] - azimuth = np.nansum(azimuth, axis=1) / (~np.isnan(azimuth)).sum(axis=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_[azimuth, dip] + return np.c_[direction, dip] diff --git a/tests/locations_test.py b/tests/locations_test.py index 526d039..2518fd5 100644 --- a/tests/locations_test.py +++ b/tests/locations_test.py @@ -228,3 +228,11 @@ def test_azimuth_dip_from_segments(tmp_path): 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]) From 33b254c44ab5daf494bfc36368dadfb16eff5952 Mon Sep 17 00:00:00 2001 From: domfournier Date: Mon, 31 Aug 2026 12:03:22 -0700 Subject: [PATCH 7/7] copilot review --- geoapps_utils/utils/locations.py | 6 ++---- geoapps_utils/utils/transformations.py | 4 ++-- 2 files changed, 4 insertions(+), 6 deletions(-) diff --git a/geoapps_utils/utils/locations.py b/geoapps_utils/utils/locations.py index 71f9431..8f9b6f6 100644 --- a/geoapps_utils/utils/locations.py +++ b/geoapps_utils/utils/locations.py @@ -223,16 +223,14 @@ def left_limits(n_tiles): return limits.tolist() -def azimuth_dip_from_segments( - curve: Curve, reverse=False -) -> tuple[np.ndarray, np.ndarray]: +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 counter-clockwise from East + :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]] diff --git a/geoapps_utils/utils/transformations.py b/geoapps_utils/utils/transformations.py index 0d93eaa..c4aa3fb 100644 --- a/geoapps_utils/utils/transformations.py +++ b/geoapps_utils/utils/transformations.py @@ -136,8 +136,8 @@ def cartesian_to_azimuth_dip(vectors: np.ndarray) -> np.ndarray: :param vectors: Array of shape (n, 3) representing x, y, z coordinates of a vector. - :return: Arrays of azimuth angles are measured in radians counterclockwise from north. - and dip angles are measured in radians positive downward from the horizontal plane. + :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)