[ENH] Improve handling and visualization of structural elements

Improved the implementation for handling structural elements by differentiating the original data from its copy and providing appropriate error messages if data is not found. Additionally, changes were made to visualization of structural elements, including adding color specifications and implementing a 3D plot of the geological model.

Author: Leguark

gempy/API/initialization_API.py

@@ -101,9 +101,9 @@ def structural_elements_from_borehole_set(
     component_lith: dict[Hashable, np.ndarray] = borehole_set.get_top_coords_for_each_lith()
 
     for name, properties in elements_dict.items():
-        top_coordinates = component_lith.get(properties['top_lith'])
+        top_coordinates = component_lith.get(properties['id'])
         if top_coordinates is None:
-            raise ValueError(f"Top lithology {properties['top_lith']} not found in borehole set.")
+            raise ValueError(f"Top lithology {properties['id']} not found in borehole set.")
 
         element = StructuralElement(
             name=name,
@@ -119,6 +119,9 @@ def structural_elements_from_borehole_set(
             orientations=OrientationsTable(np.zeros(0, dtype=OrientationsTable.dt))
         )
         elements.append(element)
+    # Reverse the list to have the oldest rocks at the bottom
+    elements.reverse()
+    
     return elements
 
 

gempy/core/data/geo_model.py

@@ -161,7 +161,7 @@ def solutions(self, value):
     @property
     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
+        surface_points_table = self.structural_frame.surface_points_copy
         if self.transform is not None:
             surface_points_table.model_transform = self.transform
         return surface_points_table
@@ -178,7 +178,7 @@ def surface_points(self, value):
     @property
     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
+        orientations_table = self.structural_frame.orientations_copy
         if self.transform is not None:
             orientations_table.model_transform = self.transform
         return orientations_table

gempy/core/data/structural_frame.py

@@ -329,11 +329,16 @@ def elements_ids(self) -> np.ndarray:
         return np.arange(len(self.structural_elements)) + 1
 
     @property
-    def surface_points(self) -> SurfacePointsTable:
+    def surface_points_copy(self) -> SurfacePointsTable:
         """Returns a SurfacePointsTable for all surface points across the structural elements. This is a copy!"""
         all_data: np.ndarray = np.concatenate([element.surface_points.data for element in self.structural_elements])
         return SurfacePointsTable(data=all_data, name_id_map=self.element_name_id_map)
 
+    @property
+    def surface_points(self):
+        raise AttributeError("This property can only be set, not read. You can access the copy with `surface_points_copy` or"
+                             "the original on the individual structural elements.")
+
     @surface_points.setter
     def surface_points(self, modified_surface_points: SurfacePointsTable) -> None:
         """Distributes the modified surface points back to the structural elements."""
@@ -344,11 +349,16 @@ def surface_points(self, modified_surface_points: SurfacePointsTable) -> None:
             start += length
 
     @property
-    def orientations(self) -> OrientationsTable:
+    def orientations_copy(self) -> OrientationsTable:
         """Returns an OrientationsTable for all orientations across the structural elements."""
         all_data: np.ndarray = np.concatenate([element.orientations.data for element in self.structural_elements])
         return OrientationsTable(data=all_data)
 
+    @property
+    def orientations(self) -> OrientationsTable:
+        raise AttributeError("This property can only be set, not read. You can access the copy with `orientations_copy` or"
+                             "the original on the individual structural elements.")
+    
     @orientations.setter
     def orientations(self, modified_orientations: OrientationsTable) -> None:
         """Distributes the modified orientations back to the structural elements."""

test/test_modules/test_pile/test_stratigraphic_pile.py

@@ -9,9 +9,10 @@
 from subsurface.core.geological_formats.boreholes.survey import Survey
 from subsurface.core.reader_helpers.readers_data import GenericReaderFilesHelper
 from subsurface.modules.reader.wells.read_borehole_interface import read_lith, read_survey, read_collar
-from subsurface.modules.visualization import to_pyvista_line, pv_plot
+from subsurface.modules.visualization import to_pyvista_line, pv_plot, to_pyvista_points
 
 import gempy as gp
+import gempy_viewer as gpv
 
 
 # @pytest.mark.skip(reason="Not implemented yet")
@@ -64,7 +65,8 @@ def borehole_set(self):
 
         return borehole_set
 
-    def test_structural_elements(self, borehole_set: BoreholeSet):
+    # TODO: Rename this to test structural elements from borehole set
+    def test_structural_elements_from_borehole_set(self, borehole_set: BoreholeSet):
         from subsurface import LineSet
         borehole_trajectory: LineSet = borehole_set.combined_trajectory
         if PLOT := False:
@@ -78,16 +80,29 @@ def test_structural_elements(self, borehole_set: BoreholeSet):
         elements = gp.structural_elements_from_borehole_set(
             borehole_set=borehole_set,
             elements_dict={
-                # "Pleistozen": {"id": 10_000, "color": "#f9f97f", "top_lith": 10_000},
-                # "Kreide": {"id": 30_000, "color": "#a6d84a", "top_lith": 30_000},
-                # "Trias": {"id": 50_000, "color": "#a4469f", "top_lith": 50_000},
-                # "Perm": {"id": 60_000, "color": "#f4a142", "top_lith": 60_000},
-                "Rotliegend": {"id": 62_000, "color": "#bb825b", "top_lith": 62_000},
-                # "Devon": {"id": 80_000, "color": "#969594", "top_lith": 80_000}
+                    "Buntsandstein"       : {
+                            "id"   : 53_300,
+                            "color": "#983999"
+                    },
+                    "Werra-Anhydrit"      : {
+                            "id"   : 61_730,
+                            "color": "#00923f"
+                    },
+                    "Kupfershiefer"       : {
+                            "id"   : 61_760,
+                            "color": "#da251d"
+                    },
+                    "Zechsteinkonglomerat": {
+                            "id"   : 61_770,
+                            "color": "#f8c300"
+                    },
+                    "Rotliegend"          : {
+                            "id"   : 62_000,
+                            "color": "#bb825b"
+                    }
             }
         )
-        
-        
+
         group = gp.data.StructuralGroup(
             name="Stratigraphic Pile",
             elements=elements,
@@ -99,11 +114,9 @@ def test_structural_elements(self, borehole_set: BoreholeSet):
         )
         print(group)
 
-
-        component_lith = borehole_set.get_top_coords_for_each_lith()
-        rotliegend_xyz = component_lith[62_000]
-        extent_from_data = rotliegend_xyz.min(axis=0), rotliegend_xyz.max(axis=0)
-
+        all_surface_points_coords: gp.data.SurfacePointsTable = structural_frame.surface_points_copy
+        extent_from_data = all_surface_points_coords.xyz.min(axis=0), all_surface_points_coords.xyz.max(axis=0)
+        
         geo_model = gp.data.GeoModel(
             name="Stratigraphic Pile",
             structural_frame=structural_frame,
@@ -117,9 +130,16 @@ def test_structural_elements(self, borehole_set: BoreholeSet):
                 mesh_extraction=True,
                 number_octree_levels=3,
             ),
-
         )
-        
-        import gempy_viewer as gpv
-        gpv.plot_3d(geo_model)
-        pass
+        gempy_plot = gpv.plot_3d(
+            model=geo_model,
+            # ve=10,
+            kwargs_pyvista_bounds={
+                    'show_xlabels': False,
+                    'show_ylabels': False,
+                    # 'show_zlabels': True,
+            },
+            show=True,
+            image=True
+        )
+