How to display a surface as a cube about gemgis HOT 6 CLOSED

@ada-hou thanks for liking GemGIS so much.

What do you mean by "as one cube" exactly? What you can do is you can represent the result of the GemPy modeling using the voxels of the GemPy lith_block. The surfaces are actually computed based on this voxel representation. So, after computing your model, geo_model.lith_block should provide you with a flattened array where the numbers are equal to the IDs of the layers. You then just have to reshape it according to your resolution.

from gemgis.

Thank you for your reply, just like the image below, each geological layer is not just a single surface.

from gemgis.

As you are using Blender, I am pinging our Blender expert @NilsChudalla to help out :)

from gemgis.

The model above was generated using a different method, while the model below was created using GemPy, GemGIS, and PyVista

This is my code
#%%

import gempy as gp
import numpy as np
import matplotlib.pyplot as plt

geo_model = gp.create_model('Model1')
geo_model = gp.init_data(geo_model, extent=[0, 791, 0, 200, -582, 0], resolution=[100, 10, 100])
print(geo_model)

#%%

gp.set_interpolator(geo_model, aesara_optimizer='fast_compile', verbose=[])
geo_model.set_default_surfaces()
geo_model.add_surface_points(X=223, Y=0.01, Z=-94, surface='surface1')
geo_model.add_surface_points(X=458, Y=0, Z=-107, surface='surface1')
geo_model.add_surface_points(X=612, Y=0, Z=-14, surface='surface1')
geo_model.add_orientations(X=350, Y=0, Z=300, surface='surface1', pole_vector=(0, 0, 1))
geo_model.add_surface_points(X=225, Y=0, Z=-269, surface='surface2')
geo_model.add_surface_points(X=459, Y=0, Z=-279, surface='surface2')
geo_model.add_surfaces(['surface3', 'basement'])
geo_model.add_surface_points(X=225, Y=0, Z=-439, surface='surface3')
geo_model.add_surface_points(X=464, Y=0, Z=-456, surface='surface3')
geo_model.add_surface_points(X=619, Y=0, Z=-433, surface='surface3')

#%%

geo_model.add_features('Fault1')
geo_model.reorder_features(['Fault1', 'Default series'])
geo_model.set_is_fault('Fault1')
geo_model.add_surfaces('fault1')
gp.map_stack_to_surfaces(geo_model, {'Fault1': 'fault1'})
geo_model.add_surface_points(X=550, Y=0, Z=-30, surface='fault1')
geo_model.add_surface_points(X=650, Y=0, Z=-200, surface='fault1')
geo_model.add_orientations(X=600, Y=0, Z=-100, surface='fault1', pole_vector=(.3, 0, .3))

#%%

gp.compute_model(geo_model)

#%%

gp.plot_3d(geo_model)

#%%

import gemgis as gg
surface1 = gg.visualization.create_depth_maps_from_gempy(geo_model=geo_model, surfaces='surface1')
surface2 = gg.visualization.create_depth_maps_from_gempy(geo_model=geo_model, surfaces='surface2')
surface3 = gg.visualization.create_depth_maps_from_gempy(geo_model=geo_model, surfaces='surface3')
fault1 = gg.visualization.create_depth_maps_from_gempy(geo_model=geo_model, surfaces='fault1')
basement = gg.visualization.create_depth_maps_from_gempy(geo_model=geo_model, surfaces='basement')
print(surface1)

#%%

import pyvista as pv
p = pv.Plotter(notebook=False)
p.add_mesh(surface1['surface1'][0])
p.add_mesh(surface2['surface2'][0])
p.add_mesh(surface3['surface3'][0])
p.add_mesh(fault1['fault1'][0])

#%%

p.export_gltf("haha.gltf")

#%%

geo_model.additional_data.kriging_data

#%%

from gemgis.

The voxels displayed in Gempy are not meshes, with edges and vertices. You can export as vtk and use Blender addons for vtk format, such as https://bvtknodes.readthedocs.io/en/latest/BVTKNodes.html. If it is onetime usage and you are not interested in a perfect watertight model because you are only looking for visualization purposes, try this maybe:

• extrude surfaces along z direction until they intersect the next surface
• extrude last layer until your desired bottom end (keep extruded surface selected)
• scale selected layer with x=1, y=1,z=0, resulting in a flat surface
• apply boolean filter to upper layer, with difference to lower layer

Be aware, that the boolean filter might be behaving a bit strange, if the mesh is either not closed (which it should be after extruding) or if some faces of you mesh are not assigned correctly. Therefore you might have to recompute normals in addition to the steps above.

Edit: especially reverse faults might make the "boolean" approach impossible if one surface is encountered at the same position horizontally in multiple depths.

from gemgis.

Thank you for your replay. It sounds quite complex, and I may need some time to verify it

from gemgis.