diff --git a/LICENSE b/LICENSE index 9c6ee51f..d60e86b7 100644 --- a/LICENSE +++ b/LICENSE @@ -1,6 +1,6 @@ MIT License -Copyright (c) 2023-2025 Mira Geoscience +Copyright (c) 2022-2026 Mira Geoscience Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal diff --git a/README-dev.rst b/README-dev.rst index e0d49c2c..94386238 100644 --- a/README-dev.rst +++ b/README-dev.rst @@ -335,4 +335,4 @@ Here is a suggestion for some plugins you can install in PyCharm. Copyright ^^^^^^^^^ -Copyright (c) 2023-2025 Mira Geoscience Ltd. +Copyright (c) 2022-2026 Mira Geoscience Ltd. diff --git a/README.rst b/README.rst index 0571164c..e6772466 100644 --- a/README.rst +++ b/README.rst @@ -149,7 +149,7 @@ License ^^^^^^^ MIT License -Copyright (c) 2023-2024 Mira Geoscience +Copyright (c) 2022-2026 Mira Geoscience Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal diff --git a/geoapps_utils-assets/__init__.py b/geoapps_utils-assets/__init__.py index 348e9ed6..02fcb1b2 100644 --- a/geoapps_utils-assets/__init__.py +++ b/geoapps_utils-assets/__init__.py @@ -1,5 +1,5 @@ # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2022-2025 Mira Geoscience Ltd. ' +# Copyright (c) 2022-2026 Mira Geoscience Ltd. ' # ' # This file is part of geoapps-utils package. ' # ' diff --git a/geoapps_utils/__init__.py b/geoapps_utils/__init__.py index 2ce78e42..d55d9265 100644 --- a/geoapps_utils/__init__.py +++ b/geoapps_utils/__init__.py @@ -1,5 +1,5 @@ # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2022-2025 Mira Geoscience Ltd. ' +# Copyright (c) 2022-2026 Mira Geoscience Ltd. ' # ' # This file is part of geoapps-utils package. ' # ' diff --git a/geoapps_utils/base.py b/geoapps_utils/base.py index bc6cf70a..f6c2c631 100644 --- a/geoapps_utils/base.py +++ b/geoapps_utils/base.py @@ -1,5 +1,5 @@ # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2023-2025 Mira Geoscience Ltd. ' +# Copyright (c) 2022-2026 Mira Geoscience Ltd. ' # ' # This file is part of geoapps-utils package. ' # ' diff --git a/geoapps_utils/driver/__init__.py b/geoapps_utils/driver/__init__.py index 76e3205e..02fcb1b2 100644 --- a/geoapps_utils/driver/__init__.py +++ b/geoapps_utils/driver/__init__.py @@ -1,5 +1,5 @@ # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2025 Mira Geoscience Ltd. ' +# Copyright (c) 2022-2026 Mira Geoscience Ltd. ' # ' # This file is part of geoapps-utils package. ' # ' diff --git a/geoapps_utils/driver/data.py b/geoapps_utils/driver/data.py index f4c83bb9..37ec4cf5 100644 --- a/geoapps_utils/driver/data.py +++ b/geoapps_utils/driver/data.py @@ -1,5 +1,5 @@ # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2023-2025 Mira Geoscience Ltd. ' +# Copyright (c) 2022-2026 Mira Geoscience Ltd. ' # ' # This file is part of geoapps-utils package. ' # ' diff --git a/geoapps_utils/driver/driver.py b/geoapps_utils/driver/driver.py index fbac38c2..3de55748 100644 --- a/geoapps_utils/driver/driver.py +++ b/geoapps_utils/driver/driver.py @@ -1,5 +1,5 @@ # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2025 Mira Geoscience Ltd. ' +# Copyright (c) 2022-2026 Mira Geoscience Ltd. ' # ' # This file is part of geoapps-utils package. ' # ' diff --git a/geoapps_utils/driver/params.py b/geoapps_utils/driver/params.py index ebdcf3aa..fd935e5b 100644 --- a/geoapps_utils/driver/params.py +++ b/geoapps_utils/driver/params.py @@ -1,5 +1,5 @@ # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2023-2025 Mira Geoscience Ltd. ' +# Copyright (c) 2022-2026 Mira Geoscience Ltd. ' # ' # This file is part of geoapps-utils package. ' # ' diff --git a/geoapps_utils/modelling/__init__.py b/geoapps_utils/modelling/__init__.py index 76e3205e..02fcb1b2 100644 --- a/geoapps_utils/modelling/__init__.py +++ b/geoapps_utils/modelling/__init__.py @@ -1,5 +1,5 @@ # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2025 Mira Geoscience Ltd. ' +# Copyright (c) 2022-2026 Mira Geoscience Ltd. ' # ' # This file is part of geoapps-utils package. ' # ' diff --git a/geoapps_utils/modelling/plates.py b/geoapps_utils/modelling/plates.py index 6c678c99..bff06e04 100644 --- a/geoapps_utils/modelling/plates.py +++ b/geoapps_utils/modelling/plates.py @@ -1,5 +1,5 @@ # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2025 Mira Geoscience Ltd. ' +# Copyright (c) 2022-2026 Mira Geoscience Ltd. ' # ' # This file is part of geoapps-utils package. ' # ' diff --git a/geoapps_utils/param_sweeps/__init__.py b/geoapps_utils/param_sweeps/__init__.py index 5197aa21..02fcb1b2 100644 --- a/geoapps_utils/param_sweeps/__init__.py +++ b/geoapps_utils/param_sweeps/__init__.py @@ -1,5 +1,5 @@ # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2023-2025 Mira Geoscience Ltd. ' +# Copyright (c) 2022-2026 Mira Geoscience Ltd. ' # ' # This file is part of geoapps-utils package. ' # ' diff --git a/geoapps_utils/param_sweeps/driver.py b/geoapps_utils/param_sweeps/driver.py index 005cec1b..e6f157ec 100644 --- a/geoapps_utils/param_sweeps/driver.py +++ b/geoapps_utils/param_sweeps/driver.py @@ -1,5 +1,5 @@ # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2023-2025 Mira Geoscience Ltd. ' +# Copyright (c) 2022-2026 Mira Geoscience Ltd. ' # ' # This file is part of geoapps-utils package. ' # ' diff --git a/geoapps_utils/param_sweeps/generate.py b/geoapps_utils/param_sweeps/generate.py index b7a0f031..6b2a9e35 100644 --- a/geoapps_utils/param_sweeps/generate.py +++ b/geoapps_utils/param_sweeps/generate.py @@ -1,5 +1,5 @@ # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2023-2025 Mira Geoscience Ltd. ' +# Copyright (c) 2022-2026 Mira Geoscience Ltd. ' # ' # This file is part of geoapps-utils package. ' # ' diff --git a/geoapps_utils/run.py b/geoapps_utils/run.py index 8e9eea92..1c3ad3ee 100644 --- a/geoapps_utils/run.py +++ b/geoapps_utils/run.py @@ -1,5 +1,5 @@ # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2023-2025 Mira Geoscience Ltd. ' +# Copyright (c) 2022-2026 Mira Geoscience Ltd. ' # ' # This file is part of geoapps-utils package. ' # ' diff --git a/geoapps_utils/utils/__init__.py b/geoapps_utils/utils/__init__.py index a23cf840..edb8f4c5 100644 --- a/geoapps_utils/utils/__init__.py +++ b/geoapps_utils/utils/__init__.py @@ -1,5 +1,5 @@ # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2024-2025 Mira Geoscience Ltd. ' +# Copyright (c) 2024-2026 Mira Geoscience Ltd. ' # ' # This file is part of geoapps-utils package. ' # ' diff --git a/geoapps_utils/utils/conversions.py b/geoapps_utils/utils/conversions.py index 6225eb7e..2f7c9629 100644 --- a/geoapps_utils/utils/conversions.py +++ b/geoapps_utils/utils/conversions.py @@ -1,5 +1,5 @@ # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2023-2025 Mira Geoscience Ltd. ' +# Copyright (c) 2022-2026 Mira Geoscience Ltd. ' # ' # This file is part of geoapps-utils package. ' # ' diff --git a/geoapps_utils/utils/formatters.py b/geoapps_utils/utils/formatters.py index 18912675..f4346cfe 100644 --- a/geoapps_utils/utils/formatters.py +++ b/geoapps_utils/utils/formatters.py @@ -1,5 +1,5 @@ # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2023-2025 Mira Geoscience Ltd. ' +# Copyright (c) 2022-2026 Mira Geoscience Ltd. ' # ' # This file is part of geoapps-utils package. ' # ' diff --git a/geoapps_utils/utils/importing.py b/geoapps_utils/utils/importing.py index 6a0dfe91..7021d20c 100644 --- a/geoapps_utils/utils/importing.py +++ b/geoapps_utils/utils/importing.py @@ -1,5 +1,5 @@ # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2023-2025 Mira Geoscience Ltd. ' +# Copyright (c) 2022-2026 Mira Geoscience Ltd. ' # ' # This file is part of geoapps-utils package. ' # ' diff --git a/geoapps_utils/utils/iterables.py b/geoapps_utils/utils/iterables.py index cb33a242..5d756571 100644 --- a/geoapps_utils/utils/iterables.py +++ b/geoapps_utils/utils/iterables.py @@ -1,5 +1,5 @@ # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2023-2025 Mira Geoscience Ltd. ' +# Copyright (c) 2022-2026 Mira Geoscience Ltd. ' # ' # This file is part of geoapps-utils package. ' # ' diff --git a/geoapps_utils/utils/locations.py b/geoapps_utils/utils/locations.py index 45ec72ce..573530ed 100644 --- a/geoapps_utils/utils/locations.py +++ b/geoapps_utils/utils/locations.py @@ -1,5 +1,5 @@ # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2023-2025 Mira Geoscience Ltd. ' +# Copyright (c) 2022-2026 Mira Geoscience Ltd. ' # ' # This file is part of geoapps-utils package. ' # ' diff --git a/geoapps_utils/utils/logger.py b/geoapps_utils/utils/logger.py index 72914742..48405e2e 100644 --- a/geoapps_utils/utils/logger.py +++ b/geoapps_utils/utils/logger.py @@ -1,5 +1,5 @@ # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2025 Mira Geoscience Ltd. ' +# Copyright (c) 2022-2026 Mira Geoscience Ltd. ' # ' # This file is part of geoapps-utils package. ' # ' diff --git a/geoapps_utils/utils/numerical.py b/geoapps_utils/utils/numerical.py index 1b2bd04a..2c351fb4 100644 --- a/geoapps_utils/utils/numerical.py +++ b/geoapps_utils/utils/numerical.py @@ -1,5 +1,5 @@ # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2023-2025 Mira Geoscience Ltd. ' +# Copyright (c) 2022-2026 Mira Geoscience Ltd. ' # ' # This file is part of geoapps-utils package. ' # ' @@ -11,6 +11,7 @@ from __future__ import annotations import numpy as np +from scipy.sparse import csr_matrix, diags from scipy.spatial import cKDTree @@ -202,3 +203,33 @@ def fit_circle(x_val: np.ndarray, y_val) -> tuple[float, float, float]: radius = (coef[0] ** 2.0 + coef[1] ** 2.0 + coef[2]) ** 0.5 return radius, coef[0], coef[1] + + +def inverse_weighted_operator( + values: np.ndarray, + col_indices: np.ndarray, + shape: tuple, + power: float, + threshold: float, +) -> csr_matrix: + """ + Create an inverse distance weighted sparse matrix. + + :param values: Distance values. + :param col_indices: Column indices for the sparse matrix. + :param shape: Shape of the sparse matrix. + :param power: Power for the inverse distance weighting. + :param threshold: Threshold to avoid singularities. + + :return: Inverse distance weighted sparse matrix. + """ + weights = (values**power + threshold) ** -1 + n_vals_row = weights.shape[0] // shape[0] + row_ids = np.repeat(np.arange(shape[0]), n_vals_row) + inv_dist_op = csr_matrix( + (weights, (row_ids, col_indices)), + shape=shape, + ) + # Normalize the rows + row_sum = np.asarray(inv_dist_op.sum(axis=1)).flatten() ** -1.0 + return diags(row_sum) @ inv_dist_op diff --git a/geoapps_utils/utils/plotting.py b/geoapps_utils/utils/plotting.py index 5957f8aa..5d601f8b 100644 --- a/geoapps_utils/utils/plotting.py +++ b/geoapps_utils/utils/plotting.py @@ -1,5 +1,5 @@ # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2023-2025 Mira Geoscience Ltd. ' +# Copyright (c) 2022-2026 Mira Geoscience Ltd. ' # ' # This file is part of geoapps-utils package. ' # ' diff --git a/geoapps_utils/utils/reference.py b/geoapps_utils/utils/reference.py index eef30edd..07d5067e 100644 --- a/geoapps_utils/utils/reference.py +++ b/geoapps_utils/utils/reference.py @@ -1,5 +1,5 @@ # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2024-2025 Mira Geoscience Ltd. ' +# Copyright (c) 2024-2026 Mira Geoscience Ltd. ' # ' # This file is part of geoapps-utils package. ' # ' diff --git a/geoapps_utils/utils/transformations.py b/geoapps_utils/utils/transformations.py index f6a9c5ac..95a37950 100644 --- a/geoapps_utils/utils/transformations.py +++ b/geoapps_utils/utils/transformations.py @@ -1,5 +1,5 @@ # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2023-2025 Mira Geoscience Ltd. ' +# Copyright (c) 2022-2026 Mira Geoscience Ltd. ' # ' # This file is part of geoapps-utils package. ' # ' @@ -8,125 +8,152 @@ # ' # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' +from collections.abc import Sequence + import numpy as np import scipy.sparse as ssp from geoh5py.objects import Surface -def z_rotation_matrix(angle: np.ndarray | float) -> ssp.csr_matrix: +def apply_rotation(operator: np.ndarray, points: np.ndarray) -> np.ndarray: """ - Sparse matrix for heterogeneous vector rotation about the z axis. + Apply heterogeneous or homogeneous rotation matrix to point set. - To be used in a matrix-vector product with an array of shape (n, 3) - where n is the number of 3-component vectors. + :param operator: Rotation matrix of shape (3, 3) or (n, 3, 3) where n is the number of + points. + :param points: Array of shape (n, 3) representing the x, y, z coordinates to be rotated. + """ - :param angle: Array of angles in radians for counterclockwise rotation - about the z-axis. + if operator.shape[0] == points.shape[1]: + return operator.dot(points.T).T + + if operator.shape[0] == np.prod(points.shape): + vec = points.copy() + vec = vec.flatten() + vec = operator.dot(vec.T).T + vec = vec.reshape(-1, 3) + return vec + + raise ValueError( + f"Shape mismatch between operator ({operator.shape}) and points " + f"({points.shape}). For n points, the rotation operator should be size n " + "if homogeneous, or 3n if heterogeneous." + ) + + +def azimuth_to_unit_vector(azimuth: float) -> np.ndarray: """ + Convert an azimuth to a unit vector. - if isinstance(angle, np.ndarray): - n = len(angle) - rza = np.c_[np.cos(angle), np.cos(angle), np.ones(n)].T - rza = rza.flatten(order="F") - rzb = np.c_[np.sin(angle), np.zeros(n), np.zeros(n)].T - rzb = rzb.flatten(order="F") - rzc = np.c_[-np.sin(angle), np.zeros(n), np.zeros(n)].T - rzc = rzc.flatten(order="F") - rot_z = ssp.diags([rzb[:-1], rza, rzc[:-1]], [-1, 0, 1]) - else: - rot_z = np.r_[ - np.c_[np.cos(angle), -np.sin(angle), 0], - np.c_[np.sin(angle), np.cos(angle), 0], - np.c_[0, 0, 1], - ] + :param azimuth: Azimuth in degrees from north (0 to 360). + :return: Unit vector in the direction of the azimuth. + """ + theta = np.deg2rad(azimuth) + mat_z = np.r_[ + np.c_[np.cos(theta), -np.sin(theta), 0.0], + np.c_[np.sin(theta), np.cos(theta), 0.0], + np.c_[0.0, 0.0, 1.0], + ] + return np.array([0.0, 1.0, 0.0]).dot(mat_z) - return rot_z +def cartesian_normal_to_direction_and_dip(normals: np.ndarray) -> np.ndarray: + """ + Convert 3D normal vectors to dip and direction. -def x_rotation_matrix(angle: np.ndarray | float) -> ssp.csr_matrix: + :param normals: Array of shape (n, 3) representing the x, y, z components of a + normal vector in 3D space. + + :returns: Array of shape (n, 2) representing azimuth from 0 to 2pi radians + clockwise from north as viewed from above and dip from -pi to pi in positive + radians below the horizon and negative above. """ - Sparse matrix for heterogeneous vector rotation about the x axis. - To be used in a matrix-vector product with an array of shape (n, 3) - where n is the number of 3-component vectors. + spherical_normals = cartesian_to_spherical(normals) + direction_and_dip = spherical_normal_to_direction_and_dip(spherical_normals[:, 1:]) - :param angle: Array of angles in radians for counterclockwise rotation - about the x-axis. + return direction_and_dip + + +def cartesian_to_polar( + locations: np.ndarray, origin: Sequence = (0, 0, 0) +) -> np.ndarray: """ + Convert Cartesian coordinates (x, y, z) to polar coordinates defined as + (distance, azimuth, height). - if isinstance(angle, np.ndarray): - n = len(angle) - rxa = np.c_[np.ones(n), np.cos(angle), np.cos(angle)].T - rxa = rxa.flatten(order="F") - rxb = np.c_[np.zeros(n), np.sin(angle), np.zeros(n)].T - rxb = rxb.flatten(order="F") - rxc = np.c_[np.zeros(n), -np.sin(angle), np.zeros(n)].T - rxc = rxc.flatten(order="F") - rot_x = ssp.diags([rxb[:-1], rxa, rxc[:-1]], [-1, 0, 1]) - else: - rot_x = np.r_[ - np.c_[1, 0, 0], - np.c_[0, np.cos(angle), -np.sin(angle)], - np.c_[0, np.sin(angle), np.cos(angle)], - ] + Distances are relative to the origin, azimuth angles are + defined clockwise from North in degree and heights are the z-coordinates. - return rot_x + :param locations: Cartesian coordinates. + :param origin: Origin for distance and azimuth calculations. + + :return: Array of polar coordinates (distance, azimuth, height). + """ + if not isinstance(origin, Sequence | np.ndarray) or len(origin) != 3: + raise ValueError("Origin must be an iterable of length 3.") + local_xyz = locations - np.asarray(origin) + distances = np.linalg.norm(local_xyz[:, :2], axis=1) + azimuths = 90 - np.rad2deg(np.arctan2(local_xyz[:, 0], local_xyz[:, 1])) -def y_rotation_matrix(angle: np.ndarray | float) -> ssp.csr_matrix: + # Deal with points close to the origin with nearest neighbor azimuth + if np.any(distances < 1e-8): + indices = np.where(distances >= 1e-8)[0] + nearest = indices[np.argsort(distances[indices])] + azimuths[distances < 1e-8] = azimuths[nearest[0]] + + return np.c_[distances, azimuths, locations[:, 2]] + + +def cartesian_to_spherical(points: np.ndarray) -> np.ndarray: """ - Sparse matrix for heterogeneous vector rotation about the y axis. + Convert cartesian to spherical coordinate system. - To be used in a matrix-vector product with an array of shape (n, 3) - where n is the number of 3-component vectors. + :param points: Array of shape (n, 3) representing x, y, z coordinates of a point + in 3D space. - :param angle: Array of angles in radians for counterclockwise rotation - about the y-axis. + :returns: Array of shape (n, 3) representing the magnitude, azimuth and inclination + in spherical coordinates. The azimuth angle is measured in radians + counterclockwise from east in the range of 0 to 2pi as viewed from above, and + inclination angle is measured in radians from the positive z-axis. """ - if isinstance(angle, np.ndarray): - n = len(angle) - rxa = np.c_[np.cos(angle), np.ones(n), np.cos(angle)].T - rxa = rxa.flatten(order="F") - rxb = np.c_[-np.sin(angle), np.zeros(n), np.zeros(n)].T - rxb = rxb.flatten(order="F") - rxc = np.c_[np.sin(angle), np.zeros(n), np.zeros(n)].T - rxc = rxc.flatten(order="F") - rot_y = ssp.diags([rxb[:-2], rxa, rxc[:-2]], [-2, 0, 2]) - else: - rot_y = np.r_[ - np.c_[np.cos(angle), 0, np.sin(angle)], - np.c_[0, 1, 0], - np.c_[-np.sin(angle), 0, np.cos(angle)], - ] + magnitude = np.linalg.norm(points, axis=1) + inclination = np.arccos(points[:, 2] / magnitude) + azimuth = np.arctan2(points[:, 1], points[:, 0]) + return np.column_stack((magnitude, azimuth, inclination)) - return rot_y +def ccw_east_to_cw_north(azimuth: np.ndarray) -> np.ndarray: + """Convert counterclockwise azimuth from east to clockwise from north.""" + return (((5 * np.pi) / 2) - azimuth) % (2 * np.pi) -def apply_rotation(operator: np.ndarray, points: np.ndarray) -> np.ndarray: + +def compute_normals(surface: Surface) -> np.ndarray: """ - Apply heterogeneous or homogeneous rotation matrix to point set. + Compute normals for each triangle in a surface. - :param operator: Rotation matrix of shape (3, 3) or (n, 3, 3) where n is the number of - points. - :param points: Array of shape (n, 3) representing the x, y, z coordinates to be rotated. + :param surface: Surface object containing vertices and cells. + + :returns: Array of shape (n, 3) representing the normals for + each cell in the mesh. """ - if operator.shape[0] == points.shape[1]: - return operator.dot(points.T).T + vertices = surface.vertices + cells = surface.cells - if operator.shape[0] == np.prod(points.shape): - vec = points.copy() - vec = vec.flatten() - vec = operator.dot(vec.T).T - vec = vec.reshape(-1, 3) - return vec + v1 = vertices[cells[:, 1]] - vertices[cells[:, 0]] + v2 = vertices[cells[:, 2]] - vertices[cells[:, 0]] + normals = np.cross(v1, v2) + normals = normals / np.linalg.norm(normals, axis=1)[:, np.newaxis] + + return normals - raise ValueError( - f"Shape mismatch between operator ({operator.shape}) and points " - f"({points.shape}). For n points, the rotation operator should be size n " - "if homogeneous, or 3n if heterogeneous." - ) + +def inclination_to_dip(inclination: np.ndarray) -> np.ndarray: + """Convert inclination from positive z-axis to dip from horizon.""" + return inclination - (np.pi / 2) def rotate_points( @@ -185,35 +212,6 @@ def rotate_xyz(xyz: np.ndarray, center: list, theta: float, phi: float = 0.0): return xyz_out -def ccw_east_to_cw_north(azimuth: np.ndarray) -> np.ndarray: - """Convert counterclockwise azimuth from east to clockwise from north.""" - return (((5 * np.pi) / 2) - azimuth) % (2 * np.pi) - - -def inclination_to_dip(inclination: np.ndarray) -> np.ndarray: - """Convert inclination from positive z-axis to dip from horizon.""" - return inclination - (np.pi / 2) - - -def cartesian_to_spherical(points: np.ndarray) -> np.ndarray: - """ - Convert cartesian to spherical coordinate system. - - :param points: Array of shape (n, 3) representing x, y, z coordinates of a point - in 3D space. - - :returns: Array of shape (n, 3) representing the magnitude, azimuth and inclination - in spherical coordinates. The azimuth angle is measured in radians - counterclockwise from east in the range of 0 to 2pi as viewed from above, and - inclination angle is measured in radians from the positive z-axis. - """ - - magnitude = np.linalg.norm(points, axis=1) - inclination = np.arccos(points[:, 2] / magnitude) - azimuth = np.arctan2(points[:, 1], points[:, 0]) - return np.column_stack((magnitude, azimuth, inclination)) - - def spherical_normal_to_direction_and_dip(angles: np.ndarray) -> np.ndarray: """ Convert normals in spherical coordinates to dip and direction of the tangent plane. @@ -239,56 +237,91 @@ def spherical_normal_to_direction_and_dip(angles: np.ndarray) -> np.ndarray: ) -def cartesian_normal_to_direction_and_dip(normals: np.ndarray) -> np.ndarray: +def x_rotation_matrix(angle: np.ndarray | float) -> ssp.csr_matrix: """ - Convert 3D normal vectors to dip and direction. + Sparse matrix for heterogeneous vector rotation about the x axis. - :param normals: Array of shape (n, 3) representing the x, y, z components of a - normal vector in 3D space. + To be used in a matrix-vector product with an array of shape (n, 3) + where n is the number of 3-component vectors. - :returns: Array of shape (n, 2) representing azimuth from 0 to 2pi radians - clockwise from north as viewed from above and dip from -pi to pi in positive - radians below the horizon and negative above. + :param angle: Array of angles in radians for counterclockwise rotation + about the x-axis. """ - spherical_normals = cartesian_to_spherical(normals) - direction_and_dip = spherical_normal_to_direction_and_dip(spherical_normals[:, 1:]) + if isinstance(angle, np.ndarray): + n = len(angle) + rxa = np.c_[np.ones(n), np.cos(angle), np.cos(angle)].T + rxa = rxa.flatten(order="F") + rxb = np.c_[np.zeros(n), np.sin(angle), np.zeros(n)].T + rxb = rxb.flatten(order="F") + rxc = np.c_[np.zeros(n), -np.sin(angle), np.zeros(n)].T + rxc = rxc.flatten(order="F") + rot_x = ssp.diags([rxb[:-1], rxa, rxc[:-1]], [-1, 0, 1]) + else: + rot_x = np.r_[ + np.c_[1, 0, 0], + np.c_[0, np.cos(angle), -np.sin(angle)], + np.c_[0, np.sin(angle), np.cos(angle)], + ] - return direction_and_dip + return rot_x -def compute_normals(surface: Surface) -> np.ndarray: +def y_rotation_matrix(angle: np.ndarray | float) -> ssp.csr_matrix: """ - Compute normals for each triangle in a surface. + Sparse matrix for heterogeneous vector rotation about the y axis. - :param surface: Surface object containing vertices and cells. + To be used in a matrix-vector product with an array of shape (n, 3) + where n is the number of 3-component vectors. - :returns: Array of shape (n, 3) representing the normals for - each cell in the mesh. + :param angle: Array of angles in radians for counterclockwise rotation + about the y-axis. """ - vertices = surface.vertices - cells = surface.cells - - v1 = vertices[cells[:, 1]] - vertices[cells[:, 0]] - v2 = vertices[cells[:, 2]] - vertices[cells[:, 0]] - normals = np.cross(v1, v2) - normals = normals / np.linalg.norm(normals, axis=1)[:, np.newaxis] + if isinstance(angle, np.ndarray): + n = len(angle) + rxa = np.c_[np.cos(angle), np.ones(n), np.cos(angle)].T + rxa = rxa.flatten(order="F") + rxb = np.c_[-np.sin(angle), np.zeros(n), np.zeros(n)].T + rxb = rxb.flatten(order="F") + rxc = np.c_[np.sin(angle), np.zeros(n), np.zeros(n)].T + rxc = rxc.flatten(order="F") + rot_y = ssp.diags([rxb[:-2], rxa, rxc[:-2]], [-2, 0, 2]) + else: + rot_y = np.r_[ + np.c_[np.cos(angle), 0, np.sin(angle)], + np.c_[0, 1, 0], + np.c_[-np.sin(angle), 0, np.cos(angle)], + ] - return normals + return rot_y -def azimuth_to_unit_vector(azimuth: float) -> np.ndarray: +def z_rotation_matrix(angle: np.ndarray | float) -> ssp.csr_matrix: """ - Convert an azimuth to a unit vector. + Sparse matrix for heterogeneous vector rotation about the z axis. - :param azimuth: Azimuth in degrees from north (0 to 360). - :return: Unit vector in the direction of the azimuth. + To be used in a matrix-vector product with an array of shape (n, 3) + where n is the number of 3-component vectors. + + :param angle: Array of angles in radians for counterclockwise rotation + about the z-axis. """ - theta = np.deg2rad(azimuth) - mat_z = np.r_[ - np.c_[np.cos(theta), -np.sin(theta), 0.0], - np.c_[np.sin(theta), np.cos(theta), 0.0], - np.c_[0.0, 0.0, 1.0], - ] - return np.array([0.0, 1.0, 0.0]).dot(mat_z) + + if isinstance(angle, np.ndarray): + n = len(angle) + rza = np.c_[np.cos(angle), np.cos(angle), np.ones(n)].T + rza = rza.flatten(order="F") + rzb = np.c_[np.sin(angle), np.zeros(n), np.zeros(n)].T + rzb = rzb.flatten(order="F") + rzc = np.c_[-np.sin(angle), np.zeros(n), np.zeros(n)].T + rzc = rzc.flatten(order="F") + rot_z = ssp.diags([rzb[:-1], rza, rzc[:-1]], [-1, 0, 1]) + else: + rot_z = np.r_[ + np.c_[np.cos(angle), -np.sin(angle), 0], + np.c_[np.sin(angle), np.cos(angle), 0], + np.c_[0, 0, 1], + ] + + return rot_z diff --git a/geoapps_utils/utils/workspace.py b/geoapps_utils/utils/workspace.py index 27b9a28e..9f78c7dd 100644 --- a/geoapps_utils/utils/workspace.py +++ b/geoapps_utils/utils/workspace.py @@ -1,5 +1,5 @@ # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2023-2025 Mira Geoscience Ltd. ' +# Copyright (c) 2022-2026 Mira Geoscience Ltd. ' # ' # This file is part of geoapps-utils package. ' # ' diff --git a/package.rst b/package.rst index 3879834e..cc969286 100644 --- a/package.rst +++ b/package.rst @@ -29,7 +29,7 @@ License ^^^^^^^ MIT License -Copyright (c) 2022-2025 Mira Geoscience +Copyright (c) 2022-2026 Mira Geoscience Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal diff --git a/tests/__init__.py b/tests/__init__.py index 348e9ed6..02fcb1b2 100644 --- a/tests/__init__.py +++ b/tests/__init__.py @@ -1,5 +1,5 @@ # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2022-2025 Mira Geoscience Ltd. ' +# Copyright (c) 2022-2026 Mira Geoscience Ltd. ' # ' # This file is part of geoapps-utils package. ' # ' diff --git a/tests/conversions_test.py b/tests/conversions_test.py index f87f1aec..b8a1b63d 100644 --- a/tests/conversions_test.py +++ b/tests/conversions_test.py @@ -1,5 +1,5 @@ # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2023-2025 Mira Geoscience Ltd. ' +# Copyright (c) 2022-2026 Mira Geoscience Ltd. ' # ' # This file is part of geoapps-utils package. ' # ' diff --git a/tests/dataclass_test.py b/tests/dataclass_test.py index be3ebab9..fe0366eb 100644 --- a/tests/dataclass_test.py +++ b/tests/dataclass_test.py @@ -1,5 +1,5 @@ # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2023-2025 Mira Geoscience Ltd. ' +# Copyright (c) 2022-2026 Mira Geoscience Ltd. ' # ' # This file is part of geoapps-utils package. ' # ' diff --git a/tests/driver_test.py b/tests/driver_test.py index f4aa9459..b174ec71 100644 --- a/tests/driver_test.py +++ b/tests/driver_test.py @@ -1,5 +1,5 @@ # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2023-2025 Mira Geoscience Ltd. ' +# Copyright (c) 2022-2026 Mira Geoscience Ltd. ' # ' # This file is part of geoapps-utils package. ' # ' diff --git a/tests/dummy_driver_test.py b/tests/dummy_driver_test.py index 441becf0..8fb22f8c 100644 --- a/tests/dummy_driver_test.py +++ b/tests/dummy_driver_test.py @@ -1,5 +1,5 @@ # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2025 Mira Geoscience Ltd. ' +# Copyright (c) 2022-2026 Mira Geoscience Ltd. ' # ' # This file is part of geoapps-utils package. ' # ' diff --git a/tests/formatters_test.py b/tests/formatters_test.py index bada8f76..63e959ef 100644 --- a/tests/formatters_test.py +++ b/tests/formatters_test.py @@ -1,5 +1,5 @@ # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2023-2025 Mira Geoscience Ltd. ' +# Copyright (c) 2022-2026 Mira Geoscience Ltd. ' # ' # This file is part of geoapps-utils package. ' # ' diff --git a/tests/importing_test.py b/tests/importing_test.py index 943ef480..f215dfe2 100644 --- a/tests/importing_test.py +++ b/tests/importing_test.py @@ -1,5 +1,5 @@ # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2023-2025 Mira Geoscience Ltd. ' +# Copyright (c) 2022-2026 Mira Geoscience Ltd. ' # ' # This file is part of geoapps-utils package. ' # ' diff --git a/tests/iterables_test.py b/tests/iterables_test.py index 09588740..925b0703 100644 --- a/tests/iterables_test.py +++ b/tests/iterables_test.py @@ -1,5 +1,5 @@ # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2023-2025 Mira Geoscience Ltd. ' +# Copyright (c) 2022-2026 Mira Geoscience Ltd. ' # ' # This file is part of geoapps-utils package. ' # ' diff --git a/tests/locations_test.py b/tests/locations_test.py index 52195b68..eb062844 100644 --- a/tests/locations_test.py +++ b/tests/locations_test.py @@ -1,5 +1,5 @@ # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2023-2025 Mira Geoscience Ltd. ' +# Copyright (c) 2022-2026 Mira Geoscience Ltd. ' # ' # This file is part of geoapps-utils package. ' # ' diff --git a/tests/numerical_test.py b/tests/numerical_test.py index a350c08b..d1764f57 100644 --- a/tests/numerical_test.py +++ b/tests/numerical_test.py @@ -1,5 +1,5 @@ # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2023-2025 Mira Geoscience Ltd. ' +# Copyright (c) 2022-2026 Mira Geoscience Ltd. ' # ' # This file is part of geoapps-utils package. ' # ' @@ -17,6 +17,7 @@ from geoapps_utils.utils.numerical import ( fibonacci_series, fit_circle, + inverse_weighted_operator, running_mean, traveling_salesman, weighted_average, @@ -196,3 +197,26 @@ def test_fit_cicle(): def test_fibonacci_series(): np.testing.assert_allclose(fibonacci_series(1), [0]) np.testing.assert_allclose(fibonacci_series(6), [0, 1, 1, 2, 3, 5]) + + +def test_inverse_weighted_operator(): + """ + Test the inverse_weighted_operator utility function. + + For a constant input, the output should be the same constant. + """ + power = 2.0 + threshold = 1e-12 + shape = (100, 1000) + values = np.random.randn(shape[0] * 2) + indices = np.c_[ + np.random.randint(0, shape[1] - 1, shape[0]), + np.random.randint(0, shape[1] - 1, shape[0]), + ].flatten() + + opt = inverse_weighted_operator(values, indices, shape, power, threshold) + test_val = np.random.randn() + interp = opt * np.full(shape[1], test_val) + + assert opt.shape == shape + np.testing.assert_allclose(interp, test_val, rtol=1e-3) diff --git a/tests/param_sweeps/__init__.py b/tests/param_sweeps/__init__.py index 5197aa21..02fcb1b2 100644 --- a/tests/param_sweeps/__init__.py +++ b/tests/param_sweeps/__init__.py @@ -1,5 +1,5 @@ # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2023-2025 Mira Geoscience Ltd. ' +# Copyright (c) 2022-2026 Mira Geoscience Ltd. ' # ' # This file is part of geoapps-utils package. ' # ' diff --git a/tests/param_sweeps/driver_test.py b/tests/param_sweeps/driver_test.py index d7b79bd0..e124c075 100644 --- a/tests/param_sweeps/driver_test.py +++ b/tests/param_sweeps/driver_test.py @@ -1,5 +1,5 @@ # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2023-2025 Mira Geoscience Ltd. ' +# Copyright (c) 2022-2026 Mira Geoscience Ltd. ' # ' # This file is part of geoapps-utils package. ' # ' diff --git a/tests/param_sweeps/generate_test.py b/tests/param_sweeps/generate_test.py index 274a079f..8fbeb9db 100644 --- a/tests/param_sweeps/generate_test.py +++ b/tests/param_sweeps/generate_test.py @@ -1,5 +1,5 @@ # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2023-2025 Mira Geoscience Ltd. ' +# Copyright (c) 2022-2026 Mira Geoscience Ltd. ' # ' # This file is part of geoapps-utils package. ' # ' diff --git a/tests/plates_test.py b/tests/plates_test.py index 9cc809ac..9e5bcb5b 100644 --- a/tests/plates_test.py +++ b/tests/plates_test.py @@ -1,5 +1,5 @@ # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2023-2025 Mira Geoscience Ltd. ' +# Copyright (c) 2022-2026 Mira Geoscience Ltd. ' # ' # This file is part of geoapps-utils package. ' # ' diff --git a/tests/plotting_test.py b/tests/plotting_test.py index 3dda4c7b..9b46abe1 100644 --- a/tests/plotting_test.py +++ b/tests/plotting_test.py @@ -1,5 +1,5 @@ # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2023-2025 Mira Geoscience Ltd. ' +# Copyright (c) 2022-2026 Mira Geoscience Ltd. ' # ' # This file is part of geoapps-utils package. ' # ' diff --git a/tests/reference_test.py b/tests/reference_test.py index 86d8d2f2..b8ca6dac 100644 --- a/tests/reference_test.py +++ b/tests/reference_test.py @@ -1,5 +1,5 @@ # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2024-2025 Mira Geoscience Ltd. ' +# Copyright (c) 2024-2026 Mira Geoscience Ltd. ' # ' # This file is part of geoapps-utils package. ' # ' diff --git a/tests/transformations_test.py b/tests/transformations_test.py index 62becebf..9104da43 100644 --- a/tests/transformations_test.py +++ b/tests/transformations_test.py @@ -1,5 +1,5 @@ # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2023-2025 Mira Geoscience Ltd. ' +# Copyright (c) 2022-2026 Mira Geoscience Ltd. ' # ' # This file is part of geoapps-utils package. ' # ' @@ -19,6 +19,7 @@ from geoapps_utils.utils.transformations import ( azimuth_to_unit_vector, cartesian_normal_to_direction_and_dip, + cartesian_to_polar, cartesian_to_spherical, compute_normals, rotate_xyz, @@ -187,3 +188,42 @@ def test_azimuth_to_unit_vector(): assert np.allclose(azimuth_to_unit_vector(180.0), np.array([0.0, -1.0, 0.0])) assert np.allclose(azimuth_to_unit_vector(270.0), np.array([-1.0, 0.0, 0.0])) assert np.allclose(azimuth_to_unit_vector(360.0), np.array([0.0, 1.0, 0.0])) + + +def test_cartesian_to_polar(): + """ + Test the xyz_to_polar utility function. + """ + + rad = np.arange(0, 100, 10) + azm = np.pi / 3.0 + + # Create x, y, z coordinates + x = rad * np.cos(azm) + y = rad * np.sin(azm) + z = np.random.randn(x.shape[0]) # Random z values + locations = np.c_[x, y, z] + + # Start reference + polar = cartesian_to_polar(locations) + np.testing.assert_almost_equal(polar[0, 0], 0.0) # First point at zero distance + np.testing.assert_allclose( + polar[:, 1], np.rad2deg(azm) + ) # All other distances positive + + with pytest.raises(ValueError, match="Origin must be an iterable of length 3."): + _ = cartesian_to_polar(locations, origin=(5.0, "abc")) + + # Mean reference locations + polar = cartesian_to_polar(locations, origin=np.mean(locations, axis=0)) + np.testing.assert_almost_equal( + polar[4::-1, 0], polar[5:, 0] + ) # Opposite side of center + np.testing.assert_almost_equal(polar[4::-1, 0], polar[5:, 0]) # Polar opposite + np.testing.assert_allclose(z, polar[:, 2]) # Preserves z + + # End reference + polar = cartesian_to_polar(locations, origin=locations[-1, :]) + np.testing.assert_allclose( + polar[:, 1], np.rad2deg(azm) + 180 + ) # All other distances positive diff --git a/tests/uijson_run_test.py b/tests/uijson_run_test.py index 6509a33f..00ae0d88 100644 --- a/tests/uijson_run_test.py +++ b/tests/uijson_run_test.py @@ -1,5 +1,5 @@ # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2025 Mira Geoscience Ltd. ' +# Copyright (c) 2022-2026 Mira Geoscience Ltd. ' # ' # This file is part of geoapps-utils package. ' # ' diff --git a/tests/version_test.py b/tests/version_test.py index 5a4e1c0d..c61d4d6a 100644 --- a/tests/version_test.py +++ b/tests/version_test.py @@ -1,5 +1,5 @@ # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2022-2025 Mira Geoscience Ltd. ' +# Copyright (c) 2022-2026 Mira Geoscience Ltd. ' # ' # This file is part of geoapps-utils package. ' # '