gempy-project / gempy

GemPy is an open-source, Python-based 3-D structural geological modeling software, which allows the implicit (i.e. automatic) creation of complex geological models from interface and orientation data. It also offers support for stochastic modeling to address parameter and model uncertainties.
https://gempy.org
European Union Public License 1.2
986 stars 235 forks source link

Export 3D-Model as XZY file #922

Closed Sebadita closed 3 months ago

Sebadita commented 4 months ago

Hi everyone!

How can I export an interpolated mesh surface as XZY. file? I have created a model, which consist only of a 3D surface, and i want to use in other software.

yangjiandendi commented 4 months ago

Hi Sebadita,

firstly, you have to export your gempy model as VTK format:

p.surface_poly['rock1'].save('rock1.vtk') # Save the vtk file for formation 1
p.surface_poly['rock2'].save('rock2.vtk') # Save the vtk file for formation 2
p.orientations_mesh.save('orientations.vtk') # Save the vtk file for the orientations
p.surface_points_mesh.save('surface_points.vtk') # Save the vtk file for the surface points
box = p.regular_grid_actor.GetMapper().GetInput() # Get the vtk file for the regular grid
box.save('box.vtk')

then, you can get the vertices of mesh surface in VTK like this:

import vtk
import pandas as pd`
def generate_normals(polydata):
    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):

    surface_mesh = mesh.extract_surface()
    polydata = surface_mesh

    # Generate normals if not present
    polydata_with_normals = generate_normals(polydata)

    # Get points (vertices)
    points = polydata_with_normals.GetPoints()
    vertices = []
    for i in range(points.GetNumberOfPoints()):
        vertices.append(points.GetPoint(i))

    # Get normals
    normals_array = polydata_with_normals.GetPointData().GetNormals()
    normals = []
    for i in range(normals_array.GetNumberOfTuples()):
        normals.append(normals_array.GetTuple(i))

    return vertices, normals

def save_to_excel(vertices, normals, vertices_file, normals_file):
    # Create DataFrames
    vertices_df = pd.DataFrame(vertices, columns=['X', 'Y', 'Z'])
    normals_df = pd.DataFrame(normals, columns=['x', 'y', 'z'])

    # Save to Excel files
    vertices_df.to_excel(vertices_file, index=False)
    normals_df.to_excel(normals_file, index=False)

# mesh = pv.read("../data/model.vtp")
mesh = p.surface_poly['rock1']
vertices, normals = get_vertices_and_normals(mesh)
vertices_df = pd.DataFrame(vertices, columns=['X', 'Y', 'Z'])
normals_df = pd.DataFrame(normals, columns=['x', 'y', 'z'])
# Save to Excel files
vertices_file = "rock1_vertices.xlsx"
normals_file = "rock1_norms.xlsx"
save_to_excel(vertices, normals, vertices_file, normals_file)

finally you can convert pandas dataframe to an XYZ file:

# Convert the DataFrame to an XYZ file
def dataframe_to_xyz(df, filename):
    with open(filename, 'w') as f:
        for index, row in df.iterrows():
            f.write(f"{row['X']} {row['Y']} {row['Z']}\n")

# Specify the filename
filename = "output.xyz"

# Call the function to write the DataFrame to an XYZ file
dataframe_to_xyz(vertices_df, filename)
AlexanderJuestel commented 4 months ago

@Sebadita this issue relates to #898

please also see this code that should allow you to export the meshes. surfaces_poly is the dict where the mesh is stored. You can access the mesh with surfaces_poly['Layer1'][0]. The XYZ coordinates of the vertices can be accessed with surfaces_poly['Layer1'][0].points. This will return irregularly spaced points. Do you need regular-spaced points like for a raster?

Here the code to export the meshes with their true XYZ coordinates


# List of Surfaces
surfaces = ['Layer1', 'Layer2']

# Getting a list of all surfaces
list_surfaces = list(geo_model.structural_frame.element_name_id_map.keys())

# Getting indices of provided surfaces
list_indices = [list_surfaces.index(surface) for surface in surfaces]

# Creating empty dict to store data
surfaces_poly = {}

for index in list_indices:

    # Extracting vertices
    vertices = geo_model.input_transform.apply_inverse(geo_model.solutions.raw_arrays.vertices[index])

    # Extracting faces
    faces = np.insert(geo_model.solutions.raw_arrays.edges[index], 0, 3, axis=1).ravel()

    # Creating PolyData from vertices and faces
    surf = pv.PolyData(vertices, faces)

    # Appending depth to PolyData
    surf['Depth [m]'] = geo_model.input_transform.apply_inverse(geo_model.solutions.raw_arrays.vertices[index])[:, 2]

    # Storing mesh, depth values and color values in dict
    surfaces_poly[list_surfaces[index]] = [surf, geo_model.structural_frame.elements_colors[index]]
javoha commented 3 months ago

Answered. Please feel free to reopen if you have additional questions.