GEOPY-2770: Consolidate code for plate model definition between different repos

simpeg_drivers-assets/uijson/plate_simulation.ui.json

@@ -33,7 +33,7 @@
         "tooltip": "Value of the basement resisitivity (ohm-m), density (g/cc) or susceptibility (SI)",
         "enabled": true
     },
-    "overburden": {
+    "overburden_property": {
         "main": true,
         "group": "Overburden",
         "label": "Physical property (SI)",
@@ -65,7 +65,7 @@
         "enabled": true,
         "tooltip": "Spacing between plates"
     },
-    "plate": {
+    "plate_property": {
         "main": true,
         "group": "Plate",
         "label": "Physical property (SI)",

simpeg_drivers/plate_simulation/driver.py

@@ -131,24 +131,29 @@ def plates(self) -> list[Plate]:
         """Generate sequence of plates."""
         if self._plates is None:
             offset = (
-                self.params.model.overburden_model.thickness
-                if self.params.model.plate_model.reference_surface == "overburden"
+                self.params.model.overburden.thickness
+                if self.params.model.plate.reference_surface == "overburden"
                 else 0.0
             )
-            center = self.params.model.plate_model.center(
+            center = self.params.model.plate.center(
                 self.survey,
                 self.topography,
                 depth_offset=-1 * offset,
             )
             plate = Plate(
-                self.params.model.plate_model,
-                center,
+                self.params.model.plate.model_copy(
+                    update={
+                        "easting": center[0],
+                        "northing": center[1],
+                        "elevation": center[2],
+                    }
+                ),
             )
             self._plates = self.replicate(
                 plate,
-                self.params.model.plate_model.number,
-                self.params.model.plate_model.spacing,
-                self.params.model.plate_model.dip_direction,
+                self.params.model.plate.number,
+                self.params.model.plate.spacing,
+                self.params.model.plate.geometry.direction,
             )
         return self._plates
 
@@ -198,12 +203,12 @@ def make_model(self) -> FloatData:
 
         overburden = Overburden(
             topography=self.simulation_parameters.active_cells.topography_object,
-            thickness=self.params.model.overburden_model.thickness,
-            value=self.params.model.overburden_model.overburden,
+            thickness=self.params.model.overburden.thickness,
+            value=self.params.model.overburden.overburden_property,
         )
 
         dikes = DikeSwarm(
-            [Anomaly(plate, plate.params.plate) for plate in self.plates],
+            [Anomaly(plate, plate.params.plate_property) for plate in self.plates],
             name="plates",
         )
 
@@ -272,10 +277,24 @@ def replicate(
 
         plates = []
         for i in range(number):
-            center = np.r_[plate.center] + azimuth_to_unit_vector(azimuth) * offsets[i]
-            new = Plate(plate.params.model_copy(), center)
-            new.params.name = f"{plate.params.name} offset {i + 1}"
-            plates.append(new)
+            center = (
+                np.r_[plate.params.geometry.origin]
+                + azimuth_to_unit_vector(azimuth) * offsets[i]
+            )
+            new_geometry = plate.params.geometry.model_copy(
+                update={
+                    "easting": center[0],
+                    "northing": center[1],
+                    "elevation": center[2],
+                }
+            )
+            new_plate = Plate(
+                plate.params.model_copy(update={"geometry": new_geometry})
+            )
+
+            new_plate.params.name = f"{plate.params.name} offset {i + 1}"
+            plates.append(new_plate)
+
         return plates
 
 

simpeg_drivers/plate_simulation/match/driver.py

@@ -192,7 +192,7 @@ def _create_plate_from_parameters(
         """
         center = self.params.survey.vertices[index_center]
         center[2] = (
-            self._drape_heights[index_center] - model_options.overburden_model.thickness
+            self._drape_heights[index_center] - model_options.overburden.thickness
         )
         indices = self.params.survey.get_segment_indices(
             index_center, self.params.max_distance
@@ -208,10 +208,10 @@ def _create_plate_from_parameters(
                     "y": center[1],
                     "z": center[2],
                 },
-                "width": model_options.plate_model.dip_length,
-                "thickness": model_options.plate_model.width,
-                "length": model_options.plate_model.strike_length,
-                "dip": model_options.plate_model.dip,
+                "width": model_options.plate.geometry.dip_length,
+                "thickness": model_options.plate.geometry.width,
+                "length": model_options.plate.geometry.strike_length,
+                "dip": model_options.plate.geometry.dip,
                 "dip_direction": (azimuth + strike_angle) % 360,
             }
         )

simpeg_drivers/plate_simulation/models/options.py

@@ -11,6 +11,7 @@
 from typing import TypeVar
 
 import numpy as np
+from geoapps_utils.modelling.plates import PlateModel
 from geoh5py.objects import Points
 from pydantic import (
     BaseModel,
@@ -29,20 +30,12 @@ class PlateOptions(BaseModel):
     Parameters describing an anomalous plate.
 
     :param plate: Value given to the plate(s).
-    :param width: V-size of the plate.
-    :param strike_length: U-size of the plate.
-    :param dip_length: W-size of the plate.
-    :param dip: Orientation of the v-axis in degree from horizontal.
-    :param dip_direction: Orientation of the u axis in degree from north.
+    :param geometry: Parameters describing the plate geometry.
     :param reference: Point of rotation to be 'center' or 'top'.
     :param number: Number of offset plates to be created.
     :param spacing: Spacing between plates.
     :param relative_locations: If True locations are relative to survey in xy and
         mean topography in z.
-    :param easting: Easting offset relative to survey.
-    :param northing: Northing offset relative to survey.
-    :param elevation: plate(s) elevation.  May be true elevation or relative to
-        overburden or topography.
     :param reference_surface: Switches between using topography and overburden as
         elevation reference of the plate.
     :param reference_type: Type of reference for plate elevation.  Can be 'mean'
@@ -53,18 +46,11 @@ class PlateOptions(BaseModel):
     model_config = ConfigDict(arbitrary_types_allowed=True)
 
     name: str = "Plate"
-    plate: float
-    width: float
-    strike_length: float
-    dip_length: float
-    dip: float = 90.0
-    dip_direction: float = 90.0
+    plate_property: float
+    geometry: PlateModel
     number: int = 1
     spacing: float = 0.0
     relative_locations: bool = False
-    easting: float = 0.0
-    northing: float = 0.0
-    elevation: float
     reference_surface: str = "topography"
     reference_type: str = "mean"
 
@@ -82,7 +68,7 @@ def single_plate(self):
     @property
     def halfplate(self):
         """Compute half the z-projection length of the plate."""
-        return 0.5 * self.dip_length * np.sin(np.deg2rad(self.dip))
+        return 0.5 * self.geometry.dip_length * np.sin(np.deg2rad(self.geometry.dip))
 
     def center(
         self,
@@ -104,11 +90,11 @@ def _get_xy(self, survey: Points) -> tuple[float, float]:
 
         if self.relative_locations:
             return (
-                survey.vertices[:, 0].mean() + self.easting,
-                survey.vertices[:, 1].mean() + self.northing,
+                survey.vertices[:, 0].mean() + self.geometry.origin[0],
+                survey.vertices[:, 1].mean() + self.geometry.origin[1],
             )
 
-        return self.easting, self.northing
+        return self.geometry.origin[0], self.geometry.origin[1]
 
     def _get_z(self, surface: Points, offset: float = 0.0) -> float:
         """
@@ -122,9 +108,9 @@ def _get_z(self, surface: Points, offset: float = 0.0) -> float:
             raise ValueError("Topography object has no vertices.")
         if self.relative_locations:
             z = getattr(surface.vertices[:, 2], self.reference_type)()
-            z += offset + self.elevation - self.halfplate
+            z += offset + self.geometry.elevation - self.halfplate
         else:
-            z = self.elevation
+            z = self.geometry.elevation
 
         return z
 
@@ -138,7 +124,7 @@ class OverburdenOptions(BaseModel):
     """
 
     thickness: float
-    overburden: float
+    overburden_property: float
 
 
 class ModelOptions(BaseModel):
@@ -153,5 +139,5 @@ class ModelOptions(BaseModel):
     model_config = ConfigDict(arbitrary_types_allowed=True)
 
     background: float
-    overburden_model: OverburdenOptions
-    plate_model: PlateOptions
+    overburden: OverburdenOptions
+    plate: PlateOptions

simpeg_drivers/plate_simulation/models/parametric.py

@@ -68,15 +68,9 @@ class Plate(Parametric):
     def __init__(
         self,
         params: PlateOptions,
-        center: tuple[float, float, float] = (
-            0.0,
-            0.0,
-            0.0,
-        ),
         workspace: Workspace | None = None,
     ):
         self.params = params
-        self.center = center
         self._workspace = workspace
         super().__init__(self._create_surface())
 
@@ -128,12 +122,16 @@ def triangles(self) -> np.ndarray:
     def vertices(self) -> np.ndarray:
         """Vertices for triangulation of a rectangular prism in 3D space."""
 
-        u_1 = self.center[0] - (self.params.strike_length / 2.0)
-        u_2 = self.center[0] + (self.params.strike_length / 2.0)
-        v_1 = self.center[1] - (self.params.dip_length / 2.0)
-        v_2 = self.center[1] + (self.params.dip_length / 2.0)
-        w_1 = self.center[2] - (self.params.width / 2.0)
-        w_2 = self.center[2] + (self.params.width / 2.0)
+        u_1 = self.params.geometry.origin[0] - (
+            self.params.geometry.strike_length / 2.0
+        )
+        u_2 = self.params.geometry.origin[0] + (
+            self.params.geometry.strike_length / 2.0
+        )
+        v_1 = self.params.geometry.origin[1] - (self.params.geometry.dip_length / 2.0)
+        v_2 = self.params.geometry.origin[1] + (self.params.geometry.dip_length / 2.0)
+        w_1 = self.params.geometry.origin[2] - (self.params.geometry.width / 2.0)
+        w_2 = self.params.geometry.origin[2] + (self.params.geometry.width / 2.0)
 
         vertices = np.array(
             [
@@ -159,29 +157,21 @@ def workspace(self) -> Workspace:
 
     def _rotate(self, vertices: np.ndarray) -> np.ndarray:
         """Rotate vertices and adjust for reference point."""
-        theta = -1 * self.params.dip_direction
-        phi = -1 * self.params.dip
-        rotated_vertices = rotate_xyz(vertices, self.center, theta, phi)
+        theta = -1 * self.params.geometry.direction
+        phi = -1 * self.params.geometry.dip
+        rotated_vertices = rotate_xyz(vertices, self.params.geometry.origin, theta, phi)
 
         return rotated_vertices
 
     def mask(self, mesh: Octree) -> np.ndarray:
-        plate = PlateModel(
-            strike_length=self.params.strike_length,
-            dip_length=self.params.dip_length,
-            width=self.params.width,
-            direction=self.params.dip_direction,
-            dip=self.params.dip,
-            origin=self.center,
-        )
         rotations = [
-            z_rotation_matrix(np.deg2rad(self.params.dip_direction)),
-            x_rotation_matrix(np.deg2rad(self.params.dip)),
+            z_rotation_matrix(np.deg2rad(self.params.geometry.direction)),
+            x_rotation_matrix(np.deg2rad(self.params.geometry.dip)),
         ]
         rotated_centers = rotate_points(
-            mesh.centroids, origin=plate.origin, rotations=rotations
+            mesh.centroids, origin=self.params.geometry.origin, rotations=rotations
         )
-        return inside_plate(rotated_centers, plate)
+        return inside_plate(rotated_centers, self.params.geometry)
 
 
 class Body(Parametric):

simpeg_drivers/utils/synthetics/meshes.py

@@ -122,19 +122,20 @@ def get_octree_mesh(
     )
 
     if plate is not None:
-        # TODO Consolidate PlateOptions and PlateModel into a single class to avoid this redundancy
         plate_options = PlateOptions(
-            plate=1.0,  # thickness
-            width=plate.width,
-            strike_length=plate.strike_length,
-            dip_length=plate.dip_length,
-            dip=plate.dip,
-            dip_direction=plate.direction,
-            elevation=0,
+            plate_property=1.0,  # thickness
+            geometry=PlateModel(
+                strike_length=plate.strike_length,
+                dip_length=plate.dip_length,
+                width=plate.width,
+                direction=plate.direction,
+                dip=plate.dip,
+                easting=plate.origin[0],
+                northing=plate.origin[1],
+                elevation=0.0,
+            ),
         )
-        center = list(plate.origin)
-
-        plate = Plate(plate_options, center=center, workspace=survey.workspace)
+        plate = Plate(plate_options, workspace=survey.workspace)
         mesh = OctreeDriver.refine_tree_from_triangulation(
             mesh, plate.surface, levels=(4,), finalize=False
         )

simpeg_drivers/utils/synthetics/options.py

@@ -56,7 +56,9 @@ class ModelOptions(BaseModel):
         strike_length=40.0,
         dip_length=40.0,
         width=40.0,
-        origin=(0.0, 0.0, 10.0),
+        easting=0.0,
+        northing=0.0,
+        elevation=10.0,
     )
     name: str = "model"
 

tests/plate_simulation/models/events_test.py

@@ -9,6 +9,7 @@
 # '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
 
 import numpy as np
+from geoapps_utils.modelling.plates import PlateModel
 from geoh5py import Workspace
 from geoh5py.objects import Surface
 
@@ -87,13 +88,19 @@ def test_anomaly(tmp_path):
         _, octree = get_topo_mesh(workspace)
         params = PlateOptions(
             name="my plate",
-            plate=10.0,
-            elevation=-1.5,
-            width=10.0,
-            strike_length=10.0,
-            dip_length=1.0,
+            plate_property=10.0,
+            geometry=PlateModel(
+                easting=5.0,
+                northing=5.0,
+                elevation=-1.5,
+                width=10.0,
+                strike_length=10.0,
+                dip_length=1.0,
+                dip_direction=90,
+                dip=90,
+            ),
         )
-        plate = Plate(params, center=(5.0, 5.0, -1.5))
+        plate = Plate(params)
 
         anomaly = Anomaly(body=plate, value=10.0)
         event_map = {1: ("Background", 1.0)}