[CLN] Refactor GeoModel.transform to input_transform

[CLN] Refactor extent to orthogonal_extent

[WIP] Trying to figure out how to apply the grid rotation

Author: Leguark

examples/tutorials/ch1_fundamentals/ch1_4_onlap_relations.py

@@ -47,7 +47,7 @@
 geo_model.structural_frame
 
 # %%
-geo_model.transform.apply_anisotropy(gp.data.GlobalAnisotropy.NONE)
+geo_model.input_transform.apply_anisotropy(gp.data.GlobalAnisotropy.NONE)
 gp.add_structural_group(
     model=geo_model,
     group_index=0,

examples/tutorials/ch1_fundamentals/ch1_5_fault_relations.py

@@ -63,7 +63,7 @@
 
 # %% 
 
-geo_model.transform.apply_anisotropy(gp.data.GlobalAnisotropy.NONE)
+geo_model.input_transform.apply_anisotropy(gp.data.GlobalAnisotropy.NONE)
 if False:
     gp.compute_model(geo_model)
     # %%

examples/tutorials/ch3-Interpolations/ch3_1_kriging_interpolation_and_simulation.py

@@ -93,7 +93,7 @@
 # creating a domain object from the gempy solution, a defined domain conditioning data
 domain = kriging.Domain(
     model_solutions=sol,
-    transform=geo_data.transform,
+    transform=geo_data.input_transform,
     domain=[2],
     data=cond_data
 )

gempy/core/data/geo_model.py

@@ -60,7 +60,7 @@ class GeoModel:
     _interpolation_options: InterpolationOptions  #: The interpolation options provided by the user.
     geophysics_input: GeophysicsInput = None  #: The geophysics input of the geological model.
 
-    transform: Transform = None  #: The transformation used in the geological model for input points.
+    input_transform: Transform = None  #: The transformation used in the geological model for input points.
 
     interpolation_grid: EngineGrid = None  #: Optional grid used for interpolation. Can be seen as a cache field.
     _interpolationInput: InterpolationInput = None  #: Input data for interpolation. Fed by the structural frame and can be seen as a cache field.
@@ -84,7 +84,7 @@ def __init__(self, name: str, structural_frame: StructuralFrame, grid: Grid, int
 
         self.grid = grid
         self._interpolation_options = interpolation_options
-        self.transform = Transform.from_input_points(
+        self.input_transform = Transform.from_input_points(
             surface_points=self.surface_points_copy,
             orientations=self.orientations_copy
         )
@@ -94,12 +94,12 @@ def __repr__(self):
         return pprint.pformat(self.__dict__)
 
     def update_transform(self, auto_anisotropy: GlobalAnisotropy = GlobalAnisotropy.NONE, anisotropy_limit: Optional[np.ndarray] = None):
-        self.transform = Transform.from_input_points(
+        self.input_transform = Transform.from_input_points(
             surface_points=self.surface_points_copy,
             orientations=self.orientations_copy
         )
 
-        self.transform.apply_anisotropy(anisotropy_type=auto_anisotropy, anisotropy_limit=anisotropy_limit)
+        self.input_transform.apply_anisotropy(anisotropy_type=auto_anisotropy, anisotropy_limit=anisotropy_limit)
 
     @property
     def interpolation_options(self) -> InterpolationOptions:
@@ -155,7 +155,7 @@ def solutions(self, value):
                 continue
 
             # TODO: These meshes are in the order of the scalar field 
-            world_coord_vertices = self.transform.apply_inverse(dc_mesh.vertices)
+            world_coord_vertices = self.input_transform.apply_inverse(dc_mesh.vertices)
 
             element.vertices = world_coord_vertices
             element.edges = (dc_mesh.edges if dc_mesh is not None else None)
@@ -170,8 +170,8 @@ def solutions(self, value):
     def surface_points_copy(self):
         """This is a copy! Returns a SurfacePointsTable for all surface points across the structural elements"""
         surface_points_table = self.structural_frame.surface_points_copy
-        if self.transform is not None:
-            surface_points_table.model_transform = self.transform
+        if self.input_transform is not None:
+            surface_points_table.model_transform = self.input_transform
         return surface_points_table
 
     @property
@@ -187,8 +187,8 @@ def surface_points(self, value):
     def orientations_copy(self) -> OrientationsTable:
         """This is a copy! Returns a OrientationsTable for all orientations across the structural elements"""
         orientations_table = self.structural_frame.orientations_copy
-        if self.transform is not None:
-            orientations_table.model_transform = self.transform
+        if self.input_transform is not None:
+            orientations_table.model_transform = self.input_transform
         return orientations_table
 
     @property
@@ -208,7 +208,7 @@ def interpolation_input_copy(self):
         self._interpolationInput = interpolation_input_from_structural_frame(
             structural_frame=self.structural_frame,
             grid=self.grid,
-            transform=self.transform
+            input_transform=self.input_transform
         )
 
         return self._interpolationInput

gempy/modules/data_manipulation/engine_factory.py

@@ -12,69 +12,75 @@
 
 
 def interpolation_input_from_structural_frame(structural_frame: StructuralFrame, grid: Grid,
-                                              transform: Transform) -> InterpolationInput:
+                                              input_transform: Transform) -> InterpolationInput:
     _legacy_factor = 0
 
     if LEGACY_COORDS := False:
         _legacy_factor = 0.5
 
+    total_transform: Transform = input_transform + grid.dense_grid.transform
+
     surface_points_copy = structural_frame.surface_points_copy
     surface_points: SurfacePoints = SurfacePoints(
-        sp_coords=transform.apply(surface_points_copy.xyz) + _legacy_factor,
+        sp_coords=total_transform.apply(surface_points_copy.xyz) + _legacy_factor,
         nugget_effect_scalar=surface_points_copy.nugget
     )
 
     orientations_copy = structural_frame.orientations_copy
     orientations: Orientations = Orientations(
-        dip_positions=transform.apply(orientations_copy.xyz) + _legacy_factor,
-        dip_gradients=transform.transform_gradient(orientations_copy.grads),
+        dip_positions=total_transform.apply(orientations_copy.xyz) + _legacy_factor,
+        dip_gradients=total_transform.transform_gradient(orientations_copy.grads),
         nugget_effect_grad=orientations_copy.nugget
     )
 
-    # region Transforming the grid
+    grid: engine_grid.EngineGrid = _apply_input_transform_to_grids(grid, input_transform)
+
 
-    transformed = transform.apply(grid.bounding_box)  # ? isn't this making the regular grid not optional?
+    interpolation_input: InterpolationInput = InterpolationInput(
+        surface_points=surface_points,
+        orientations=orientations,
+        grid=grid,
+        unit_values=structural_frame.elements_ids  # TODO: Here we will need to pass densities etc.
+    )
+
+    return interpolation_input
+
+
+def _apply_input_transform_to_grids(grid: Grid, input_transform: Transform) -> engine_grid.EngineGrid:
+    transformed = input_transform.apply(grid.bounding_box)  # ! grid already has the grid transform applied
     new_extents = np.array([transformed[:, 0].min(), transformed[:, 0].max(),
                             transformed[:, 1].min(), transformed[:, 1].max(),
                             transformed[:, 2].min(), transformed[:, 2].max()])
-
+    
     # Initialize all variables to None
     octree_grid: Optional[engine_grid.RegularGrid] = None
     regular_grid: Optional[engine_grid.RegularGrid] = None
     custom_values: Optional[engine_grid.GenericGrid] = None
     topography_values: Optional[engine_grid.GenericGrid] = None
     section_values: Optional[engine_grid.GenericGrid] = None
     centered_grid: Optional[engine_grid.CenteredGrid] = None
-
+    
     if grid.GridTypes.DENSE in grid.active_grids:
         regular_grid = engine_grid.RegularGrid(
-            extent=new_extents,
+            orthogonal_extent=new_extents,
             regular_grid_shape=grid.dense_grid.resolution,
-            transform=grid.dense_grid.transform
         )
-
     if grid.GridTypes.CUSTOM in grid.active_grids and grid.custom_grid is not None:
-        custom_values = engine_grid.GenericGrid(values=transform.apply(grid.custom_grid.values))
-
+        custom_values = engine_grid.GenericGrid(values=input_transform.apply(grid.custom_grid.values))
     if grid.GridTypes.TOPOGRAPHY in grid.active_grids and grid.topography is not None:
-        topography_values = engine_grid.GenericGrid(values=transform.apply(grid.topography.values))
-
+        topography_values = engine_grid.GenericGrid(values=input_transform.apply(grid.topography.values))
     if grid.GridTypes.SECTIONS in grid.active_grids and grid.sections is not None:
-        section_values = engine_grid.GenericGrid(values=transform.apply(grid.sections.values))
-
+        section_values = engine_grid.GenericGrid(values=input_transform.apply(grid.sections.values))
     if grid.GridTypes.CENTERED in grid.active_grids and grid.centered_grid is not None:
         centered_grid = engine_grid.CenteredGrid(
-            centers=transform.apply(grid.centered_grid.centers),
-            radius=transform.scale_points(np.atleast_2d(grid.centered_grid.radius))[0],
+            centers=input_transform.apply(grid.centered_grid.centers),
+            radius=input_transform.scale_points(np.atleast_2d(grid.centered_grid.radius))[0],
             resolution=grid.centered_grid.resolution
         )
-
     octree_grid = engine_grid.RegularGrid(
-        extent=new_extents,
-        regular_grid_shape=np.array([2, 2, 2]),
-        transform=grid.dense_grid.transform # ! That here we need to grab the dense grid smells
+        orthogonal_extent=new_extents,
+        regular_grid_shape=np.array([2, 2, 2])
     )
-
     grid: engine_grid.EngineGrid = engine_grid.EngineGrid(  # * Here we convert the GemPy grid to the
         octree_grid=octree_grid,  # BUG: Adapt the engine to deal with this
         dense_grid=regular_grid,
@@ -83,14 +89,4 @@ def interpolation_input_from_structural_frame(structural_frame: StructuralFrame,
         custom_grid=custom_values,
         geophysics_grid=centered_grid
     )
-
-    # endregion
-
-    interpolation_input: InterpolationInput = InterpolationInput(
-        surface_points=surface_points,
-        orientations=orientations,
-        grid=grid,
-        unit_values=structural_frame.elements_ids  # TODO: Here we will need to pass densities etc.
-    )
-
-    return interpolation_input
+    return grid

test/test_core/test_transfoms.py

@@ -7,8 +7,8 @@
 
 def test_transform_1():
     geo_data: GeoModel = _create_data()
-    print(geo_data.transform)
-    transformed_xyz = geo_data.transform.apply(geo_data.surface_points_copy.xyz)
+    print(geo_data.input_transform)
+    transformed_xyz = geo_data.input_transform.apply(geo_data.surface_points_copy.xyz)
     print(transformed_xyz)
     return
 

test/test_modules/test_faults/test_finite_faults.py

@@ -25,8 +25,8 @@ def test_finite_fault_scalar_field_on_fault():
     regular_grid = geo_model.grid.regular_grid
 
     # TODO: Extract grid from the model
-    scaled_center = geo_model.transform.apply(center.reshape(1, -1))[0]
-    scaled_radius = geo_model.transform.scale_points(radius.reshape(1, -1))[0]
+    scaled_center = geo_model.input_transform.apply(center.reshape(1, -1))[0]
+    scaled_radius = geo_model.input_transform.scale_points(radius.reshape(1, -1))[0]
     scalar_funtion: callable = gp.implicit_functions.ellipsoid_3d_factory(  # * This paints the 3d regular grid
         center=scaled_center,
         radius=scaled_radius,