Merge pull request #286 from Loop3D/export/gocad-vo

fix: adding gocad export for structured grid support

Author: lachlangrose

.github/workflows/tester.yml

@@ -33,7 +33,7 @@ jobs:
       - name: Installing dependencies
         shell: bash -l {0}
         run: |
-          conda install -c conda-forge numpy scipy scikit-image scikit-learn pytest networkx osqp matplotlib -y
+          conda install -c conda-forge numpy scipy scikit-image scikit-learn pyvista pandas pytest networkx osqp matplotlib -y
       - name: Building and install
         shell: bash -l {0}
         run: |

LoopStructural/datatypes/_structured_grid.py

@@ -28,6 +28,7 @@ class StructuredGrid:
     name : str, optional
         Name of the grid, by default "default_grid"
     """
+
     origin: np.ndarray = field(default_factory=lambda: np.array([0, 0, 0]))
     step_vector: np.ndarray = field(default_factory=lambda: np.array([1, 1, 1]))
     nsteps: np.ndarray = field(default_factory=lambda: np.array([10, 10, 10]))
@@ -60,9 +61,9 @@ def maximum(self):
         Returns
         -------
         np.ndarray
-            Maximum coordinates (origin + nsteps * step_vector)
+            Maximum coordinates (origin + (nsteps - 1) * step_vector)
         """
-        return self.origin + self.nsteps * self.step_vector
+        return self.origin + (self.nsteps-1) * self.step_vector
 
     def vtk(self):
         """Convert the structured grid to a PyVista RectilinearGrid.
@@ -134,15 +135,16 @@ def cell_centres(self):
             self.nsteps[2] - 1,
         )
         x, y, z = np.meshgrid(x, y, z, indexing="ij")
-        return np.vstack([x.flatten(order='f'), y.flatten(order='f'), z.flatten(order='f')]).T
+        return np.vstack([x.flatten(order="f"), y.flatten(order="f"), z.flatten(order="f")]).T
 
     @property
     def nodes(self):
+
         x = np.linspace(self.origin[0], self.maximum[0], self.nsteps[0])
         y = np.linspace(self.origin[1], self.maximum[1], self.nsteps[1])
         z = np.linspace(self.origin[2], self.maximum[2], self.nsteps[2])
         x, y, z = np.meshgrid(x, y, z, indexing="ij")
-        return np.vstack([x.flatten(order='f'), y.flatten(order='f'), z.flatten(order='f')]).T
+        return np.vstack([x.flatten(order="f"), y.flatten(order="f"), z.flatten(order="f")]).T
 
     def merge(self, other):
         if not np.all(np.isclose(self.origin, other.origin)):
@@ -157,36 +159,33 @@ def merge(self, other):
         for name, data in other.properties.items():
             self.properties[name] = data
 
-    def save(self, filename, *,group='Loop'):
+    def save(self, filename, *, group="Loop"):
         filename = str(filename)
-        ext = filename.split('.')[-1]
-        if ext == 'json':
+        ext = filename.split(".")[-1].lower()
+        if ext == "json":
             import json
 
-            with open(filename, 'w') as f:
+            with open(filename, "w") as f:
                 json.dump(self.to_dict(), f)
-        elif ext == 'vtk':
+        elif ext == "vtk":
             self.vtk().save(filename)
 
-        elif ext == 'geoh5':
+        elif ext == "geoh5":
             from LoopStructural.export.geoh5 import add_structured_grid_to_geoh5
 
             add_structured_grid_to_geoh5(filename, self, groupname=group)
-        elif ext == 'pkl':
+        elif ext == "pkl":
             import pickle
 
-            with open(filename, 'wb') as f:
+            with open(filename, "wb") as f:
                 pickle.dump(self, f)
-        elif ext == 'omf':
+        elif ext == "omf":
             from LoopStructural.export.omf_wrapper import add_structured_grid_to_omf
 
             add_structured_grid_to_omf(self, filename)
-        elif ext == 'vs':
-            raise NotImplementedError(
-                "Saving structured grids in gocad format is not yet implemented"
-            )
-            # from LoopStructural.export.gocad import _write_structued_grid
+        elif ext == "vo" or ext == "gocad":
+            from LoopStructural.export.gocad import _write_structured_grid_gocad
 
-            # _write_pointset(self, filename)
+            _write_structured_grid_gocad(self, filename)
         else:
-            raise ValueError(f'Unknown file extension {ext}')
+            raise ValueError(f"Unknown file extension {ext}")

LoopStructural/export/gocad.py

@@ -1,5 +1,153 @@
+import re
+from pathlib import Path
+
 import numpy as np
 
+from LoopStructural.utils import getLogger
+
+logger = getLogger(__name__)
+
+
+def _normalise_voxet_property(values, property_name, nsteps):
+    array = np.asarray(values)
+    expected_shape = tuple(int(step) for step in nsteps)
+    expected_size = int(np.prod(expected_shape))
+
+    if array.shape == expected_shape:
+        flat_values = array.reshape(-1, order="F")
+    else:
+        flat_values = np.squeeze(array)
+        if flat_values.shape == expected_shape:
+            flat_values = flat_values.reshape(-1, order="F")
+        elif flat_values.ndim == 1 and flat_values.size == expected_size:
+            flat_values = flat_values
+        else:
+            raise ValueError(
+                f"Property '{property_name}' must have shape {expected_shape} or size {expected_size}"
+            )
+
+    if np.issubdtype(flat_values.dtype, np.integer):
+        if flat_values.size == 0:
+            export_dtype = np.int8
+            storage_type = "Octet"
+            element_size = 1
+        elif flat_values.min() >= np.iinfo(np.int8).min and flat_values.max() <= np.iinfo(np.int8).max:
+            export_dtype = np.int8
+            storage_type = "Octet"
+            element_size = 1
+        else:
+            export_dtype = np.dtype(">i4")
+            storage_type = "Integer"
+            element_size = 4
+        no_data_value = None
+    elif np.issubdtype(flat_values.dtype, np.floating):
+        export_dtype = np.dtype(">f4")
+        storage_type = "Float"
+        element_size = 4
+        no_data_value = -999999.0
+        flat_values = np.nan_to_num(flat_values, nan=no_data_value)
+    else:
+        raise ValueError(f"Property '{property_name}' has unsupported dtype {flat_values.dtype}")
+
+    return {
+        "values": np.asarray(flat_values, dtype=export_dtype),
+        "storage_type": storage_type,
+        "element_size": element_size,
+        "no_data_value": no_data_value,
+    }
+
+
+def _write_structured_grid_gocad(grid, file_name):
+    """Write a StructuredGrid to GOCAD VOXET format."""
+    vo_path = Path(file_name).with_suffix(".vo")
+    axis_n = np.asarray(grid.nsteps, dtype=int)
+    property_source = grid.properties
+    axis_min = np.min(grid.nodes, axis=0)
+    axis_max = np.max(grid.nodes, axis=0)
+
+    if property_source:
+        if grid.cell_properties:
+            logger.warning(
+                "StructuredGrid GOCAD export uses point properties; cell_properties were not exported"
+            )
+    elif grid.cell_properties:
+        axis_n = np.asarray(grid.nsteps, dtype=int) - 1
+        if np.any(axis_n <= 0):
+            raise ValueError("StructuredGrid cell_properties require at least two grid nodes per axis")
+        property_source = grid.cell_properties
+        axis_min = np.min(grid.cell_centres, axis=0)
+        axis_max = np.max(grid.cell_centres, axis=0)
+    else:
+        raise ValueError("StructuredGrid has no properties to export to GOCAD")
+
+    export_properties = []
+    for index, (property_name, values) in enumerate(property_source.items(), start=1):
+        export_info = _normalise_voxet_property(values, property_name, axis_n)
+        safe_name = re.sub(r"[^0-9A-Za-z_-]+", "_", property_name).strip("_") or f"property_{index}"
+        data_path = vo_path.with_name(f"{vo_path.stem}_{safe_name}@@")
+        export_properties.append(
+            {
+                "index": index,
+                "name": safe_name,
+                "data_path": data_path,
+                **export_info,
+            }
+        )
+
+    with open(vo_path, "w") as fp:
+        fp.write(
+            f"""GOCAD Voxet 1
+HEADER {{
+name: {grid.name}
+}}
+GOCAD_ORIGINAL_COORDINATE_SYSTEM
+NAME Default
+AXIS_NAME \"X\" \"Y\" \"Z\"
+AXIS_UNIT \"m\" \"m\" \"m\"
+ZPOSITIVE Elevation
+END_ORIGINAL_COORDINATE_SYSTEM
+AXIS_O 0.000000 0.000000 0.000000
+AXIS_U 1.000000 0.000000 0.000000
+AXIS_V 0.000000 1.000000 0.000000
+AXIS_W 0.000000 0.000000 1.000000
+AXIS_MIN {axis_min[0]} {axis_min[1]} {axis_min[2]}
+AXIS_MAX {axis_max[0]} {axis_max[1]} {axis_max[2]}
+AXIS_N {axis_n[0]} {axis_n[1]} {axis_n[2]}
+AXIS_NAME \"X\" \"Y\" \"Z\"
+AXIS_UNIT \"m\" \"m\" \"m\"
+AXIS_TYPE even even even
+"""
+        )
+        for export_property in export_properties:
+            fp.write(
+                f"""PROPERTY {export_property['index']} {export_property['name']}
+PROPERTY_CLASS {export_property['index']} {export_property['name']}
+PROP_UNIT {export_property['index']} {export_property['name']}
+PROPERTY_CLASS_HEADER {export_property['index']} {export_property['name']} {{
+}}
+PROPERTY_SUBCLASS {export_property['index']} QUANTITY {export_property['storage_type']}
+"""
+            )
+            if export_property["no_data_value"] is not None:
+                fp.write(
+                    f"PROP_NO_DATA_VALUE {export_property['index']} {export_property['no_data_value']}\n"
+                )
+            fp.write(
+                f"""PROP_ETYPE {export_property['index']} IEEE
+PROP_FORMAT {export_property['index']} RAW
+PROP_ESIZE {export_property['index']} {export_property['element_size']}
+PROP_OFFSET {export_property['index']} 0
+PROP_FILE {export_property['index']} {export_property['data_path'].name}
+"""
+            )
+        fp.write("END\n")
+
+    for export_property in export_properties:
+        with open(export_property["data_path"], "wb") as fp:
+            export_property["values"].tofile(fp)
+
+    return True
+
 
 def _write_feat_surfs_gocad(surf, file_name):
     """
@@ -23,8 +171,6 @@ def _write_feat_surfs_gocad(surf, file_name):
     True if successful
 
     """
-    from pathlib import Path
-
     properties_header = None
     if surf.properties:
 
@@ -112,8 +258,6 @@ def _write_feat_surfs_gocad(surf, file_name):
 #     True if successful
 
 #     """
-#     from pathlib import Path
-
 #     file_name = Path(file_name).with_suffix(".vs")
 #     with open(f"{file_name}", "w") as fd:
 #         fd.write(
@@ -123,4 +267,4 @@ def _write_feat_surfs_gocad(surf, file_name):
 # }}
 # GOCAD_ORIGINAL_COORDINATE_SYSTEM
 # NAME Default
-# PROJECTION Unknown
+# PROJECTION Unknown

LoopStructural/interpolators/supports/_2d_structured_grid.py

@@ -378,12 +378,12 @@ def position_to_cell_corners(self, pos):
         ----------
         pos : np.array
             (N, 2) array of xy coordinates representing the positions of N points.
-    
+
         Returns
         -------
         globalidx : np.array
-            (N, 4) array of global indices corresponding to the 4 corner nodes of the cell 
-            each point lies in. If a point lies outside the support, its corresponding entry 
+            (N, 4) array of global indices corresponding to the 4 corner nodes of the cell
+            each point lies in. If a point lies outside the support, its corresponding entry
             will be set to -1.
         inside : np.array
             (N,) boolean array indicating whether each point is inside the support domain.
@@ -490,9 +490,36 @@ def position_to_cell_vertices(self, pos):
     def onGeometryChange(self):
         pass
 
-    def vtk(self, node_properties={}, cell_properties={}):
-        raise NotImplementedError("VTK output not implemented for structured grid")
-        pass
+    def vtk(self, node_properties=None, cell_properties=None, z=0.0):
+        """
+        Create a vtk unstructured grid from the mesh
+        """
+        if node_properties is None:
+            node_properties = {}
+        if cell_properties is None:
+            cell_properties = {}
+
+        try:
+            import pyvista as pv
+        except ImportError:
+            raise ImportError("pyvista is required for this functionality")
+
+        from pyvista import CellType
+
+        points = np.zeros((self.n_nodes, 3))
+        points[:, :2] = self.nodes
+        points[:, 2] = z
+        celltype = np.full(self.n_elements, CellType.QUAD, dtype=np.uint8)
+        vtk_elements = self.elements[:, [0, 1, 3, 2]]
+        elements = np.hstack(
+            [np.full((vtk_elements.shape[0], 1), 4, dtype=int), vtk_elements.astype(np.int64)]
+        ).ravel()
+        grid = pv.UnstructuredGrid(elements, celltype, points)
+        for key, value in node_properties.items():
+            grid.point_data[key] = value
+        for key, value in cell_properties.items():
+            grid.cell_data[key] = value
+        return grid
 
     def get_operators(self, weights: Dict[str, float]) -> Dict[str, Tuple[np.ndarray, float]]:
         """Get
@@ -509,8 +536,8 @@ def get_operators(self, weights: Dict[str, float]) -> Dict[str, Tuple[np.ndarray
         """
         # in a map we only want the xy operators
         operators = {
-            'dxy': (Operator.Dxy_mask[1, :, :], weights['dxy'] * 2),
-            'dxx': (Operator.Dxx_mask[1, :, :], weights['dxx']),
-            'dyy': (Operator.Dyy_mask[1, :, :], weights['dyy']),
+            "dxy": (Operator.Dxy_mask[1, :, :], weights["dxy"] * 2),
+            "dxx": (Operator.Dxx_mask[1, :, :], weights["dxx"]),
+            "dyy": (Operator.Dyy_mask[1, :, :], weights["dyy"]),
         }
         return operators