update simpeg drivers for new plate reference (top center)

Author: benk-mira

simpeg_drivers/plate_simulation/models/options.py

@@ -12,6 +12,7 @@
 
 import numpy as np
 from geoapps_utils.modelling.plates import PlateModel
+from geoapps_utils.utils.locations import topo_drape_elevation
 from geoh5py.objects import Points
 from pydantic import (
     AliasChoices,
@@ -22,6 +23,7 @@
     field_validator,
     model_validator,
 )
+from scipy.spatial import cKDTree
 
 
 T = TypeVar("T")
@@ -64,16 +66,20 @@ def center(
 
         :param survey: geoh5py survey object for plate simulation.
         :param surface: Points-like object to reference plate depth from.
-        :param depth_offset: Additional offset to be added to the depth of the plate.
         """
 
-        xy = (
-            survey.vertices[:, 0].mean() + self.geometry.origin[0],
-            survey.vertices[:, 1].mean() + self.geometry.origin[1],
+        xyz = np.atleast_2d(
+            [
+                survey.vertices[:, 0].mean(),
+                survey.vertices[:, 1].mean(),
+                0,
+            ]
+        )
+        topo_at_center = topo_drape_elevation(
+            xyz, surface.vertices, method="linear", triangulation=surface.cells
         )
-        z_topo = topography_above_point(topography=surface, point=xy)  # TODO
 
-        return xy + (z_topo - self.geometry.elevation,)
+        return xyz[0, 0], xyz[0, 1], topo_at_center[0, 2] - self.geometry.elevation
 
 
 class OverburdenOptions(BaseModel):

simpeg_drivers/utils/synthetics/options.py

@@ -141,7 +141,9 @@ class ModelOptions(BaseModel):
         width=40.0,
         easting=0.0,
         northing=0.0,
-        elevation=10.0,
+        elevation=0.0,
+        dip=90.0,
+        direction=0.0,
     )
     name: str = "model"
 

tests/plate_simulation/models/params_test.py

@@ -14,6 +14,44 @@
 from geoh5py.objects import Points, Surface
 
 from simpeg_drivers.plate_simulation.models.options import PlateOptions
+from simpeg_drivers.utils.synthetics.driver import SyntheticsComponents
+from simpeg_drivers.utils.synthetics.options import (
+    SurveyOptions,
+    SyntheticsComponentsOptions,
+)
+
+
+def test_plate_options_center(tmp_path):
+    with Workspace(tmp_path / "test.geoh5") as workspace:
+        components = SyntheticsComponents(
+            geoh5=workspace,
+            options=SyntheticsComponentsOptions(
+                method="gravity",
+                refine_plate=True,
+                survey=SurveyOptions(
+                    center=(0.0, 0.0), n_stations=10, n_lines=10, drape=15.0
+                ),
+            ),
+        )
+
+        params = PlateOptions(
+            name="my plate",
+            plate_property=1.0,
+            geometry=PlateModel(
+                strike_length=40.0,
+                dip_length=80.0,
+                width=5.0,
+                easting=10.0,
+                northing=-10.0,
+                elevation=30.0,
+                direction=0.0,
+                dip=90.0,
+            ),
+            number=1,
+            spacing=10.0,
+        )
+        center = params.center(components.survey, components.topography)
+        assert np.allclose(center, [0.0, 0.0, 20], atol=7e-1)
 
 
 def test_plate_params(tmp_path):
@@ -27,16 +65,12 @@ def test_plate_params(tmp_path):
             width=20.0,
             easting=10.0,
             northing=10.0,
-            elevation=-100.0,
+            elevation=100.0,
             direction=0.0,
             dip=90.0,
         ),
-        reference="center",
         number=1,
         spacing=10.0,
-        relative_locations=True,
-        reference_surface="topography",
-        reference_type="mean",
     )
     assert params.spacing == 0.0
 
@@ -53,33 +87,4 @@ def test_plate_params(tmp_path):
     )
 
     center = params.center(survey, topography)
-    assert np.allclose(center, [0, 0, -300])
-
-    params.relative_locations = False
-    center = params.center(survey, topography)
-    assert np.allclose(center, [10, 10, -100])
-
-
-def test_plate_params_empty_reference():
-    params = PlateOptions(
-        name="my plate",
-        plate_property=1.0,
-        geometry=PlateModel(
-            strike_length=1500.0,
-            dip_length=400.0,
-            width=20.0,
-            easting=10.0,
-            northing=10.0,
-            elevation=-100.0,
-            direction=0.0,
-            dip=90.0,
-        ),
-        reference="center",
-        number=1,
-        spacing=10.0,
-        relative_locations=True,
-        reference_surface=None,
-        reference_type=None,
-    )
-    assert params.reference_surface == "topography"
-    assert params.reference_type == "mean"
+    assert np.allclose(center, [-10.0, -10.0, -100])

tests/plate_simulation/models/plates_test.py

@@ -58,14 +58,17 @@ def test_vertical_east_striking_plate(tmp_path):
             vertical_east_striking.extent[1, 2] - vertical_east_striking.extent[0, 2],
             500.0,
         )
-        assert (
-            vertical_east_striking.vertices[:, 0].mean() == 0.0  # pylint: disable=no-member
+        assert np.isclose(
+            vertical_east_striking.vertices[:, 0].mean(),
+            0.0,  # pylint: disable=no-member
         )
-        assert (
-            vertical_east_striking.vertices[:, 1].mean() == 0.0  # pylint: disable=no-member
+        assert np.isclose(
+            vertical_east_striking.vertices[:, 1].mean(),
+            0.0,  # pylint: disable=no-member
         )
-        assert (
-            vertical_east_striking.vertices[:, 2].mean() == 0.0  # pylint: disable=no-member
+        assert np.isclose(
+            vertical_east_striking.vertices[:, 2].mean(),
+            -250.0,  # pylint: disable=no-member
         )
 
 
@@ -88,7 +91,9 @@ def test_dipping_plates_all_quadrants(tmp_path):
                 )
 
                 plate = Plate(params)
-                surface = plate.surface(workspace)
+                surface = plate.surface(
+                    workspace, name=f"Plate (dip: {dip}, dir: {dip_direction}"
+                )
                 locs = rotate_xyz(surface.vertices, [0.0, 0.0, 0.0], dip_direction, 0.0)
                 locs = rotate_xyz(locs, [0.0, 0.0, 0.0], 0.0, dip - 90.0)
                 assert np.allclose(locs, reference.vertices)
@@ -97,13 +102,13 @@ def test_dipping_plates_all_quadrants(tmp_path):
 def test_replicate_even(tmp_path):
     with Workspace(tmp_path / "test.geoh5") as workspace:
         options = PlateModel(
-            strike_length=1.0,
-            dip_length=1.0,
-            width=1.0,
+            strike_length=2.0,
+            dip_length=2.0,
+            width=2.0,
             direction=0.0,
             dip=0.0,
             easting=0.0,
-            northing=0.0,
+            northing=-1.0,
             elevation=0.0,
         )
         plate = Plate(options)
@@ -123,13 +128,13 @@ def test_replicate_even(tmp_path):
 def test_replicate_odd(tmp_path):
     with Workspace(tmp_path / "test.geoh5") as workspace:
         options = PlateModel(
-            strike_length=1.0,
-            dip_length=1.0,
-            width=1.0,
+            strike_length=2.0,
+            dip_length=2.0,
+            width=2.0,
             direction=0.0,
             dip=0.0,
             easting=0.0,
-            northing=0.0,
+            northing=-1.0,
             elevation=0.0,
         )
         plate = Plate(options)

tests/plate_simulation/runtest/driver_test.py

@@ -85,12 +85,7 @@ def test_plate_simulation_params_from_input_file(tmp_path):
         ifile.data["dip_direction"] = 0.0
         ifile.data["number"] = 9
         ifile.data["spacing"] = 10.0
-        ifile.data["relative_locations"] = True
-        ifile.data["easting"] = 10.0
-        ifile.data["northing"] = 10.0
         ifile.data["elevation"] = -250
-        ifile.data["reference_surface"] = "topography"
-        ifile.data["reference_type"] = "mean"
 
     params = PlateSimulationOptions.build(ifile)
     assert isinstance(params.simulation, SimPEGGroup)
@@ -130,7 +125,4 @@ def test_plate_simulation_params_from_input_file(tmp_path):
 
     assert params.model.plate_options.number == 9
     assert params.model.plate_options.spacing == 10.0
-    assert params.model.plate_options.relative_locations
-    assert params.model.plate_options.geometry.easting == 10.0
-    assert params.model.plate_options.geometry.northing == 10.0
     assert params.model.plate_options.geometry.elevation == -250.0

tests/plate_simulation/runtest/gravity_test.py

@@ -60,17 +60,14 @@ def test_gravity_plate_simulation(tmp_path):
         plate_params = PlateOptions(
             name="plate",
             geometry=PlateModel(
-                easting=0.0,
-                northing=0.0,
-                elevation=-250.0,
+                elevation=100.0,
                 width=100.0,
                 strike_length=100.0,
                 dip_length=100.0,
                 dip=0.0,
                 direction=0.0,
             ),
             plate_property=0.5,
-            reference="center",
         )
 
         model_params = ModelOptions(

tests/plate_simulation/runtest/sweep_test.py

@@ -56,12 +56,7 @@ def setup_plate_sweep(workspace) -> SimPEGGroup:
     options_dict["padding_distance"]["value"] = 1500.0
     options_dict["dip_direction"]["value"] = 0.0
     options_dict["number"]["value"] = 1
-    options_dict["relative_locations"]["value"] = True
-    options_dict["easting"]["value"] = 10.0
-    options_dict["northing"]["value"] = 10.0
-    options_dict["elevation"]["value"] = -250.0
-    options_dict["reference_surface"]["value"] = "topography"
-    options_dict["reference_type"]["value"] = "mean"
+    options_dict["elevation"]["value"] = 100.0
     options_dict["out_group"]["value"] = str(simulation.uid)
     simulation.options = options_dict