Thickness and Depth Maps for geological Interpretations

[AlexanderJuestel]

Is your feature request related to a problem? Please describe. For a geologist, it is vital to describe the depth of a layer in the underground. This could be relevant for geothermal exploration (depth of the target layers) or in general to describe the setting of a layer. In addition, the thickness of a layer may provide insights about local depocenters or just how thick the layer is across the whole model area (where does my reservoir have the highest thickness).

Describe the solution you'd like It would be great to have the possibility to plot the surface and its color as a function of the depth or thickness like commercial software is capable of. This could go in hand with a better method to select single layers for plotting next plot_all_surfaces.

Describe alternatives you've considered None but it should be fairly simple if one can access the different meshes and their attributes...

Add Option to display topography as well

[AlexanderJuestel]

A first depth map could be plotted by using the following code. This could be nicely implemented as a standalone function.

sargs = dict(fmt="%.0f",
            color='black')

gpv = vista.GemPyToVista(modelA5, extent = modelA5.grid.regular_grid.extent, plotter_type='basic', notebook = True)
gpv.p.camera_position =[(32302161.6978571, 5607748.0884169, 29527.965933430853),
 (32307438.614038907, 5629764.632521507, -2521.139653086674),
 (0.042944366888203855, 0.8201741150031122, 0.5704999583097559)]
surfaces_df = gpv._select_surfaces_data(modelA5.surfaces.df, surfaces=  ['ObererKohlenkalkGP'])
select_active = surfaces_df['isActive']
for idx, val in surfaces_df[select_active][['vertices', 'edges', 'color', 'surface', 'id']].dropna().iterrows():
    surf = pv.PolyData(val['vertices'], np.insert(val['edges'], 0, 3, axis=1).ravel())
    gpv.surface_poly[val['surface']] = surf
    array = surfaces_df['vertices'][16][:,2]   cgre-aachen/gempy#16 is the index of the surface in the df
    gpv.surface_actors[val['surface']] = gpv.p.add_mesh(
        surf,scalars = array, show_scalar_bar=True, cmap='gist_earth', clim = [-3000, 500], scalar_bar_args=sargs,stitle="Altitude [m]", smooth_shading=True)

    contours = surf.contour()
    gpv.p.add_mesh(contours, color="white", line_width=1)

gpv.p.add_point_labels(poly_elements, "Labels", point_size=10, font_size=10)

gpv.p.show_grid(color='black')
gpv.p.show()

[AlexanderJuestel]

To the more PyVista affine people among us (@Leguark, @alex-schaaf @banesullivan):

I have a surface (shown below) with vertices and edges and would like to create a mesh that only shows the top part (synclinal structure). Since only the top part is relevant, I would like to filter the data and remove the data below a certaint depth. However, removing the data and plotting it again results in a messed up mesh.

I tried filtering values below -750 m by using:

surfaces_df['vertices'][11] = surfaces_df['vertices'][11][np.where(surfaces_df['vertices'][11][:,2]>-750)]
surfaces_df['edges'][11] = surfaces_df['edges'][11][np.where(surfaces_df['vertices'][11][:,2]>-750)]

and applying it to the following code adapted from GemPy:

surfaces_df = gpv._select_surfaces_data(modelA5.surfaces.df, surfaces=  ['BreitgangFM'])

surfaces_df['vertices'][11] = surfaces_df['vertices'][11][np.where(surfaces_df['vertices'][11][:,2]>-500)]
surfaces_df['edges'][11] = surfaces_df['edges'][11][np.where(surfaces_df['vertices'][11][:,2]>-500)]

select_active = surfaces_df['isActive']
for idx, val in surfaces_df[select_active][['vertices', 'edges', 'color', 'surface', 'id']].dropna().iterrows():

    surf = pv.PolyData(val['vertices'], np.insert(val['edges'], 0, 3, axis=1).ravel())
    gpv.surface_poly[val['surface']] = surf
    array = surfaces_df['vertices'][11][:,2]
    gpv.surface_actors[val['surface']] = gpv.p.add_mesh(
        surf,scalars = array, show_scalar_bar=True, cmap='gist_earth', clim = [-1000, 500], scalar_bar_args=sargs,stitle="Altitude [m]", smooth_shading=True)

    contours = surf.contour()
    gpv.p.add_mesh(contours, color="white", line_width=1)

gpv.p.add_point_labels(poly_elements, "Labels", point_size=10, font_size=10)

gpv.p.show_grid(color='black')
gpv.p.show()

[banesullivan]

I have a surface (shown below) with vertices and edges and would like to create a mesh that only shows the top part (synclinal structure). Since only the top part is relevant, I would like to filter the data and remove the data below a certaint depth. However, removing the data and plotting it again results in a messed up mesh.

@AlexanderJuestel, I think I can help here but could you please simplify the problem down as much as possible? (i.e. give me the mesh in a VTK file)

Have you tried using a PyVista filter to do the extraction of the top surface rather than several np.wheres?

Also, is the mesh messed up, or are the colormap limits just not right?

---

Also, when I first read this, I immediately thought of https://github.com/pyvista/pyvista-support/issues/80 - perhaps this might be useful to you

[AlexanderJuestel]

Hello @banesullivan,

thanks for getting back to me! I have not tried filters before but will do that.

The colormap limits were adapted so that the synclinal structure in the centre is colored properly. The black part is just not needed. But the edges look a little off. I will have a look at that again.

Also thanks for mentioning the pyvista issue, I will certainly have a look and get back to you then!

[AlexanderJuestel]

Will be moved to GemGIS and integrated as part of the post processing

[bluetyson]

Depth definitely relevant for drilling!

[AlexanderJuestel]

Agree, as you can see I figured out how to create the depth maps already. However, vertices and edges are not returned in the current versions of GemPy 2.2.3 for some reason (https://github.com/cgre-aachen/gempy/issues/494).

As soon as that is fixed, we should be able to quickly implement it and apply it to the current example models ;)

[AlexanderJuestel]

cgre-aachen/gempy#494 was fixed by updating skimage to version 0.17.2

[AlexanderJuestel]

Closing this as depth and thickness maps seem to work now :)