From 896b18655fc690a27c354dace49c67904f05c8c0 Mon Sep 17 00:00:00 2001 From: benjamink Date: Thu, 31 Jul 2025 14:32:03 -0700 Subject: [PATCH 1/5] refactor for simpeg-drivers --- geoapps_utils/modelling/halfspace.py | 9 + .../modelling/plate_in_a_halfspace.py | 222 ++++++++++++++++++ geoapps_utils/modelling/surveys/__init__.py | 9 + geoapps_utils/modelling/surveys/dcip.py | 75 ++++++ .../surveys/frequency_domain/__init__.py | 15 ++ .../surveys/frequency_domain/fdem.py | 59 +++++ .../modelling/surveys/time_domain/__init__.py | 23 ++ .../surveys/time_domain/airborne_tdem.py | 42 ++++ .../surveys/time_domain/ground_tdem.py | 100 ++++++++ geoapps_utils/utils/locations.py | 35 ++- 10 files changed, 588 insertions(+), 1 deletion(-) create mode 100644 geoapps_utils/modelling/halfspace.py create mode 100644 geoapps_utils/modelling/plate_in_a_halfspace.py create mode 100644 geoapps_utils/modelling/surveys/__init__.py create mode 100644 geoapps_utils/modelling/surveys/dcip.py create mode 100644 geoapps_utils/modelling/surveys/frequency_domain/__init__.py create mode 100644 geoapps_utils/modelling/surveys/frequency_domain/fdem.py create mode 100644 geoapps_utils/modelling/surveys/time_domain/__init__.py create mode 100644 geoapps_utils/modelling/surveys/time_domain/airborne_tdem.py create mode 100644 geoapps_utils/modelling/surveys/time_domain/ground_tdem.py diff --git a/geoapps_utils/modelling/halfspace.py b/geoapps_utils/modelling/halfspace.py new file mode 100644 index 00000000..76e3205e --- /dev/null +++ b/geoapps_utils/modelling/halfspace.py @@ -0,0 +1,9 @@ +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' +# Copyright (c) 2025 Mira Geoscience Ltd. ' +# ' +# This file is part of geoapps-utils package. ' +# ' +# geoapps-utils is distributed under the terms and conditions of the MIT License ' +# (see LICENSE file at the root of this source code package). ' +# ' +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' diff --git a/geoapps_utils/modelling/plate_in_a_halfspace.py b/geoapps_utils/modelling/plate_in_a_halfspace.py new file mode 100644 index 00000000..b53ddaf4 --- /dev/null +++ b/geoapps_utils/modelling/plate_in_a_halfspace.py @@ -0,0 +1,222 @@ +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' +# Copyright (c) 2025 Mira Geoscience Ltd. ' +# ' +# This file is part of geoapps-utils package. ' +# ' +# geoapps-utils is distributed under the terms and conditions of the MIT License ' +# (see LICENSE file at the root of this source code package). ' +# ' +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' + +import numpy as np +from pathlib import Path +from scipy.spatial import Delaunay + +from geoh5py import Workspace +from geoh5py.objects import Points, Surface + + +from geoapps_utils.modelling.plates import PlateModel, make_plate +from tests.testing_utils.terrain import gaussian_topo_drape + + + + + + +plate_model_default = PlateModel( + strike_length=40.0, + dip_length=40.0, + width=40.0, + origin=(0.0, 0.0, 10.0), +) + + +def setup_inversion_workspace( + work_dir, + plate_model: PlateModel = plate_model_default, + background=None, + anomaly=None, + cell_size=(5.0, 5.0, 5.0), + center=(0.0, 0.0, 0.0), + n_electrodes=20, + n_lines=5, + refinement=(4, 6), + x_limits=(-100.0, 100.0), + y_limits=(-100.0, 100.0), + padding_distance=100, + drape_height=5.0, + inversion_type="other", + flatten=False, + geoh5=None, +): + """ + Creates a workspace populated with objects to simulate/invert a simple model. + + rot_xyz: Rotation angles in degrees about the x, y, and z axes. + """ + filepath = Path(work_dir) / "inversion_test.ui.geoh5" + if geoh5 is None: + if filepath.is_file(): + filepath.unlink() + geoh5 = Workspace.create(filepath) + # Topography + xx, yy = np.meshgrid(np.linspace(-200.0, 200.0, 50), np.linspace(-200.0, 200.0, 50)) + + zz = gaussian_topo_drape(xx, yy, flatten=flatten) + topo = np.c_[ + utils.mkvc(xx) + center[0], + utils.mkvc(yy) + center[1], + utils.mkvc(zz) + center[2], + ] + triang = Delaunay(topo[:, :2]) + topography = Surface.create( + geoh5, + vertices=topo, + cells=triang.simplices, # pylint: disable=E1101 + name="topography", + ) + + # Observation points + n_electrodes = ( + 4 + if ("polarization" in inversion_type or "current" in inversion_type) + & (n_electrodes < 4) + else n_electrodes + ) + xr = np.linspace(x_limits[0], x_limits[1], int(n_electrodes)) + yr = np.linspace(y_limits[0], y_limits[1], int(n_lines)) + X, Y = np.meshgrid(xr, yr) + Z = gaussian_topo_drape(X, Y, flatten=flatten) + drape_height + + vertices = np.c_[ + utils.mkvc(X.T) + center[0], + utils.mkvc(Y.T) + center[1], + utils.mkvc(Z.T) + center[2], + ] + + if "polarization" in inversion_type or "current" in inversion_type: + survey = generate_dc_survey(geoh5, X, Y, Z) + + elif "magnetotellurics" in inversion_type: + survey = MTReceivers.create( + geoh5, + vertices=vertices, + name="survey", + components=[ + "Zxx (real)", + "Zxx (imag)", + "Zxy (real)", + "Zxy (imag)", + "Zyx (real)", + "Zyx (imag)", + "Zyy (real)", + "Zyy (imag)", + ], + channels=[10.0, 100.0, 1000.0], + ) + + elif "tipper" in inversion_type: + survey = TipperReceivers.create( + geoh5, + vertices=vertices, + name="survey", + components=[ + "Txz (real)", + "Txz (imag)", + "Tyz (real)", + "Tyz (imag)", + ], + ) + survey.base_stations = TipperBaseStations.create( + geoh5, vertices=np.c_[vertices[0, :]].T + ) + survey.channels = [10.0, 100.0, 1000.0] + dist = np.linalg.norm( + survey.vertices[survey.cells[:, 0], :] + - survey.vertices[survey.cells[:, 1], :], + axis=1, + ) + # survey.cells = survey.cells[dist < 100.0, :] + survey.remove_cells(np.where(dist > 200)[0]) + + elif inversion_type in ["fdem", "fem", "fdem 1d"]: + survey = generate_fdem_survey(geoh5, vertices) + + elif "tdem" in inversion_type: + survey = generate_tdem_survey( + geoh5, + vertices, + n_lines, + flatten=flatten, + airborne="airborne" in inversion_type, + ) + + else: + survey = Points.create( + geoh5, + vertices=vertices, + name="survey", + ) + + # Create a mesh + if "2d" in inversion_type: + lines = survey.get_entity("line_ids")[0].values + entity, mesh, _ = get_drape_model( # pylint: disable=unbalanced-tuple-unpacking + geoh5, + "Models", + survey.vertices[np.unique(survey.cells[lines == 101, :]), :], + [cell_size[0], cell_size[2]], + 100.0, + [padding_distance] * 2 + [padding_distance] * 2, + 1.1, + parent=None, + return_colocated_mesh=True, + return_sorting=True, + ) + active = active_from_xyz(entity, topography.vertices, grid_reference="top") + + else: + padDist = np.ones((3, 2)) * padding_distance + mesh = mesh_builder_xyz( + vertices - np.r_[cell_size] / 2.0, + cell_size, + depth_core=100.0, + padding_distance=padDist, + mesh_type="TREE", + tree_diagonal_balance=False, + ) + mesh = OctreeDriver.refine_tree_from_surface( + mesh, + topography, + levels=refinement, + finalize=False, + ) + + if inversion_type in ["fdem", "airborne tdem"]: + mesh = OctreeDriver.refine_tree_from_points( + mesh, + vertices, + levels=[2], + finalize=False, + ) + + mesh.finalize() + entity = treemesh_2_octree(geoh5, mesh, name="mesh") + active = active_from_xyz(entity, topography.vertices, grid_reference="top") + + # Create the Model + cell_centers = entity.centroids.copy() + + model = make_plate( + cell_centers, plate_model, background=background, anomaly=anomaly + ) + + if "1d" in inversion_type: + model = background * np.ones(mesh.nC) + model[(mesh.cell_centers[:, 2] < 0) & (mesh.cell_centers[:, 2] > -20)] = anomaly + + model[~active] = np.nan + model = entity.add_data({"model": {"values": model}}) + geoh5.close() + return geoh5, entity, model, survey, topography diff --git a/geoapps_utils/modelling/surveys/__init__.py b/geoapps_utils/modelling/surveys/__init__.py new file mode 100644 index 00000000..76e3205e --- /dev/null +++ b/geoapps_utils/modelling/surveys/__init__.py @@ -0,0 +1,9 @@ +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' +# Copyright (c) 2025 Mira Geoscience Ltd. ' +# ' +# This file is part of geoapps-utils package. ' +# ' +# geoapps-utils is distributed under the terms and conditions of the MIT License ' +# (see LICENSE file at the root of this source code package). ' +# ' +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' diff --git a/geoapps_utils/modelling/surveys/dcip.py b/geoapps_utils/modelling/surveys/dcip.py new file mode 100644 index 00000000..ea5b60c7 --- /dev/null +++ b/geoapps_utils/modelling/surveys/dcip.py @@ -0,0 +1,75 @@ +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' +# Copyright (c) 2025 Mira Geoscience Ltd. ' +# ' +# This file is part of geoapps-utils package. ' +# ' +# geoapps-utils is distributed under the terms and conditions of the MIT License ' +# (see LICENSE file at the root of this source code package). ' +# ' +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' + +import numpy as np +from geoh5py import Workspace +from geoh5py.objects import CurrentElectrode, PotentialElectrode + + +def generate_dc_survey( + workspace: Workspace, + x_loc: np.ndarray, + y_loc: np.ndarray, + z_loc: np.ndarray | None = None, +) -> PotentialElectrode: + """ + Utility function to generate a DC survey. + """ + # Create sources along line + if z_loc is None: + z_loc = np.zeros_like(x_loc) + + vertices = np.c_[x_loc.ravel(), y_loc.ravel(), z_loc.ravel()] + parts = np.kron(np.arange(x_loc.shape[0]), np.ones(x_loc.shape[1])).astype("int") + currents = CurrentElectrode.create(workspace, vertices=vertices, parts=parts) + currents.add_default_ab_cell_id() + n_dipoles = 9 + dipoles = [] + current_id = [] + + for val in currents.ab_cell_id.values: + cell_id = int(currents.ab_map[val]) - 1 + + for dipole in range(n_dipoles): + dipole_ids = currents.cells[cell_id, :] + 2 + dipole + + if ( + any(dipole_ids > (currents.n_vertices - 1)) + or len( + np.unique(parts[np.r_[currents.cells[cell_id, 0], dipole_ids[1]]]) + ) + > 1 + ): + continue + + dipoles += [dipole_ids] + current_id += [val] + + potentials = PotentialElectrode.create( + workspace, vertices=vertices, cells=np.vstack(dipoles).astype("uint32") + ) + line_id = potentials.vertices[potentials.cells[:, 0], 1] + line_id = (line_id - np.min(line_id) + 1).astype(np.int32) + line_reference = {0: "Unknown"} + line_reference.update({k: str(k) for k in np.unique(line_id)}) + potentials.add_data( + { + "line_ids": { + "values": line_id, + "type": "REFERENCED", + "value_map": line_reference, + } + } + ) + potentials.ab_cell_id = np.hstack(current_id).astype("int32") + potentials.current_electrodes = currents + currents.potential_electrodes = potentials + + return potentials \ No newline at end of file diff --git a/geoapps_utils/modelling/surveys/frequency_domain/__init__.py b/geoapps_utils/modelling/surveys/frequency_domain/__init__.py new file mode 100644 index 00000000..c50ec18f --- /dev/null +++ b/geoapps_utils/modelling/surveys/frequency_domain/__init__.py @@ -0,0 +1,15 @@ +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' +# Copyright (c) 2025 Mira Geoscience Ltd. ' +# ' +# This file is part of geoapps-utils package. ' +# ' +# geoapps-utils is distributed under the terms and conditions of the MIT License ' +# (see LICENSE file at the root of this source code package). ' +# ' +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' + +frequency_config = [ + {"Coaxial data": False, "Frequency": 900.0, "Offset": 7.86}, + {"Coaxial data": False, "Frequency": 7200.0, "Offset": 7.86}, + {"Coaxial data": False, "Frequency": 56000.0, "Offset": 6.3}, +] \ No newline at end of file diff --git a/geoapps_utils/modelling/surveys/frequency_domain/fdem.py b/geoapps_utils/modelling/surveys/frequency_domain/fdem.py new file mode 100644 index 00000000..4c2f4971 --- /dev/null +++ b/geoapps_utils/modelling/surveys/frequency_domain/fdem.py @@ -0,0 +1,59 @@ +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' +# Copyright (c) 2025 Mira Geoscience Ltd. ' +# ' +# This file is part of geoapps-utils package. ' +# ' +# geoapps-utils is distributed under the terms and conditions of the MIT License ' +# (see LICENSE file at the root of this source code package). ' +# ' +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' + +import numpy as np +from geoh5py import Workspace +from geoh5py.objects import AirborneFEMReceivers, AirborneFEMTransmitters +from geoapps_utils.utils.locations import mask_large_connections + +from . import frequency_config + +def generate_fdem_survey( + geoh5: Workspace, + vertices: np.ndarray +): + survey = AirborneFEMReceivers.create(geoh5, vertices=vertices, name="Airborne_rx") + + survey.metadata["EM Dataset"]["Frequency configurations"] = frequency_config + + tx_locs = [] + freqs = [] + for config in frequency_config: + tx_vertices = vertices.copy() + tx_vertices[:, 0] -= config["Offset"] + tx_locs.append(tx_vertices) + freqs.append([[config["Frequency"]] * len(vertices)]) + tx_locs = np.vstack(tx_locs) + freqs = np.hstack(freqs).flatten() + + transmitters = AirborneFEMTransmitters.create( + geoh5, vertices=tx_locs, name="Airborne_tx" + ) + survey.transmitters = transmitters + survey.channels = [k["Frequency"] for k in frequency_config] + + transmitters.add_data( + { + "Tx frequency": { + "values": freqs, + "association": "VERTEX", + "primitive_type": "REFERENCED", + "value_map": {k: str(k) for k in freqs}, + } + } + ) + + survey.remove_cells(mask_large_connections(survey, 200.0)) + transmitters.remove_cells(mask_large_connections(transmitters, 200.0)) + + + return survey + + diff --git a/geoapps_utils/modelling/surveys/time_domain/__init__.py b/geoapps_utils/modelling/surveys/time_domain/__init__.py new file mode 100644 index 00000000..cb618213 --- /dev/null +++ b/geoapps_utils/modelling/surveys/time_domain/__init__.py @@ -0,0 +1,23 @@ +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' +# Copyright (c) 2025 Mira Geoscience Ltd. ' +# ' +# This file is part of geoapps-utils package. ' +# ' +# geoapps-utils is distributed under the terms and conditions of the MIT License ' +# (see LICENSE file at the root of this source code package). ' +# ' +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' + +import numpy as np + +channels = np.r_[3e-04, 6e-04, 1.2e-03] * 1e3 +waveform = np.c_[ + np.r_[ + np.arange(-0.002, -0.0001, 5e-4), + np.arange(-0.0004, 0.0, 1e-4), + np.arange(0.0, 0.002, 5e-4), + ] + * 1e3 + + 2.0, + np.r_[np.linspace(0, 1, 4), np.linspace(0.9, 0.0, 4), np.zeros(4)], +] \ No newline at end of file diff --git a/geoapps_utils/modelling/surveys/time_domain/airborne_tdem.py b/geoapps_utils/modelling/surveys/time_domain/airborne_tdem.py new file mode 100644 index 00000000..e6236bf6 --- /dev/null +++ b/geoapps_utils/modelling/surveys/time_domain/airborne_tdem.py @@ -0,0 +1,42 @@ +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' +# Copyright (c) 2025 Mira Geoscience Ltd. ' +# ' +# This file is part of geoapps-utils package. ' +# ' +# geoapps-utils is distributed under the terms and conditions of the MIT License ' +# (see LICENSE file at the root of this source code package). ' +# ' +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' + +import numpy as np +from geoh5py import Workspace +from geoh5py.objects import ( + AirborneTEMReceivers, + AirborneTEMTransmitters, +) + +from geoapps_utils.utils.locations import mask_large_connections +from . import channels, waveform + +def generate_airborne_tdem_survey( + geoh5: Workspace, + vertices: np.ndarray, +): + + survey = AirborneTEMReceivers.create( + geoh5, vertices=vertices, name="Airborne_rx" + ) + transmitters = AirborneTEMTransmitters.create( + geoh5, vertices=vertices, name="Airborne_tx" + ) + mask = mask_large_connections(survey, 200.0) + survey.remove_cells(mask) + transmitters.remove_cells(mask) + + survey.transmitters = transmitters + survey.channels = channels + survey.waveform = waveform + survey.timing_mark = 2.0 + survey.unit = "Milliseconds (ms)" + + return survey \ No newline at end of file diff --git a/geoapps_utils/modelling/surveys/time_domain/ground_tdem.py b/geoapps_utils/modelling/surveys/time_domain/ground_tdem.py new file mode 100644 index 00000000..d3a3ead8 --- /dev/null +++ b/geoapps_utils/modelling/surveys/time_domain/ground_tdem.py @@ -0,0 +1,100 @@ +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' +# Copyright (c) 2025 Mira Geoscience Ltd. ' +# ' +# This file is part of geoapps-utils package. ' +# ' +# geoapps-utils is distributed under the terms and conditions of the MIT License ' +# (see LICENSE file at the root of this source code package). ' +# ' +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' + +import numpy as np + +from geoh5py import Workspace +from geoh5py.objects import ( + LargeLoopGroundTEMReceivers, + LargeLoopGroundTEMTransmitters, +) + +from tests.testing_utils.terrain import gaussian_topo_drape +from . import channels, waveform + +def generate_tdem_survey( + geoh5: Workspace, + vertices: np.ndarray, + n_lines: int, + flatten: bool = False, +): + + + X = vertices[:, 0].reshape((n_lines, -1)) + Y = vertices[:, 1].reshape((n_lines, -1)) + Z = vertices[:, 2].reshape((n_lines, -1)) + center = np.mean(vertices, axis=0) + arrays = [ + np.c_[ + X[: int(n_lines / 2), :].flatten(), + Y[: int(n_lines / 2), :].flatten(), + Z[: int(n_lines / 2), :].flatten(), + ], + np.c_[ + X[int(n_lines / 2) :, :].flatten(), + Y[int(n_lines / 2) :, :].flatten(), + Z[int(n_lines / 2) :, :].flatten(), + ], + ] + loops = [] + loop_cells = [] + loop_id = [] + count = 0 + for ind, array in enumerate(arrays): + loop_id += [np.ones(array.shape[0]) * (ind + 1)] + min_loc = np.min(array, axis=0) + max_loc = np.max(array, axis=0) + loop = np.vstack( + [ + np.c_[ + np.ones(5) * min_loc[0], + np.linspace(min_loc[1], max_loc[1], 5), + ], + np.c_[ + np.linspace(min_loc[0], max_loc[0], 5)[1:], + np.ones(4) * max_loc[1], + ], + np.c_[ + np.ones(4) * max_loc[0], + np.linspace(max_loc[1], min_loc[1], 5)[1:], + ], + np.c_[ + np.linspace(max_loc[0], min_loc[0], 5)[1:-1], + np.ones(3) * min_loc[1], + ], + ] + ) + loop = (loop - np.mean(loop, axis=0)) * 1.5 + np.mean(loop, axis=0) + loop = np.c_[ + loop, gaussian_topo_drape(loop[:, 0], loop[:, 1], flatten=flatten) + ] + loops += [loop + np.asarray(center)] + loop_cells += [np.c_[np.arange(15) + count, np.arange(15) + count + 1]] + loop_cells += [np.c_[count + 15, count]] + count += 16 + + transmitters = LargeLoopGroundTEMTransmitters.create( + geoh5, + vertices=np.vstack(loops), + cells=np.vstack(loop_cells), + ) + transmitters.tx_id_property = transmitters.parts + 1 + survey = LargeLoopGroundTEMReceivers.create(geoh5, vertices=np.vstack(vertices)) + survey.transmitters = transmitters + survey.tx_id_property = np.hstack(loop_id) + # survey.parts = np.repeat(np.arange(n_lines), n_electrodes) + + survey.channels = channels + + survey.waveform = waveform + survey.timing_mark = 2.0 + survey.unit = "Milliseconds (ms)" + + return survey \ No newline at end of file diff --git a/geoapps_utils/utils/locations.py b/geoapps_utils/utils/locations.py index a8a9d958..da09384f 100644 --- a/geoapps_utils/utils/locations.py +++ b/geoapps_utils/utils/locations.py @@ -13,13 +13,46 @@ from uuid import UUID import numpy as np +from typing import Callable from geoh5py import Workspace from geoh5py.data import Data -from geoh5py.objects import Grid2D, Points +from geoh5py.objects import Grid2D, Points, CellObject from geoh5py.objects.grid_object import GridObject from scipy.interpolate import LinearNDInterpolator from scipy.spatial import Delaunay, cKDTree +def gaussian( + x: np.ndarray, y: np.ndarray, amplitude: float, width: float +) -> np.ndarray: + """ + Gaussian function for 2D data. + + :param x: X-coordinates. + :param y: Y-coordinates. + :param amplitude: Amplitude of the Gaussian. + :param width: Width parameter of the Gaussian. + """ + + return amplitude * np.exp(-0.5 * ((x / width) ** 2.0 + (y / width) ** 2.0)) + +def mask_large_connections(cell_object: CellObject, distance_threshold: float): + """ + Trim connections in cell based objects. + + :param cell_object: Cell object containing segments with small vertex spacing + along-line, but large spacing between segments. + + :return: Cleaned object without cells exceeding the distance threshold. + """ + + dist = np.linalg.norm( + cell_object.vertices[cell_object.cells[:, 0], :] + - cell_object.vertices[cell_object.cells[:, 1], :], + axis=1, + ) + + return np.where(dist > distance_threshold)[0] + def mask_under_horizon(locations: np.ndarray, horizon: np.ndarray) -> np.ndarray: """ From 3777833d2b656d68fe6efcbaf8c83ab6763d1c37 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Thu, 31 Jul 2025 21:33:54 +0000 Subject: [PATCH 2/5] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- geoapps_utils/modelling/plate_in_a_halfspace.py | 9 ++------- geoapps_utils/modelling/surveys/dcip.py | 2 +- .../modelling/surveys/frequency_domain/__init__.py | 2 +- .../modelling/surveys/frequency_domain/fdem.py | 10 +++------- .../modelling/surveys/time_domain/__init__.py | 3 ++- .../modelling/surveys/time_domain/airborne_tdem.py | 9 ++++----- .../modelling/surveys/time_domain/ground_tdem.py | 11 ++++------- geoapps_utils/utils/locations.py | 5 +++-- 8 files changed, 20 insertions(+), 31 deletions(-) diff --git a/geoapps_utils/modelling/plate_in_a_halfspace.py b/geoapps_utils/modelling/plate_in_a_halfspace.py index b53ddaf4..9943db5d 100644 --- a/geoapps_utils/modelling/plate_in_a_halfspace.py +++ b/geoapps_utils/modelling/plate_in_a_halfspace.py @@ -8,22 +8,17 @@ # ' # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -import numpy as np from pathlib import Path -from scipy.spatial import Delaunay +import numpy as np from geoh5py import Workspace from geoh5py.objects import Points, Surface - +from scipy.spatial import Delaunay from geoapps_utils.modelling.plates import PlateModel, make_plate from tests.testing_utils.terrain import gaussian_topo_drape - - - - plate_model_default = PlateModel( strike_length=40.0, dip_length=40.0, diff --git a/geoapps_utils/modelling/surveys/dcip.py b/geoapps_utils/modelling/surveys/dcip.py index ea5b60c7..e297b522 100644 --- a/geoapps_utils/modelling/surveys/dcip.py +++ b/geoapps_utils/modelling/surveys/dcip.py @@ -72,4 +72,4 @@ def generate_dc_survey( potentials.current_electrodes = currents currents.potential_electrodes = potentials - return potentials \ No newline at end of file + return potentials diff --git a/geoapps_utils/modelling/surveys/frequency_domain/__init__.py b/geoapps_utils/modelling/surveys/frequency_domain/__init__.py index c50ec18f..abe19e1f 100644 --- a/geoapps_utils/modelling/surveys/frequency_domain/__init__.py +++ b/geoapps_utils/modelling/surveys/frequency_domain/__init__.py @@ -12,4 +12,4 @@ {"Coaxial data": False, "Frequency": 900.0, "Offset": 7.86}, {"Coaxial data": False, "Frequency": 7200.0, "Offset": 7.86}, {"Coaxial data": False, "Frequency": 56000.0, "Offset": 6.3}, -] \ No newline at end of file +] diff --git a/geoapps_utils/modelling/surveys/frequency_domain/fdem.py b/geoapps_utils/modelling/surveys/frequency_domain/fdem.py index 4c2f4971..0e124f37 100644 --- a/geoapps_utils/modelling/surveys/frequency_domain/fdem.py +++ b/geoapps_utils/modelling/surveys/frequency_domain/fdem.py @@ -11,14 +11,13 @@ import numpy as np from geoh5py import Workspace from geoh5py.objects import AirborneFEMReceivers, AirborneFEMTransmitters + from geoapps_utils.utils.locations import mask_large_connections from . import frequency_config -def generate_fdem_survey( - geoh5: Workspace, - vertices: np.ndarray -): + +def generate_fdem_survey(geoh5: Workspace, vertices: np.ndarray): survey = AirborneFEMReceivers.create(geoh5, vertices=vertices, name="Airborne_rx") survey.metadata["EM Dataset"]["Frequency configurations"] = frequency_config @@ -53,7 +52,4 @@ def generate_fdem_survey( survey.remove_cells(mask_large_connections(survey, 200.0)) transmitters.remove_cells(mask_large_connections(transmitters, 200.0)) - return survey - - diff --git a/geoapps_utils/modelling/surveys/time_domain/__init__.py b/geoapps_utils/modelling/surveys/time_domain/__init__.py index cb618213..3a8c6bb6 100644 --- a/geoapps_utils/modelling/surveys/time_domain/__init__.py +++ b/geoapps_utils/modelling/surveys/time_domain/__init__.py @@ -10,6 +10,7 @@ import numpy as np + channels = np.r_[3e-04, 6e-04, 1.2e-03] * 1e3 waveform = np.c_[ np.r_[ @@ -20,4 +21,4 @@ * 1e3 + 2.0, np.r_[np.linspace(0, 1, 4), np.linspace(0.9, 0.0, 4), np.zeros(4)], -] \ No newline at end of file +] diff --git a/geoapps_utils/modelling/surveys/time_domain/airborne_tdem.py b/geoapps_utils/modelling/surveys/time_domain/airborne_tdem.py index e6236bf6..d127408c 100644 --- a/geoapps_utils/modelling/surveys/time_domain/airborne_tdem.py +++ b/geoapps_utils/modelling/surveys/time_domain/airborne_tdem.py @@ -16,16 +16,15 @@ ) from geoapps_utils.utils.locations import mask_large_connections + from . import channels, waveform + def generate_airborne_tdem_survey( geoh5: Workspace, vertices: np.ndarray, ): - - survey = AirborneTEMReceivers.create( - geoh5, vertices=vertices, name="Airborne_rx" - ) + survey = AirborneTEMReceivers.create(geoh5, vertices=vertices, name="Airborne_rx") transmitters = AirborneTEMTransmitters.create( geoh5, vertices=vertices, name="Airborne_tx" ) @@ -39,4 +38,4 @@ def generate_airborne_tdem_survey( survey.timing_mark = 2.0 survey.unit = "Milliseconds (ms)" - return survey \ No newline at end of file + return survey diff --git a/geoapps_utils/modelling/surveys/time_domain/ground_tdem.py b/geoapps_utils/modelling/surveys/time_domain/ground_tdem.py index d3a3ead8..a4da0e92 100644 --- a/geoapps_utils/modelling/surveys/time_domain/ground_tdem.py +++ b/geoapps_utils/modelling/surveys/time_domain/ground_tdem.py @@ -9,7 +9,6 @@ # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' import numpy as np - from geoh5py import Workspace from geoh5py.objects import ( LargeLoopGroundTEMReceivers, @@ -17,16 +16,16 @@ ) from tests.testing_utils.terrain import gaussian_topo_drape + from . import channels, waveform + def generate_tdem_survey( geoh5: Workspace, vertices: np.ndarray, n_lines: int, flatten: bool = False, ): - - X = vertices[:, 0].reshape((n_lines, -1)) Y = vertices[:, 1].reshape((n_lines, -1)) Z = vertices[:, 2].reshape((n_lines, -1)) @@ -72,9 +71,7 @@ def generate_tdem_survey( ] ) loop = (loop - np.mean(loop, axis=0)) * 1.5 + np.mean(loop, axis=0) - loop = np.c_[ - loop, gaussian_topo_drape(loop[:, 0], loop[:, 1], flatten=flatten) - ] + loop = np.c_[loop, gaussian_topo_drape(loop[:, 0], loop[:, 1], flatten=flatten)] loops += [loop + np.asarray(center)] loop_cells += [np.c_[np.arange(15) + count, np.arange(15) + count + 1]] loop_cells += [np.c_[count + 15, count]] @@ -97,4 +94,4 @@ def generate_tdem_survey( survey.timing_mark = 2.0 survey.unit = "Milliseconds (ms)" - return survey \ No newline at end of file + return survey diff --git a/geoapps_utils/utils/locations.py b/geoapps_utils/utils/locations.py index da09384f..97bf3e6b 100644 --- a/geoapps_utils/utils/locations.py +++ b/geoapps_utils/utils/locations.py @@ -13,14 +13,14 @@ from uuid import UUID import numpy as np -from typing import Callable from geoh5py import Workspace from geoh5py.data import Data -from geoh5py.objects import Grid2D, Points, CellObject +from geoh5py.objects import CellObject, Grid2D, Points from geoh5py.objects.grid_object import GridObject from scipy.interpolate import LinearNDInterpolator from scipy.spatial import Delaunay, cKDTree + def gaussian( x: np.ndarray, y: np.ndarray, amplitude: float, width: float ) -> np.ndarray: @@ -35,6 +35,7 @@ def gaussian( return amplitude * np.exp(-0.5 * ((x / width) ** 2.0 + (y / width) ** 2.0)) + def mask_large_connections(cell_object: CellObject, distance_threshold: float): """ Trim connections in cell based objects. From 2076221fd14e9d3bde2d239c92f1289c227eeafc Mon Sep 17 00:00:00 2001 From: benjamink Date: Fri, 8 Aug 2025 08:42:17 -0700 Subject: [PATCH 3/5] remove code that moved to simpeg-drivers --- geoapps_utils/modelling/halfspace.py | 9 - .../modelling/plate_in_a_halfspace.py | 222 ------------------ geoapps_utils/modelling/surveys/dcip.py | 2 +- .../surveys/frequency_domain/__init__.py | 2 +- .../surveys/frequency_domain/fdem.py | 10 +- .../modelling/surveys/time_domain/__init__.py | 3 +- .../surveys/time_domain/airborne_tdem.py | 9 +- .../surveys/time_domain/ground_tdem.py | 11 +- geoapps_utils/utils/locations.py | 5 +- 9 files changed, 18 insertions(+), 255 deletions(-) delete mode 100644 geoapps_utils/modelling/halfspace.py delete mode 100644 geoapps_utils/modelling/plate_in_a_halfspace.py diff --git a/geoapps_utils/modelling/halfspace.py b/geoapps_utils/modelling/halfspace.py deleted file mode 100644 index 76e3205e..00000000 --- a/geoapps_utils/modelling/halfspace.py +++ /dev/null @@ -1,9 +0,0 @@ -# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2025 Mira Geoscience Ltd. ' -# ' -# This file is part of geoapps-utils package. ' -# ' -# geoapps-utils is distributed under the terms and conditions of the MIT License ' -# (see LICENSE file at the root of this source code package). ' -# ' -# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' diff --git a/geoapps_utils/modelling/plate_in_a_halfspace.py b/geoapps_utils/modelling/plate_in_a_halfspace.py deleted file mode 100644 index b53ddaf4..00000000 --- a/geoapps_utils/modelling/plate_in_a_halfspace.py +++ /dev/null @@ -1,222 +0,0 @@ -# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2025 Mira Geoscience Ltd. ' -# ' -# This file is part of geoapps-utils package. ' -# ' -# geoapps-utils is distributed under the terms and conditions of the MIT License ' -# (see LICENSE file at the root of this source code package). ' -# ' -# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' - -import numpy as np -from pathlib import Path -from scipy.spatial import Delaunay - -from geoh5py import Workspace -from geoh5py.objects import Points, Surface - - -from geoapps_utils.modelling.plates import PlateModel, make_plate -from tests.testing_utils.terrain import gaussian_topo_drape - - - - - - -plate_model_default = PlateModel( - strike_length=40.0, - dip_length=40.0, - width=40.0, - origin=(0.0, 0.0, 10.0), -) - - -def setup_inversion_workspace( - work_dir, - plate_model: PlateModel = plate_model_default, - background=None, - anomaly=None, - cell_size=(5.0, 5.0, 5.0), - center=(0.0, 0.0, 0.0), - n_electrodes=20, - n_lines=5, - refinement=(4, 6), - x_limits=(-100.0, 100.0), - y_limits=(-100.0, 100.0), - padding_distance=100, - drape_height=5.0, - inversion_type="other", - flatten=False, - geoh5=None, -): - """ - Creates a workspace populated with objects to simulate/invert a simple model. - - rot_xyz: Rotation angles in degrees about the x, y, and z axes. - """ - filepath = Path(work_dir) / "inversion_test.ui.geoh5" - if geoh5 is None: - if filepath.is_file(): - filepath.unlink() - geoh5 = Workspace.create(filepath) - # Topography - xx, yy = np.meshgrid(np.linspace(-200.0, 200.0, 50), np.linspace(-200.0, 200.0, 50)) - - zz = gaussian_topo_drape(xx, yy, flatten=flatten) - topo = np.c_[ - utils.mkvc(xx) + center[0], - utils.mkvc(yy) + center[1], - utils.mkvc(zz) + center[2], - ] - triang = Delaunay(topo[:, :2]) - topography = Surface.create( - geoh5, - vertices=topo, - cells=triang.simplices, # pylint: disable=E1101 - name="topography", - ) - - # Observation points - n_electrodes = ( - 4 - if ("polarization" in inversion_type or "current" in inversion_type) - & (n_electrodes < 4) - else n_electrodes - ) - xr = np.linspace(x_limits[0], x_limits[1], int(n_electrodes)) - yr = np.linspace(y_limits[0], y_limits[1], int(n_lines)) - X, Y = np.meshgrid(xr, yr) - Z = gaussian_topo_drape(X, Y, flatten=flatten) + drape_height - - vertices = np.c_[ - utils.mkvc(X.T) + center[0], - utils.mkvc(Y.T) + center[1], - utils.mkvc(Z.T) + center[2], - ] - - if "polarization" in inversion_type or "current" in inversion_type: - survey = generate_dc_survey(geoh5, X, Y, Z) - - elif "magnetotellurics" in inversion_type: - survey = MTReceivers.create( - geoh5, - vertices=vertices, - name="survey", - components=[ - "Zxx (real)", - "Zxx (imag)", - "Zxy (real)", - "Zxy (imag)", - "Zyx (real)", - "Zyx (imag)", - "Zyy (real)", - "Zyy (imag)", - ], - channels=[10.0, 100.0, 1000.0], - ) - - elif "tipper" in inversion_type: - survey = TipperReceivers.create( - geoh5, - vertices=vertices, - name="survey", - components=[ - "Txz (real)", - "Txz (imag)", - "Tyz (real)", - "Tyz (imag)", - ], - ) - survey.base_stations = TipperBaseStations.create( - geoh5, vertices=np.c_[vertices[0, :]].T - ) - survey.channels = [10.0, 100.0, 1000.0] - dist = np.linalg.norm( - survey.vertices[survey.cells[:, 0], :] - - survey.vertices[survey.cells[:, 1], :], - axis=1, - ) - # survey.cells = survey.cells[dist < 100.0, :] - survey.remove_cells(np.where(dist > 200)[0]) - - elif inversion_type in ["fdem", "fem", "fdem 1d"]: - survey = generate_fdem_survey(geoh5, vertices) - - elif "tdem" in inversion_type: - survey = generate_tdem_survey( - geoh5, - vertices, - n_lines, - flatten=flatten, - airborne="airborne" in inversion_type, - ) - - else: - survey = Points.create( - geoh5, - vertices=vertices, - name="survey", - ) - - # Create a mesh - if "2d" in inversion_type: - lines = survey.get_entity("line_ids")[0].values - entity, mesh, _ = get_drape_model( # pylint: disable=unbalanced-tuple-unpacking - geoh5, - "Models", - survey.vertices[np.unique(survey.cells[lines == 101, :]), :], - [cell_size[0], cell_size[2]], - 100.0, - [padding_distance] * 2 + [padding_distance] * 2, - 1.1, - parent=None, - return_colocated_mesh=True, - return_sorting=True, - ) - active = active_from_xyz(entity, topography.vertices, grid_reference="top") - - else: - padDist = np.ones((3, 2)) * padding_distance - mesh = mesh_builder_xyz( - vertices - np.r_[cell_size] / 2.0, - cell_size, - depth_core=100.0, - padding_distance=padDist, - mesh_type="TREE", - tree_diagonal_balance=False, - ) - mesh = OctreeDriver.refine_tree_from_surface( - mesh, - topography, - levels=refinement, - finalize=False, - ) - - if inversion_type in ["fdem", "airborne tdem"]: - mesh = OctreeDriver.refine_tree_from_points( - mesh, - vertices, - levels=[2], - finalize=False, - ) - - mesh.finalize() - entity = treemesh_2_octree(geoh5, mesh, name="mesh") - active = active_from_xyz(entity, topography.vertices, grid_reference="top") - - # Create the Model - cell_centers = entity.centroids.copy() - - model = make_plate( - cell_centers, plate_model, background=background, anomaly=anomaly - ) - - if "1d" in inversion_type: - model = background * np.ones(mesh.nC) - model[(mesh.cell_centers[:, 2] < 0) & (mesh.cell_centers[:, 2] > -20)] = anomaly - - model[~active] = np.nan - model = entity.add_data({"model": {"values": model}}) - geoh5.close() - return geoh5, entity, model, survey, topography diff --git a/geoapps_utils/modelling/surveys/dcip.py b/geoapps_utils/modelling/surveys/dcip.py index ea5b60c7..e297b522 100644 --- a/geoapps_utils/modelling/surveys/dcip.py +++ b/geoapps_utils/modelling/surveys/dcip.py @@ -72,4 +72,4 @@ def generate_dc_survey( potentials.current_electrodes = currents currents.potential_electrodes = potentials - return potentials \ No newline at end of file + return potentials diff --git a/geoapps_utils/modelling/surveys/frequency_domain/__init__.py b/geoapps_utils/modelling/surveys/frequency_domain/__init__.py index c50ec18f..abe19e1f 100644 --- a/geoapps_utils/modelling/surveys/frequency_domain/__init__.py +++ b/geoapps_utils/modelling/surveys/frequency_domain/__init__.py @@ -12,4 +12,4 @@ {"Coaxial data": False, "Frequency": 900.0, "Offset": 7.86}, {"Coaxial data": False, "Frequency": 7200.0, "Offset": 7.86}, {"Coaxial data": False, "Frequency": 56000.0, "Offset": 6.3}, -] \ No newline at end of file +] diff --git a/geoapps_utils/modelling/surveys/frequency_domain/fdem.py b/geoapps_utils/modelling/surveys/frequency_domain/fdem.py index 4c2f4971..0e124f37 100644 --- a/geoapps_utils/modelling/surveys/frequency_domain/fdem.py +++ b/geoapps_utils/modelling/surveys/frequency_domain/fdem.py @@ -11,14 +11,13 @@ import numpy as np from geoh5py import Workspace from geoh5py.objects import AirborneFEMReceivers, AirborneFEMTransmitters + from geoapps_utils.utils.locations import mask_large_connections from . import frequency_config -def generate_fdem_survey( - geoh5: Workspace, - vertices: np.ndarray -): + +def generate_fdem_survey(geoh5: Workspace, vertices: np.ndarray): survey = AirborneFEMReceivers.create(geoh5, vertices=vertices, name="Airborne_rx") survey.metadata["EM Dataset"]["Frequency configurations"] = frequency_config @@ -53,7 +52,4 @@ def generate_fdem_survey( survey.remove_cells(mask_large_connections(survey, 200.0)) transmitters.remove_cells(mask_large_connections(transmitters, 200.0)) - return survey - - diff --git a/geoapps_utils/modelling/surveys/time_domain/__init__.py b/geoapps_utils/modelling/surveys/time_domain/__init__.py index cb618213..3a8c6bb6 100644 --- a/geoapps_utils/modelling/surveys/time_domain/__init__.py +++ b/geoapps_utils/modelling/surveys/time_domain/__init__.py @@ -10,6 +10,7 @@ import numpy as np + channels = np.r_[3e-04, 6e-04, 1.2e-03] * 1e3 waveform = np.c_[ np.r_[ @@ -20,4 +21,4 @@ * 1e3 + 2.0, np.r_[np.linspace(0, 1, 4), np.linspace(0.9, 0.0, 4), np.zeros(4)], -] \ No newline at end of file +] diff --git a/geoapps_utils/modelling/surveys/time_domain/airborne_tdem.py b/geoapps_utils/modelling/surveys/time_domain/airborne_tdem.py index e6236bf6..d127408c 100644 --- a/geoapps_utils/modelling/surveys/time_domain/airborne_tdem.py +++ b/geoapps_utils/modelling/surveys/time_domain/airborne_tdem.py @@ -16,16 +16,15 @@ ) from geoapps_utils.utils.locations import mask_large_connections + from . import channels, waveform + def generate_airborne_tdem_survey( geoh5: Workspace, vertices: np.ndarray, ): - - survey = AirborneTEMReceivers.create( - geoh5, vertices=vertices, name="Airborne_rx" - ) + survey = AirborneTEMReceivers.create(geoh5, vertices=vertices, name="Airborne_rx") transmitters = AirborneTEMTransmitters.create( geoh5, vertices=vertices, name="Airborne_tx" ) @@ -39,4 +38,4 @@ def generate_airborne_tdem_survey( survey.timing_mark = 2.0 survey.unit = "Milliseconds (ms)" - return survey \ No newline at end of file + return survey diff --git a/geoapps_utils/modelling/surveys/time_domain/ground_tdem.py b/geoapps_utils/modelling/surveys/time_domain/ground_tdem.py index d3a3ead8..a4da0e92 100644 --- a/geoapps_utils/modelling/surveys/time_domain/ground_tdem.py +++ b/geoapps_utils/modelling/surveys/time_domain/ground_tdem.py @@ -9,7 +9,6 @@ # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' import numpy as np - from geoh5py import Workspace from geoh5py.objects import ( LargeLoopGroundTEMReceivers, @@ -17,16 +16,16 @@ ) from tests.testing_utils.terrain import gaussian_topo_drape + from . import channels, waveform + def generate_tdem_survey( geoh5: Workspace, vertices: np.ndarray, n_lines: int, flatten: bool = False, ): - - X = vertices[:, 0].reshape((n_lines, -1)) Y = vertices[:, 1].reshape((n_lines, -1)) Z = vertices[:, 2].reshape((n_lines, -1)) @@ -72,9 +71,7 @@ def generate_tdem_survey( ] ) loop = (loop - np.mean(loop, axis=0)) * 1.5 + np.mean(loop, axis=0) - loop = np.c_[ - loop, gaussian_topo_drape(loop[:, 0], loop[:, 1], flatten=flatten) - ] + loop = np.c_[loop, gaussian_topo_drape(loop[:, 0], loop[:, 1], flatten=flatten)] loops += [loop + np.asarray(center)] loop_cells += [np.c_[np.arange(15) + count, np.arange(15) + count + 1]] loop_cells += [np.c_[count + 15, count]] @@ -97,4 +94,4 @@ def generate_tdem_survey( survey.timing_mark = 2.0 survey.unit = "Milliseconds (ms)" - return survey \ No newline at end of file + return survey diff --git a/geoapps_utils/utils/locations.py b/geoapps_utils/utils/locations.py index da09384f..97bf3e6b 100644 --- a/geoapps_utils/utils/locations.py +++ b/geoapps_utils/utils/locations.py @@ -13,14 +13,14 @@ from uuid import UUID import numpy as np -from typing import Callable from geoh5py import Workspace from geoh5py.data import Data -from geoh5py.objects import Grid2D, Points, CellObject +from geoh5py.objects import CellObject, Grid2D, Points from geoh5py.objects.grid_object import GridObject from scipy.interpolate import LinearNDInterpolator from scipy.spatial import Delaunay, cKDTree + def gaussian( x: np.ndarray, y: np.ndarray, amplitude: float, width: float ) -> np.ndarray: @@ -35,6 +35,7 @@ def gaussian( return amplitude * np.exp(-0.5 * ((x / width) ** 2.0 + (y / width) ** 2.0)) + def mask_large_connections(cell_object: CellObject, distance_threshold: float): """ Trim connections in cell based objects. From 9dbd979523340b2b0808a7af6e51568f16cc2e6e Mon Sep 17 00:00:00 2001 From: benjamink Date: Mon, 11 Aug 2025 09:52:49 -0700 Subject: [PATCH 4/5] remove surveys module --- geoapps_utils/modelling/surveys/__init__.py | 9 -- geoapps_utils/modelling/surveys/dcip.py | 75 -------------- .../surveys/frequency_domain/__init__.py | 15 --- .../surveys/frequency_domain/fdem.py | 55 ----------- .../modelling/surveys/time_domain/__init__.py | 24 ----- .../surveys/time_domain/airborne_tdem.py | 41 -------- .../surveys/time_domain/ground_tdem.py | 97 ------------------- 7 files changed, 316 deletions(-) delete mode 100644 geoapps_utils/modelling/surveys/__init__.py delete mode 100644 geoapps_utils/modelling/surveys/dcip.py delete mode 100644 geoapps_utils/modelling/surveys/frequency_domain/__init__.py delete mode 100644 geoapps_utils/modelling/surveys/frequency_domain/fdem.py delete mode 100644 geoapps_utils/modelling/surveys/time_domain/__init__.py delete mode 100644 geoapps_utils/modelling/surveys/time_domain/airborne_tdem.py delete mode 100644 geoapps_utils/modelling/surveys/time_domain/ground_tdem.py diff --git a/geoapps_utils/modelling/surveys/__init__.py b/geoapps_utils/modelling/surveys/__init__.py deleted file mode 100644 index 76e3205e..00000000 --- a/geoapps_utils/modelling/surveys/__init__.py +++ /dev/null @@ -1,9 +0,0 @@ -# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2025 Mira Geoscience Ltd. ' -# ' -# This file is part of geoapps-utils package. ' -# ' -# geoapps-utils is distributed under the terms and conditions of the MIT License ' -# (see LICENSE file at the root of this source code package). ' -# ' -# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' diff --git a/geoapps_utils/modelling/surveys/dcip.py b/geoapps_utils/modelling/surveys/dcip.py deleted file mode 100644 index e297b522..00000000 --- a/geoapps_utils/modelling/surveys/dcip.py +++ /dev/null @@ -1,75 +0,0 @@ -# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2025 Mira Geoscience Ltd. ' -# ' -# This file is part of geoapps-utils package. ' -# ' -# geoapps-utils is distributed under the terms and conditions of the MIT License ' -# (see LICENSE file at the root of this source code package). ' -# ' -# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' - -import numpy as np -from geoh5py import Workspace -from geoh5py.objects import CurrentElectrode, PotentialElectrode - - -def generate_dc_survey( - workspace: Workspace, - x_loc: np.ndarray, - y_loc: np.ndarray, - z_loc: np.ndarray | None = None, -) -> PotentialElectrode: - """ - Utility function to generate a DC survey. - """ - # Create sources along line - if z_loc is None: - z_loc = np.zeros_like(x_loc) - - vertices = np.c_[x_loc.ravel(), y_loc.ravel(), z_loc.ravel()] - parts = np.kron(np.arange(x_loc.shape[0]), np.ones(x_loc.shape[1])).astype("int") - currents = CurrentElectrode.create(workspace, vertices=vertices, parts=parts) - currents.add_default_ab_cell_id() - n_dipoles = 9 - dipoles = [] - current_id = [] - - for val in currents.ab_cell_id.values: - cell_id = int(currents.ab_map[val]) - 1 - - for dipole in range(n_dipoles): - dipole_ids = currents.cells[cell_id, :] + 2 + dipole - - if ( - any(dipole_ids > (currents.n_vertices - 1)) - or len( - np.unique(parts[np.r_[currents.cells[cell_id, 0], dipole_ids[1]]]) - ) - > 1 - ): - continue - - dipoles += [dipole_ids] - current_id += [val] - - potentials = PotentialElectrode.create( - workspace, vertices=vertices, cells=np.vstack(dipoles).astype("uint32") - ) - line_id = potentials.vertices[potentials.cells[:, 0], 1] - line_id = (line_id - np.min(line_id) + 1).astype(np.int32) - line_reference = {0: "Unknown"} - line_reference.update({k: str(k) for k in np.unique(line_id)}) - potentials.add_data( - { - "line_ids": { - "values": line_id, - "type": "REFERENCED", - "value_map": line_reference, - } - } - ) - potentials.ab_cell_id = np.hstack(current_id).astype("int32") - potentials.current_electrodes = currents - currents.potential_electrodes = potentials - - return potentials diff --git a/geoapps_utils/modelling/surveys/frequency_domain/__init__.py b/geoapps_utils/modelling/surveys/frequency_domain/__init__.py deleted file mode 100644 index abe19e1f..00000000 --- a/geoapps_utils/modelling/surveys/frequency_domain/__init__.py +++ /dev/null @@ -1,15 +0,0 @@ -# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2025 Mira Geoscience Ltd. ' -# ' -# This file is part of geoapps-utils package. ' -# ' -# geoapps-utils is distributed under the terms and conditions of the MIT License ' -# (see LICENSE file at the root of this source code package). ' -# ' -# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' - -frequency_config = [ - {"Coaxial data": False, "Frequency": 900.0, "Offset": 7.86}, - {"Coaxial data": False, "Frequency": 7200.0, "Offset": 7.86}, - {"Coaxial data": False, "Frequency": 56000.0, "Offset": 6.3}, -] diff --git a/geoapps_utils/modelling/surveys/frequency_domain/fdem.py b/geoapps_utils/modelling/surveys/frequency_domain/fdem.py deleted file mode 100644 index 0e124f37..00000000 --- a/geoapps_utils/modelling/surveys/frequency_domain/fdem.py +++ /dev/null @@ -1,55 +0,0 @@ -# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2025 Mira Geoscience Ltd. ' -# ' -# This file is part of geoapps-utils package. ' -# ' -# geoapps-utils is distributed under the terms and conditions of the MIT License ' -# (see LICENSE file at the root of this source code package). ' -# ' -# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' - -import numpy as np -from geoh5py import Workspace -from geoh5py.objects import AirborneFEMReceivers, AirborneFEMTransmitters - -from geoapps_utils.utils.locations import mask_large_connections - -from . import frequency_config - - -def generate_fdem_survey(geoh5: Workspace, vertices: np.ndarray): - survey = AirborneFEMReceivers.create(geoh5, vertices=vertices, name="Airborne_rx") - - survey.metadata["EM Dataset"]["Frequency configurations"] = frequency_config - - tx_locs = [] - freqs = [] - for config in frequency_config: - tx_vertices = vertices.copy() - tx_vertices[:, 0] -= config["Offset"] - tx_locs.append(tx_vertices) - freqs.append([[config["Frequency"]] * len(vertices)]) - tx_locs = np.vstack(tx_locs) - freqs = np.hstack(freqs).flatten() - - transmitters = AirborneFEMTransmitters.create( - geoh5, vertices=tx_locs, name="Airborne_tx" - ) - survey.transmitters = transmitters - survey.channels = [k["Frequency"] for k in frequency_config] - - transmitters.add_data( - { - "Tx frequency": { - "values": freqs, - "association": "VERTEX", - "primitive_type": "REFERENCED", - "value_map": {k: str(k) for k in freqs}, - } - } - ) - - survey.remove_cells(mask_large_connections(survey, 200.0)) - transmitters.remove_cells(mask_large_connections(transmitters, 200.0)) - - return survey diff --git a/geoapps_utils/modelling/surveys/time_domain/__init__.py b/geoapps_utils/modelling/surveys/time_domain/__init__.py deleted file mode 100644 index 3a8c6bb6..00000000 --- a/geoapps_utils/modelling/surveys/time_domain/__init__.py +++ /dev/null @@ -1,24 +0,0 @@ -# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2025 Mira Geoscience Ltd. ' -# ' -# This file is part of geoapps-utils package. ' -# ' -# geoapps-utils is distributed under the terms and conditions of the MIT License ' -# (see LICENSE file at the root of this source code package). ' -# ' -# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' - -import numpy as np - - -channels = np.r_[3e-04, 6e-04, 1.2e-03] * 1e3 -waveform = np.c_[ - np.r_[ - np.arange(-0.002, -0.0001, 5e-4), - np.arange(-0.0004, 0.0, 1e-4), - np.arange(0.0, 0.002, 5e-4), - ] - * 1e3 - + 2.0, - np.r_[np.linspace(0, 1, 4), np.linspace(0.9, 0.0, 4), np.zeros(4)], -] diff --git a/geoapps_utils/modelling/surveys/time_domain/airborne_tdem.py b/geoapps_utils/modelling/surveys/time_domain/airborne_tdem.py deleted file mode 100644 index d127408c..00000000 --- a/geoapps_utils/modelling/surveys/time_domain/airborne_tdem.py +++ /dev/null @@ -1,41 +0,0 @@ -# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2025 Mira Geoscience Ltd. ' -# ' -# This file is part of geoapps-utils package. ' -# ' -# geoapps-utils is distributed under the terms and conditions of the MIT License ' -# (see LICENSE file at the root of this source code package). ' -# ' -# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' - -import numpy as np -from geoh5py import Workspace -from geoh5py.objects import ( - AirborneTEMReceivers, - AirborneTEMTransmitters, -) - -from geoapps_utils.utils.locations import mask_large_connections - -from . import channels, waveform - - -def generate_airborne_tdem_survey( - geoh5: Workspace, - vertices: np.ndarray, -): - survey = AirborneTEMReceivers.create(geoh5, vertices=vertices, name="Airborne_rx") - transmitters = AirborneTEMTransmitters.create( - geoh5, vertices=vertices, name="Airborne_tx" - ) - mask = mask_large_connections(survey, 200.0) - survey.remove_cells(mask) - transmitters.remove_cells(mask) - - survey.transmitters = transmitters - survey.channels = channels - survey.waveform = waveform - survey.timing_mark = 2.0 - survey.unit = "Milliseconds (ms)" - - return survey diff --git a/geoapps_utils/modelling/surveys/time_domain/ground_tdem.py b/geoapps_utils/modelling/surveys/time_domain/ground_tdem.py deleted file mode 100644 index a4da0e92..00000000 --- a/geoapps_utils/modelling/surveys/time_domain/ground_tdem.py +++ /dev/null @@ -1,97 +0,0 @@ -# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2025 Mira Geoscience Ltd. ' -# ' -# This file is part of geoapps-utils package. ' -# ' -# geoapps-utils is distributed under the terms and conditions of the MIT License ' -# (see LICENSE file at the root of this source code package). ' -# ' -# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' - -import numpy as np -from geoh5py import Workspace -from geoh5py.objects import ( - LargeLoopGroundTEMReceivers, - LargeLoopGroundTEMTransmitters, -) - -from tests.testing_utils.terrain import gaussian_topo_drape - -from . import channels, waveform - - -def generate_tdem_survey( - geoh5: Workspace, - vertices: np.ndarray, - n_lines: int, - flatten: bool = False, -): - X = vertices[:, 0].reshape((n_lines, -1)) - Y = vertices[:, 1].reshape((n_lines, -1)) - Z = vertices[:, 2].reshape((n_lines, -1)) - center = np.mean(vertices, axis=0) - arrays = [ - np.c_[ - X[: int(n_lines / 2), :].flatten(), - Y[: int(n_lines / 2), :].flatten(), - Z[: int(n_lines / 2), :].flatten(), - ], - np.c_[ - X[int(n_lines / 2) :, :].flatten(), - Y[int(n_lines / 2) :, :].flatten(), - Z[int(n_lines / 2) :, :].flatten(), - ], - ] - loops = [] - loop_cells = [] - loop_id = [] - count = 0 - for ind, array in enumerate(arrays): - loop_id += [np.ones(array.shape[0]) * (ind + 1)] - min_loc = np.min(array, axis=0) - max_loc = np.max(array, axis=0) - loop = np.vstack( - [ - np.c_[ - np.ones(5) * min_loc[0], - np.linspace(min_loc[1], max_loc[1], 5), - ], - np.c_[ - np.linspace(min_loc[0], max_loc[0], 5)[1:], - np.ones(4) * max_loc[1], - ], - np.c_[ - np.ones(4) * max_loc[0], - np.linspace(max_loc[1], min_loc[1], 5)[1:], - ], - np.c_[ - np.linspace(max_loc[0], min_loc[0], 5)[1:-1], - np.ones(3) * min_loc[1], - ], - ] - ) - loop = (loop - np.mean(loop, axis=0)) * 1.5 + np.mean(loop, axis=0) - loop = np.c_[loop, gaussian_topo_drape(loop[:, 0], loop[:, 1], flatten=flatten)] - loops += [loop + np.asarray(center)] - loop_cells += [np.c_[np.arange(15) + count, np.arange(15) + count + 1]] - loop_cells += [np.c_[count + 15, count]] - count += 16 - - transmitters = LargeLoopGroundTEMTransmitters.create( - geoh5, - vertices=np.vstack(loops), - cells=np.vstack(loop_cells), - ) - transmitters.tx_id_property = transmitters.parts + 1 - survey = LargeLoopGroundTEMReceivers.create(geoh5, vertices=np.vstack(vertices)) - survey.transmitters = transmitters - survey.tx_id_property = np.hstack(loop_id) - # survey.parts = np.repeat(np.arange(n_lines), n_electrodes) - - survey.channels = channels - - survey.waveform = waveform - survey.timing_mark = 2.0 - survey.unit = "Milliseconds (ms)" - - return survey From e62f330f7a702964c7583b1ad882ba9fe0a70376 Mon Sep 17 00:00:00 2001 From: benjamink Date: Mon, 11 Aug 2025 10:55:13 -0700 Subject: [PATCH 5/5] Add unit tests --- tests/locations_test.py | 24 +++++++++++++++++++++++- 1 file changed, 23 insertions(+), 1 deletion(-) diff --git a/tests/locations_test.py b/tests/locations_test.py index 1309bb25..9277e7a1 100644 --- a/tests/locations_test.py +++ b/tests/locations_test.py @@ -12,17 +12,39 @@ import numpy as np from geoh5py import Workspace -from geoh5py.objects import Grid2D, Points +from geoh5py.objects import Curve, Grid2D, Points from geoapps_utils.utils.locations import ( + gaussian, get_locations, get_overlapping_limits, map_indices_to_coordinates, + mask_large_connections, mask_under_horizon, ) from geoapps_utils.utils.transformations import rotate_points, z_rotation_matrix +def test_gaussian(): + x = np.linspace(-10, 10, 100) + y = np.linspace(-10, 10, 100) + x_grid, y_grid = np.meshgrid(x, y) + z_grid = gaussian(x_grid, y_grid, 10, 5) + assert np.isclose(z_grid.max(), 10, rtol=1e-3) + + +def test_mask_large_connections(tmp_path): + with Workspace(tmp_path / "test.geoh5") as ws: + x = np.linspace(0, 100, 11) + y = np.linspace(0, 300, 4) + x_grid, y_grid = np.meshgrid(x, y) + z_grid = np.zeros_like(x_grid) + vertices = np.column_stack([x_grid.ravel(), y_grid.ravel(), z_grid.ravel()]) + crv = Curve.create(ws, name="test_curve", vertices=vertices) + mask = mask_large_connections(crv, distance_threshold=50.0) + assert len(mask) == 3 + + def test_rotate_points(): points = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1], [2, 3, 4]]) validation = rotate_points(