Merge branch 'main' into fork/flohorovicic/feature/json_io

Author: Leguark

.github/PULL_REQUEST_TEMPLATE.md

@@ -4,12 +4,9 @@ Please include a summary of the changes.
 Relates to <issue>
 
 # Checklist
-- [ ] My code follows the [PEP 8 style guidelines](https://www.python.org/dev/peps/pep-0008/).
 - [ ] My code uses type hinting for function and method arguments and return values.
-- [ ] My code contains descriptive and helpful docstrings 
-which are formatted per the [Google Python Style Guidelines](http://google.github.io/styleguide/pyguide.html).
-- [ ] I have created tests which entirely cover my code.
+- [ ] I have created tests which cover my code.
 - [ ] The test code either 1. demonstrates at least one valuable use case (e.g. integration tests) 
 or 2. verifies that outputs are as expected for given inputs (e.g. unit tests).
-- [ ] New and existing tests pass locally with my changes.
+- [ ] New tests pass locally with my changes.
 

README.md

@@ -2,10 +2,11 @@
 
 [![GitHub Stars](https://img.shields.io/github/stars/cgre-aachen/gempy.svg)](https://github.com/cgre-aachen/gempy/stargazers)
 [![GitHub Forks](https://img.shields.io/github/forks/cgre-aachen/gempy.svg)](https://github.com/cgre-aachen/gempy/network)
+[![Documentation Status](https://assets.readthedocs.org/static/projects/badges/passing-flat.svg)](http://docs.gempy.org)
+[![Build Status](http://terranigma-solutions.teamcity.com/app/rest/builds/buildType:(id:Gempy_TestingGemPy)/statusIcon)](http://terranigma-solutions.teamcity.com/viewType.html?buildTypeId=Gempy_TestingGemPy&guest=1)
 [![PyPI](https://img.shields.io/badge/python-3.10%20%7C%203.11%20%7C%203.12-blue)](https://www.python.org/downloads/)
 [![PyPI](https://img.shields.io/badge/pypi-1.0-blue.svg)](https://pypi.org/project/gempy/)
 [![license: EUPL v1.2](https://img.shields.io/badge/license-EUPL%20v1.2-blue.svg)](https://github.com/cgre-aachen/gempy/blob/master/LICENSE)
-[![Documentation Status](https://assets.readthedocs.org/static/projects/badges/passing-flat.svg)](http://docs.gempy.org)
 [![DOI](https://zenodo.org/badge/96211155.svg)](https://zenodo.org/badge/latestdoi/96211155)
 
 

examples/examples/geometries/a01_horizontal_stratigraphic.py

@@ -1,27 +1,39 @@
 """
-Model 1 - Horizontal stratigraphy
-==================================
+Model 1 - Horizontal Stratigraphy
+===================================
 
-A simple example with horizontal layers
+This example demonstrates how to build a basic geological model with horizontally stacked layers using GemPy,
+an open-source Python library for implicit geological modeling. In this script, we will:
 
-This script demonstrates how to create a basic model with horizontally stacked layers using GemPy,
-a Python-based, open-source library for implicit geological modeling.
-"""
+- Set up and compute a simple geological model using input data (orientations and surface points).
+- Visualize the model in 2D and 3D.
+- Export model components (formations, orientations, surface points, and the regular grid) to VTK files.
+- Extract mesh vertices and normals, save them to Excel files, and convert the vertex data to an XYZ file.
 
+Each section includes detailed comments to explain its purpose and functionality.
+"""
 
 # %%
 # Import necessary libraries
 import gempy as gp
 import gempy_viewer as gpv
-
+import pyvista as pv
+import vtk
+import pandas as pd
 
 # sphinx_gallery_thumbnail_number = 2
 
 # %%
-# Generate the model
-# Define the path to data
+# Generate the Geological Model
+# -----------------------------
+# In this section, we define the data source and create a GeoModel instance.
+# The model is defined with a specific spatial extent and resolution (refinement).
+# Input data (orientations and surface points) are loaded from a remote repository.
+
+# Define the path to the dataset hosted online
 data_path = 'https://raw.githubusercontent.com/cgre-aachen/gempy_data/master/'
-# Create a GeoModel instance
+
+# Create a GeoModel instance with the given project name, extent, and refinement level.
 data = gp.create_geomodel(
     project_name='horizontal',
     extent=[0, 1000, 0, 1000, 0, 1000],
@@ -31,19 +43,174 @@
         path_to_surface_points=data_path + "/data/input_data/jan_models/model1_surface_points.csv"
     )
 )
-# Map geological series to surfaces
+
+# Map a geological series to the corresponding surfaces.
+# Here, the "Strat_Series" is associated with two formations ('rock2' and 'rock1'),
+# which defines the stacking order of horizontal layers.
 gp.map_stack_to_surfaces(
     gempy_model=data,
     mapping_object={"Strat_Series": ('rock2', 'rock1')}
 )
-# Compute the geological model
+
+# Compute the geological model using the provided data and mappings.
 gp.compute_model(data)
 geo_data = data
 
 # %%
-# Plot the initial geological model in the y direction without results
+# 2D Visualization of the Geological Model
+# -----------------------------------------
+# The following plots provide 2D views of the model.
+# - The first plot shows the initial model without computed results.
+# - The subsequent plots display the computed geological data in the x and y directions.
+# Note: The boundaries are hidden in the latter plots for a cleaner visualization.
+
+# Plot the initial model in the y direction (without displaying computed results)
 gpv.plot_2d(geo_data, direction=['y'], show_results=False)
 
-# Plot the result of the model in the x and y direction with data and without boundaries
+# Plot the computed model results in the x and y directions, including data but excluding boundaries
 gpv.plot_2d(geo_data, direction=['x'], show_data=True, show_boundaries=False)
 gpv.plot_2d(geo_data, direction=['y'], show_data=True, show_boundaries=False)
+
+# %%
+# 3D Visualization and Export to VTK
+# ----------------------------------
+# Visualize the model in 3D with data, results, and boundaries.
+# The generated 3D plot is also used to export various components of the model to VTK files.
+p = gpv.plot_3d(geo_data, show_data=True, show_results=True, show_boundaries=True, image=True)
+
+# Export different components of the model to VTK files:
+p.surface_poly['rock1'].save('rock1.vtk')  # Save formation 'rock1'
+p.surface_poly['rock2'].save('rock2.vtk')  # Save formation 'rock2'
+p.orientations_mesh.save('orientations.vtk')  # Save the orientations mesh
+p.surface_points_mesh.save('surface_points.vtk')  # Save the surface points mesh
+
+# Retrieve and export the regular grid (volume representation) of the model
+box = p.regular_grid_actor.GetMapper().GetInput()
+box.save('box.vtk')
+
+# %%
+# Display the Exported VTK Files with PyVista
+# -------------------------------------------
+# Load and plot each exported VTK file to verify the export.
+pv.read('rock1.vtk').plot(show_edges=False)
+pv.read('rock2.vtk').plot(show_edges=False)
+pv.read('orientations.vtk').plot(show_edges=False)
+pv.read('surface_points.vtk').plot(show_edges=False)
+pv.read('box.vtk').plot(show_edges=False)
+
+
+# %%
+# Extract and Save Mesh Vertices and Normals
+# ------------------------------------------
+# The following functions extract vertex coordinates and corresponding normal vectors
+# from a mesh (vtkPolyData), and then save this data into Excel files for further analysis.
+
+def generate_normals(polydata):
+    """
+    Generate normals for the given polydata if they are not already computed.
+
+    Parameters:
+        polydata (vtk.vtkPolyData): Input polydata for which normals are to be computed.
+
+    Returns:
+        vtk.vtkPolyData: The polydata with computed point normals.
+    """
+    normal_generator = vtk.vtkPolyDataNormals()
+    normal_generator.SetInputData(polydata)
+    normal_generator.ComputePointNormalsOn()
+    normal_generator.ComputeCellNormalsOff()
+    normal_generator.Update()
+    return normal_generator.GetOutput()
+
+
+def get_vertices_and_normals(mesh):
+    """
+    Extract vertices and normals from a mesh.
+
+    Parameters:
+        mesh: The mesh object (typically from a formation) to extract surface vertices and normals.
+
+    Returns:
+        tuple: Two lists containing vertices and normals respectively.
+    """
+    # Extract the surface of the mesh
+    surface_mesh = mesh.extract_surface()
+    polydata = surface_mesh
+
+    # Ensure normals are computed for the polydata
+    polydata_with_normals = generate_normals(polydata)
+
+    # Extract points (vertices)
+    points = polydata_with_normals.GetPoints()
+    vertices = [points.GetPoint(i) for i in range(points.GetNumberOfPoints())]
+
+    # Extract normals associated with the points
+    normals_array = polydata_with_normals.GetPointData().GetNormals()
+    normals = [normals_array.GetTuple(i) for i in range(normals_array.GetNumberOfTuples())]
+
+    return vertices, normals
+
+
+def save_to_excel(vertices, normals, vertices_file, normals_file):
+    """
+    Save the vertices and normals to separate Excel files.
+
+    Parameters:
+        vertices (list): List of vertex coordinates.
+        normals (list): List of normal vectors.
+        vertices_file (str): File name for saving vertices.
+        normals_file (str): File name for saving normals.
+    """
+    # Create DataFrames from the lists
+    vertices_df = pd.DataFrame(vertices, columns=['X', 'Y', 'Z'])
+    normals_df = pd.DataFrame(normals, columns=['x', 'y', 'z'])
+
+    # Write the DataFrames to Excel files
+    vertices_df.to_excel(vertices_file, index=False)
+    normals_df.to_excel(normals_file, index=False)
+
+
+# Extract vertices and normals from the mesh of formation 'rock1'
+mesh = p.surface_poly['rock1']
+vertices, normals = get_vertices_and_normals(mesh)
+
+# Create DataFrames for later processing or verification
+vertices_df = pd.DataFrame(vertices, columns=['X', 'Y', 'Z'])
+normals_df = pd.DataFrame(normals, columns=['x', 'y', 'z'])
+
+# Define file names for the Excel outputs
+vertices_file = "rock1_vertices.xlsx"
+normals_file = "rock1_norms.xlsx"
+
+# Save the extracted data to Excel files
+save_to_excel(vertices, normals, vertices_file, normals_file)
+
+# Optionally, read back the Excel files to verify their content
+pd.read_excel(vertices_file)
+pd.read_excel(normals_file)
+
+
+# %%
+# Convert Vertices DataFrame to an XYZ File
+# -----------------------------------------
+# This function converts a DataFrame containing vertex coordinates to an XYZ format file,
+# which is a simple text file with one point per line.
+
+def dataframe_to_xyz(df, filename):
+    """
+    Write vertex coordinates from a DataFrame to an XYZ file.
+
+    Parameters:
+        df (pandas.DataFrame): DataFrame containing 'X', 'Y', and 'Z' columns.
+        filename (str): Output filename for the XYZ file.
+    """
+    with open(filename, 'w') as f:
+        for index, row in df.iterrows():
+            f.write(f"{row['X']} {row['Y']} {row['Z']}\n")
+
+
+# Specify the output file name for the XYZ file
+filename = "output.xyz"
+
+# Convert the vertices DataFrame to an XYZ file
+dataframe_to_xyz(vertices_df, filename)

examples/examples/real/Perth_basin.py

@@ -25,8 +25,7 @@
 geo_model: gp.data.GeoModel = gp.create_geomodel(
     project_name='Perth_Basin',
     extent=[337000, 400000, 6640000, 6710000, -18000, 1000],
-    resolution=[100, 100, 100],
-    refinement=4,
+    refinement=6,
     importer_helper=gp.data.ImporterHelper(
         path_to_orientations=data_path + "/data/input_data/perth_basin/Paper_GU2F_sc_faults_topo_Foliations.csv",
         path_to_surface_points=data_path + "/data/input_data/perth_basin/Paper_GU2F_sc_faults_topo_Points.csv",

gempy/core/data/geo_model.py

@@ -147,6 +147,8 @@ def solutions(self) -> Solutions:
 
     @solutions.setter
     def solutions(self, value):
+        # * This is set  from the gempy engine
+        
         self._solutions = value
 
         # * Set solutions per group

gempy/core/data/grid.py

@@ -254,7 +254,11 @@ def _update_values(self):
         if self.GridTypes.CENTERED in self.active_grids:
             if self.centered_grid is None: raise AttributeError('Centered grid is active but not defined')
             values.append(self.centered_grid.values)
-
+        
+        # make sure values is not empty
+        if len(values) == 0:
+            return self.values
+        
         self.values = np.concatenate(values)
 
         return self.values

gempy/modules/mesh_extranction/__init__.py

@@ -0,0 +1,88 @@
+import numpy as np
+from typing import Optional
+from skimage import measure
+from gempy_engine.core.data.interp_output import InterpOutput
+from gempy_engine.core.data.raw_arrays_solution import RawArraysSolution
+
+from gempy.core.data import GeoModel, StructuralElement, StructuralGroup
+from gempy.core.data.grid_modules import RegularGrid
+
+
+def set_meshes_with_marching_cubes(model: GeoModel) -> None:
+    """Extract meshes for all structural elements using the marching cubes algorithm.
+    
+    Parameters
+    ----------
+    model : GeoModel
+        The geological model containing solutions and structural elements.
+    
+    Raises
+    ------
+    ValueError
+        If the model solutions do not contain dense grid data.
+    """
+    # Verify that solutions contain dense grid data
+    if (model.solutions is None or 
+            model.solutions.block_solution_type != RawArraysSolution.BlockSolutionType.DENSE_GRID):
+        raise ValueError("Model solutions must contain dense grid data for mesh extraction.")
+    
+    regular_grid: RegularGrid = model.grid.regular_grid
+    structural_groups: list[StructuralGroup] = model.structural_frame.structural_groups
+    
+    if not model.solutions.octrees_output or not model.solutions.octrees_output[0].outputs_centers:
+        raise ValueError("No interpolation outputs available for mesh extraction.")
+    
+    output_lvl0: list[InterpOutput] = model.solutions.octrees_output[0].outputs_centers
+
+    for e, structural_group in enumerate(structural_groups):
+        if e >= len(output_lvl0):
+            continue
+            
+        output_group: InterpOutput = output_lvl0[e]
+        scalar_field_matrix = output_group.exported_fields_dense_grid.scalar_field
+        
+        for element in structural_group.elements:
+            extract_mesh_for_element(
+                structural_element=element,
+                regular_grid=regular_grid,
+                scalar_field=scalar_field_matrix
+            )
+
+
+def extract_mesh_for_element(structural_element: StructuralElement, 
+                             regular_grid: RegularGrid, 
+                             scalar_field: np.ndarray, 
+                             mask: Optional[np.ndarray] = None) -> None:
+    """Extract a mesh for a single structural element using marching cubes.
+    
+    Parameters
+    ----------
+    structural_element : StructuralElement
+        The structural element for which to extract a mesh.
+    regular_grid : RegularGrid
+        The regular grid defining the spatial discretization.
+    scalar_field : np.ndarray
+        The scalar field used for isosurface extraction.
+    mask : np.ndarray, optional
+        Optional mask to restrict the mesh extraction to specific regions.
+    """
+    # Apply mask if provided
+    volume = scalar_field.reshape(regular_grid.resolution)
+    if mask is not None:
+        volume = volume * mask
+    
+    # Extract mesh using marching cubes
+    verts, faces, _, _ = measure.marching_cubes(
+        volume=volume,
+        level=structural_element.scalar_field,
+        spacing=(regular_grid.dx, regular_grid.dy, regular_grid.dz)
+    )
+    
+    # Adjust vertices to correct coordinates in the model's extent
+    verts = (verts + [regular_grid.extent[0],
+                      regular_grid.extent[2],
+                      regular_grid.extent[4]])
+
+    # Store mesh in the structural element
+    structural_element.vertices = verts
+    structural_element.edges = faces

gempy/modules/mesh_extranction/marching_cubes.py

@@ -1,7 +1,8 @@
 -r optional-requirements.txt
 
 pytest
-subsurface~=2024.2.0
+pyvista
+git+https://github.com/terranigma-solutions/[email protected]
 
 # Testing
 pytest-approvaltests