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
941 stars 232 forks source link

After the model is exported to pyvista, the coordinates are incorrect #898

Closed Lian-Po-s-Study-Diary closed 2 months ago

Lian-Po-s-Study-Diary commented 2 months ago

After I built the model using gempy3.0, I used the create_depth_maps_from_gempy method of gemgis to output the model to pyvista and found that the coordinates were wrong. It may be because gempy uses the method of normalizing coordinates. How should I restore the model to the actual size and actual position? Or, is there a way to directly output model files in gempy3.0?

javoha commented 2 months ago

Hi, so if understand you correctly the coordinates within gempy are correct, but using the gemgis function does not work? As far as I know gemgis has not been updated to work with gempy v3, so it might take a while until this functionality is available. I would suggest creating an issue in gemgis to ask about the current state. Gempy v3 as of now does not have the functionality to save a full model. You can however export all the raw arrays of the result using .solutions. Just some examples: geo_model.solutions.dc_meshes to retrieve the meshes plotted in vista, or data.grid.regular_grid.values and data.solutions.raw_arrays.lith_block for getting the vertices and lithologies of the grid. Does this help?

Lian-Po-s-Study-Diary commented 2 months ago

Thank you for your help. I successfully restored the model using geo_model. transforme

AlexanderJuestel commented 2 days ago

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]]