fix: adding omf export

Author: lachlangrose

LoopStructural/datatypes/_point.py

@@ -67,6 +67,10 @@ def save(self, filename: Union[str, io.StringIO], ext=None):
                 for k, v in self.properties.items():
                     df[k] = v
             df.to_csv(filename, index=False)
+        elif ext == 'omf':
+            from LoopStructural.export.omf_wrapper import add_pointset_to_omf
+
+            add_pointset_to_omf(self, filename)
         else:
             raise ValueError(f'Unknown file extension {ext}')
 
@@ -154,5 +158,9 @@ def save(self, filename):
                 for k, v in self.properties.items():
                     df[k] = v
             df.to_csv(filename)
+        elif ext == 'omf':
+            from LoopStructural.export.omf_wrapper import add_pointset_to_omf
+
+            add_pointset_to_omf(self, filename)
         else:
             raise ValueError(f'Unknown file extension {ext}')

LoopStructural/datatypes/_surface.py

@@ -2,6 +2,7 @@
 from typing import Optional
 import numpy as np
 import io
+from pathlib import Path
 
 
 @dataclass
@@ -128,9 +129,8 @@ def from_dict(cls, d, flatten=False):
             d.get('cell_properties', None),
         )
 
-    def save(self, filename, ext=None):
-
-        print(filename, ext)
+    def save(self, filename, replace_spaces=True, ext=None):
+        filename = filename.replace(' ', '_') if replace_spaces else filename
         if isinstance(filename, (io.StringIO, io.BytesIO)):
             if ext is None:
                 raise ValueError('Please provide an extension for StringIO')
@@ -139,14 +139,12 @@ def save(self, filename, ext=None):
             filename = str(filename)
             if ext is None:
                 ext = filename.split('.')[-1].lower()
-
         if ext == 'json':
             import json
 
             with open(filename, 'w') as f:
                 json.dump(self.to_dict(), f)
         elif ext == 'vtk':
-            print(filename)
             self.vtk().save(filename)
         elif ext == 'obj':
             import meshio
@@ -179,5 +177,9 @@ def save(self, filename, ext=None):
                 for k, v in self.properties.items():
                     df[k] = v
             df.to_csv(filename, index=False)
+        elif ext == 'omf':
+            from LoopStructural.export.omf_wrapper import add_surface_to_omf
+
+            add_surface_to_omf(self, filename)
         else:
             raise ValueError(f"Extension {ext} not supported")

LoopStructural/export/omf_wrapper.py

@@ -0,0 +1,89 @@
+try:
+    import omf
+except ImportError:
+    raise ImportError(
+        "You need to install the omf package to use this feature. "
+        "You can install it with: pip install --pre omf"
+    )
+import numpy as np
+
+
+def get_project(filename):
+    try:
+        project = omf.load(filename)
+    except FileNotFoundError:
+        project = omf.Project(name='LoopStructural Model')
+    return project
+
+
+def get_cell_attributes(loopobject):
+    attributes = []
+    if loopobject.cell_properties:
+        attributes += [
+            omf.NumericAttribute(name=k, array=v, location="faces")
+            for k, v in loopobject.cell_properties.items()
+        ]
+    return attributes
+
+
+def get_point_attributed(loopobject):
+    attributes = []
+    if loopobject.properties:
+        attributes += [
+            omf.NumericAttribute(name=k, array=v, location="vertices")
+            for k, v in loopobject.properties.items()
+        ]
+    return attributes
+
+
+def add_surface_to_omf(surface, filename):
+
+    attributes = []
+    attributes += get_cell_attributes(surface)
+    attributes += get_point_attributed(surface)
+    surface = omf.Surface(
+        vertices=surface.vertices,
+        triangles=surface.triangles,
+        attributes=attributes,
+        name=surface.name,
+    )
+    project = get_project(filename)
+
+    project.elements += [surface]
+    omf.save(project, filename, mode='w')
+
+
+def add_pointset_to_omf(points, filename):
+
+    attributes = []
+    attributes += get_point_attributed(points)
+
+    points = omf.PointSet(
+        vertices=points.locations,
+        attributes=attributes,
+        name=points.name,
+    )
+
+    project = get_project(filename)
+    project.elements += [points]
+    omf.save(project, filename, mode='w')
+
+
+def add_structured_grid_to_omf(grid, filename):
+    print('Open Mining Format cannot store structured grids')
+    return
+    # attributes = []
+    # attributes += get_cell_attributes(grid)
+    # attributes += get_point_attributed(grid)
+
+    # vol = omf.TensorGridBlockModel(
+    #     name=grid.name,
+    #     tensor_u=np.ones(grid.nsteps[0]).astype(float),
+    #     tensor_v=np.ones(grid.nsteps[0]).astype(float),
+    #     tensor_w=np.ones(grid.nsteps[0]).astype(float),
+    #     origin=grid.origin,
+    #     attributes=attributes,
+    # )
+    # project = get_project(filename)
+    # project.elements += [vol]
+    # omf.save(project, filename, mode='w')

LoopStructural/modelling/core/geological_model.py

@@ -1810,7 +1810,7 @@ def get_stratigraphic_surfaces(self, units: List[str] = [], bottoms: bool = True
     def get_block_model(self, name='block model'):
         grid = self.bounding_box.structured_grid(name=name)
 
-        grid.properties_cell['stratigraphy'] = self.evaluate_model(
+        grid.cell_properties['stratigraphy'] = self.evaluate_model(
             self.bounding_box.cell_centers(), scale=False
         )
         return grid, self.stratigraphic_ids()
@@ -1831,19 +1831,19 @@ def save(
         if fault_surfaces:
             for s in self.get_fault_surfaces():
                 ## geoh5 can save everything into the same file
-                if extension == ".geoh5":
+                if extension == ".geoh5" or extension == '.omf':
                     s.save(filename)
                 else:
                     s.save(f'{name}_{s.name}.{extension}')
         if stratigraphic_surfaces:
             for s in self.get_stratigraphic_surfaces():
-                if extension == ".geoh5":
+                if extension == ".geoh5" or extension == '.omf':
                     s.save(filename)
                 else:
                     s.save(f'{name}_{s.name}.{extension}')
         if block_model:
             grid, ids = self.get_block_model()
-            if extension == ".geoh5":
+            if extension == ".geoh5" or extension == '.omf':
                 grid.save(filename)
             else:
                 grid.save(f'{name}_block_model.{extension}')
@@ -1852,15 +1852,16 @@ def save(
                 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)
+                    for data in self.__getitem__(group).get_data():
+                        if extension == ".geoh5" or extension == '.omf':
+                            data.save(filename)
                         else:
-                            self.__getitem__(series).save(f'{name}_{series}.{extension}')
+                            data.save(f'{name}_{group}_data.{extension}')
         if fault_data:
             for f in self.fault_names():
                 for d in self.__getitem__(f).get_data():
-                    if extension == ".geoh5":
+                    if extension == ".geoh5" or extension == '.omf':
+
                         d.save(filename)
                     else:
                         d.save(f'{name}_{group}.{extension}')

LoopStructural/modelling/features/_geological_feature.py

@@ -90,7 +90,7 @@ def __setitem__(self, key, item):
     def set_model(self, model):
         self.model = model
 
-    def evaluate_value(self, pos: np.ndarray, ignore_regions=False) -> np.ndarray:
+    def evaluate_value(self, pos: np.ndarray, ignore_regions=False, fillnan=None) -> np.ndarray:
         """
         Evaluate the scalar field value of the geological feature at the locations
         specified
@@ -122,6 +122,13 @@ def evaluate_value(self, pos: np.ndarray, ignore_regions=False) -> np.ndarray:
             logger.error(f"Unable to evaluate value for {self.name}")
         else:
             v[mask] = self.interpolator.evaluate_value(evaluation_points[mask, :])
+        if fillnan == 'nearest':
+            import scipy.spatial as spatial
+
+            nanmask = np.isnan(v)
+            tree = spatial.cKDTree(evaluation_points[~nanmask, :])
+            _d, i = tree.query(evaluation_points[nanmask, :])
+            v[nanmask] = v[~nanmask][i]
         return v
 
     def evaluate_gradient(self, pos: np.ndarray, ignore_regions=False) -> np.ndarray: