fix: adding export methods

Author: lachlangrose

LoopStructural/datatypes/_structured_grid.py

@@ -39,3 +39,6 @@ def vtk(self):
         )
         grid[self.name] = self.data.flatten(order="F")
         return grid
+
+    def save(self, filename):
+        raise NotImplementedError("Saving structured grids is not yet implemented")

LoopStructural/datatypes/_surface.py

@@ -95,7 +95,7 @@ def save(self, filename):
             with open(filename, 'w') as f:
                 json.dump(self.to_dict(), f)
         elif ext == 'vtk':
-            self.vtk.save(filename)
+            self.vtk().save(filename)
         elif ext == 'obj':
             import meshio
 

LoopStructural/export/geoh5.py

@@ -1,4 +1,5 @@
 import geoh5py
+import numpy as np
 
 
 def add_surface_to_geoh5(filename, surface, overwrite=True, groupname="Loop"):
@@ -27,4 +28,56 @@ def add_surface_to_geoh5(filename, surface, overwrite=True, groupname="Loop"):
         surface.add_data(data)
 
 
-# def add_points_to_geoh5(filename, points, overwrite=True, groupname="Loop"):
+def add_points_to_geoh5(filename, points, overwrite=True, groupname="Loop"):
+    with geoh5py.workspace.Workspace(filename) as workspace:
+        group = workspace.get_entity(groupname)[0]
+        if not group:
+            group = geoh5py.groups.ContainerGroup.create(
+                workspace, name=groupname, allow_delete=True
+            )
+        for point in points:
+            if point.name in workspace.list_entities_name.values():
+                existing_point = workspace.get_entity(point.name)
+                existing_point[0].allow_delete = True
+                if overwrite:
+                    workspace.remove_entity(existing_point[0])
+            data = {}
+            if point.properties is not None:
+                for k, v in point.properties.items():
+                    data[k] = {'association': "VERTEX", "values": v}
+            point = geoh5py.objects.Points.create(
+                workspace,
+                name=point.name,
+                vertices=point.vertices,
+                parent=group,
+            )
+            point.add_data(data)
+
+
+def add_block_model_to_geoh5py(filename, block_model, overwrite=True, groupname="Loop"):
+    with geoh5py.workspace.Workspace(filename) as workspace:
+        group = workspace.get_entity(groupname)[0]
+        if not group:
+            group = geoh5py.groups.ContainerGroup.create(
+                workspace, name=groupname, allow_delete=True
+            )
+        if block_model.name in workspace.list_entities_name.values():
+            existing_block = workspace.get_entity(block_model.name)
+            existing_block[0].allow_delete = True
+            if overwrite:
+                workspace.remove_entity(existing_block[0])
+        data = {}
+        if block_model.properties is not None:
+            for k, v in block_model.properties.items():
+                data[k] = {'association': "CELL", "values": v}
+        block = geoh5py.objects.BlockModel.create(
+            workspace,
+            name=block_model.name,
+            origin=[25, -100, 50],
+            u_cell_delimiters=np.cumsum(np.ones(16) * 5),  # Offsets along u
+            v_cell_delimiters=np.cumsum(np.ones(32) * 5),  # Offsets along v
+            z_cell_delimiters=np.cumsum(np.ones(16) * -2.5),  # Offsets along z (down)
+            rotation=30.0,
+            parent=group,
+        )
+        block.add_data(data)

LoopStructural/modelling/core/geological_model.py

@@ -7,7 +7,7 @@
 import numpy as np
 import pandas as pd
 from typing import List
-
+import pathlib
 from ...modelling.features.fault import FaultSegment
 
 from ...modelling.features.builders import (
@@ -37,6 +37,7 @@
 from ...modelling.intrusions import IntrusionBuilder
 
 from ...modelling.intrusions import IntrusionFrameBuilder
+from LoopStructural.modelling.features import fault
 
 
 logger = getLogger(__name__)
@@ -1807,3 +1808,51 @@ def get_block_model(self):
             self.bounding_box.cell_centers(), scale=False
         )
         return grid, self.stratigraphic_ids()
+
+    def save(
+        self,
+        filename: str,
+        block_model: bool = False,
+        stratigraphic_surfaces=True,
+        fault_surfaces=True,
+        stratigraphic_data=True,
+        fault_data=True,
+    ):
+        path = pathlib.Path(filename)
+        extension = path.suffix
+        name = path.stem
+        stratigraphic_surfaces = self.get_stratigraphic_surfaces()
+        if fault_surfaces:
+            for s in self.get_fault_surfaces():
+                ## geoh5 can save everything into the same file
+                if extension == ".geoh5":
+                    s.save(filename)
+                else:
+                    s.save(f'{name}_{s.name}.{extension}')
+        if stratigraphic_surfaces:
+            for s in self.get_stratigraphic_surfaces():
+                if extension == ".geoh5":
+                    s.save(filename)
+                else:
+                    s.save(f'{name}_{s.name}.{extension}')
+        if block_model:
+            grid, ids = self.get_block_model()
+            if extension == ".geoh5":
+                grid.save(filename)
+            else:
+                grid.save(f'{name}_block_model.{extension}')
+        if stratigraphic_data:
+            for group in self.stratigraphic_column.keys():
+                if group == "faults":
+                    continue
+                for series in self.stratigraphic_column[group].keys():
+                    if extension == ".geoh5":
+                        self.__getitem__(series).save(filename)
+                    else:
+                        self.__getitem__(series).save(f'{name}_{series}.{extension}')
+        if fault_data:
+            for f in self.fault_names():
+                if extension == ".geoh5":
+                    self.__getitem__(f).save(filename)
+                else:
+                    self.__getitem__(f).save(f'{name}_{f}.{extension}')

LoopStructural/utils/_surface.py

@@ -102,7 +102,13 @@ def fit(
                 values,
             )
         logger.info(f'Isosurfacing at values: {isovalues}')
-        for isovalue in isovalues:
+        if name is None:
+            names = ["surface"] * len(isovalues)
+        if isinstance(name, str):
+            names = [name] * len(isovalues)
+        if isinstance(name, list):
+            names = name
+        for name, isovalue in zip(names, isovalues):
             try:
                 step_vector = (self.bounding_box.maximum - self.bounding_box.origin) / (
                     np.array(self.bounding_box.nsteps) - 1

examples/1_basic/plot_7_exporting.py

@@ -0,0 +1,46 @@
+"""
+
+1j. Exporting models
+===============================
+
+Models can be exported to vtk, gocad and geoh5 formats.
+"""
+
+from LoopStructural import GeologicalModel
+from LoopStructural.datasets import load_claudius
+
+data, bb = load_claudius()
+
+model = GeologicalModel(bb[0, :], bb[1, :])
+model.data = data
+model.create_and_add_foliation("strati")
+
+
+######################################################################
+# Export surfaces to vtk
+# ~~~~~~~~~~~~~~~~~~~~~~
+# Isosurfaces can be extracted from a geological feature by calling
+# the `.surfaces` method on the feature. The argument for this method
+# is the value, values or number of surfaces that are extracted.
+# This returns a list of `LoopStructural.datatypes.Surface` objects
+# These objects can be interrogated to return the triangles, vertices
+# and normals. Or can be exported into another format using the `save`
+# method. The supported file formats are `vtk`, `ts` and `geoh5`.
+#
+
+surfaces = model['strati'].surfaces(value=0.0)
+
+print(surfaces)
+
+print(surfaces[0].vtk)
+
+surfaces[0].save('text.geoh5')
+
+######################################################################
+# Export the model to geoh5
+# ~~~~~~~~~~~~~~~~~~~~~~~~~
+# The entire model can be exported to a geoh5 file using the `save_model`
+# method. This will save all the data, foliations, faults and other objects
+# in the model to a geoh5 file. This file can be loaded into LoopStructural
+
+model.save_model('model.geoh5')