[WIP/ENH] Dealing with the RegularGrid transform during interpolation and output

Author: Leguark

gempy/core/data/xsolution.py renamed to gempy/core/data/DEP_xsolution.py

@@ -148,9 +148,16 @@ def solutions(self, value):
 
         # * Set solutions per element
         for e, element in enumerate(self.structural_frame.structural_elements[:-1]):  # * Ignore basement
-            dc_mesh = self._solutions.dc_meshes[e] if self._solutions.dc_meshes is not None else None
-            # TODO: These meshes are in the order of the scalar field
-            element.vertices = (self.transform.apply_inverse(dc_mesh.vertices) if dc_mesh is not None else None)
+            if self._solutions.dc_meshes is None:
+                continue
+            dc_mesh = self._solutions.dc_meshes[e]
+            if dc_mesh is None:
+                continue
+                
+            # TODO: These meshes are in the order of the scalar field 
+            world_coord_vertices = self.transform.apply_inverse(dc_mesh.vertices)
+            
+            element.vertices = world_coord_vertices
             element.edges = (dc_mesh.edges if dc_mesh is not None else None)
 
         # * Reordering the elements according to the scalar field

gempy/core/data/geo_model.py

@@ -1,3 +1,5 @@
+import warnings
+
 import dataclasses
 import enum
 import numpy as np
@@ -11,13 +13,13 @@
 @dataclasses.dataclass
 class Grid:
     class GridTypes(enum.Flag):
-        OCTREE = 2**0
-        DENSE = 2**1
-        CUSTOM = 2**2
-        TOPOGRAPHY = 2**3
-        SECTIONS = 2**4
-        CENTERED = 2**5
-        NONE = 2**10
+        OCTREE = 2 ** 0
+        DENSE = 2 ** 1
+        CUSTOM = 2 ** 2
+        TOPOGRAPHY = 2 ** 3
+        SECTIONS = 2 ** 4
+        CENTERED = 2 ** 5
+        NONE = 2 ** 10
 
     # ? What should we do with the extent?
     _extent: Optional[np.ndarray]  # * Model extent should be cross grid
@@ -31,11 +33,11 @@ class GridTypes(enum.Flag):
 
     values: np.ndarray = np.empty((0, 3))
     length: np.ndarray = np.empty(0)
-    
+
     _active_grids = GridTypes.NONE
-    
+
     _octree_levels: int = -1
-    
+
     def __init__(self, extent=None, resolution=None):
         self.extent = extent
         # Init basic grid empty
@@ -62,11 +64,11 @@ def extent(self):
                 raise AttributeError('Extent is not defined')
         else:
             return self._extent
-        
+
     @extent.setter
     def extent(self, value):
         self._extent = value
-        
+
     @property
     def bounding_box(self):
         extents = self.extent
@@ -80,81 +82,87 @@ def bounding_box(self):
                                         [extents[1], extents[3], extents[4]],  # max x, max y, min z
                                         [extents[1], extents[3], extents[5]]])  # max x, max y, max z
         return bounding_box_points
-        
-    
+
     @property
     def active_grids(self):
         return self._active_grids
-    
+
     @active_grids.setter
     def active_grids(self, value):
         self._active_grids = value
         self._update_values()
 
     @property
-    def dense_grid(self):
+    def dense_grid(self) -> RegularGrid:
         return self._dense_grid
-    
+
     @dense_grid.setter
     def dense_grid(self, value):
         self._dense_grid = value
         self.active_grids |= self.GridTypes.DENSE
         self._update_values()
-    
+
     @property
     def octree_grid(self):
         return self._octree_grid
-    
+
     @octree_grid.setter
     def octree_grid(self, value):
         self._octree_grid = value
         self.active_grids |= self.GridTypes.OCTREE
         self._update_values()
-        
+
     @property
     def custom_grid(self):
         return self._custom_grid
-    
+
     @custom_grid.setter
     def custom_grid(self, value):
         self._custom_grid = value
         self.active_grids |= self.GridTypes.CUSTOM
         self._update_values()
-        
+
     @property
     def topography(self):
         return self._topography
-    
+
     @topography.setter
     def topography(self, value):
         self._topography = value
         self.active_grids |= self.GridTypes.TOPOGRAPHY
         self._update_values()
-        
+
     @property
     def sections(self):
         return self._sections
-    
+
     @sections.setter
     def sections(self, value):
         self._sections = value
         self.active_grids |= self.GridTypes.SECTIONS
         self._update_values()
-        
+
     @property
     def centered_grid(self):
         return self._centered_grid
-    
+
     @centered_grid.setter
     def centered_grid(self, value):
         self._centered_grid = value
         self.active_grids |= self.GridTypes.CENTERED
         self._update_values()
-    
-    
+
     @property
     def regular_grid(self):
-        raise AttributeError('This property is deprecated. Use the dense_grid or octree_grid instead')
+        warnings.warn('This property is deprecated. Use the dense_grid or octree_grid instead', DeprecationWarning)
+        if self.dense_grid is not None and self.octree_grid is not None:
+            raise AttributeError('Both dense_grid and octree_grid are active. This is not possible.')
+        elif self.dense_grid is not None:
+            return self.dense_grid
+        elif self.octree_grid is not None:
+            return self.octree_grid
+        else:
+            return None
 
     @regular_grid.setter
     def regular_grid(self, value):
@@ -164,6 +172,10 @@ def regular_grid(self, value):
     def octree_levels(self):
         return self._octree_levels
 
+    @property
+    def octree_levels(self):
+        return self._octree_levels
+
     @octree_levels.setter
     def octree_levels(self, value):
         self._octree_levels = value
@@ -175,7 +187,7 @@ def octree_levels(self, value):
 
     def _update_values(self):
         values = []
-        
+
         if self.GridTypes.OCTREE in self.active_grids:
             values.append(self.octree_grid.values)
         if self.GridTypes.DENSE in self.active_grids:
@@ -188,12 +200,11 @@ def _update_values(self):
             values.append(self.sections.values)
         if self.GridTypes.CENTERED in self.active_grids:
             values.append(self.centered_grid.values)
-        
+
         self.values = np.concatenate(values)
 
         return self.values
 
-  
     def get_section_args(self, section_name: str):
         # TODO: This method should be part of the sections
         # assert type(section_name) is str, 'Only one section type can be retrieved'

gempy/core/data/grid.py

@@ -17,7 +17,7 @@ class RegularGrid:
 
     """
     resolution: np.ndarray
-    extent: np.ndarray  #: this is the orthogonal extent. If the grid is rotated, the extent will be different
+    extent: np.ndarray  #: this is the ORTHOGONAL extent. If the grid is rotated, the extent will be different
     values: np.ndarray
     mask_topo: np.ndarray
     x: Optional[np.ndarray]

gempy/core/data/grid_modules/grid_types.py

@@ -49,7 +49,8 @@ def interpolation_input_from_structural_frame(structural_frame: StructuralFrame,
     if grid.GridTypes.DENSE in grid.active_grids:
         regular_grid = engine_grid.RegularGrid(
             extent=new_extents,
-            regular_grid_shape=grid.dense_grid.resolution
+            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:
@@ -70,7 +71,8 @@ def interpolation_input_from_structural_frame(structural_frame: StructuralFrame,
 
     octree_grid = engine_grid.RegularGrid(
         extent=new_extents,
-        regular_grid_shape=np.array([2, 2, 2])
+        regular_grid_shape=np.array([2, 2, 2]),
+        transform=grid.dense_grid.transform # ! That here we need to grab the dense grid smells
     )
 
     grid: engine_grid.EngineGrid = engine_grid.EngineGrid(  # * Here we convert the GemPy grid to the