From ceeee0277e05e4db9f558d5fabbef33cf7120707 Mon Sep 17 00:00:00 2001 From: Lachlan Grose Date: Wed, 5 Jun 2024 10:19:44 +1000 Subject: [PATCH] fix: adding export methods --- LoopStructural/datatypes/_structured_grid.py | 3 + LoopStructural/datatypes/_surface.py | 2 +- LoopStructural/export/geoh5.py | 55 ++++++++++++++++++- .../modelling/core/geological_model.py | 51 ++++++++++++++++- LoopStructural/utils/_surface.py | 8 ++- examples/1_basic/plot_7_exporting.py | 46 ++++++++++++++++ 6 files changed, 161 insertions(+), 4 deletions(-) create mode 100644 examples/1_basic/plot_7_exporting.py diff --git a/LoopStructural/datatypes/_structured_grid.py b/LoopStructural/datatypes/_structured_grid.py index 11cdf10c..2487650e 100644 --- a/LoopStructural/datatypes/_structured_grid.py +++ b/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") diff --git a/LoopStructural/datatypes/_surface.py b/LoopStructural/datatypes/_surface.py index 1d91c049..cd3e6969 100644 --- a/LoopStructural/datatypes/_surface.py +++ b/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 diff --git a/LoopStructural/export/geoh5.py b/LoopStructural/export/geoh5.py index a4c869c4..c801a4fc 100644 --- a/LoopStructural/export/geoh5.py +++ b/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) diff --git a/LoopStructural/modelling/core/geological_model.py b/LoopStructural/modelling/core/geological_model.py index 2460dc6f..6b911eb2 100644 --- a/LoopStructural/modelling/core/geological_model.py +++ b/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}') diff --git a/LoopStructural/utils/_surface.py b/LoopStructural/utils/_surface.py index a88f605e..8dbe520b 100644 --- a/LoopStructural/utils/_surface.py +++ b/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 diff --git a/examples/1_basic/plot_7_exporting.py b/examples/1_basic/plot_7_exporting.py new file mode 100644 index 00000000..c58067e4 --- /dev/null +++ b/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')