Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Export 3D-Model as XZY file #922

Closed
Sebadita opened this issue Jul 16, 2024 · 3 comments
Closed

Export 3D-Model as XZY file #922

Sebadita opened this issue Jul 16, 2024 · 3 comments

Comments

@Sebadita
Copy link

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
Copy link

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
Copy link
Contributor

@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
Copy link
Member

javoha commented Aug 12, 2024

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

@javoha javoha closed this as completed Aug 12, 2024
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

4 participants