In this first example, we will show how to construct a first basic model and the main objects and functions. First we import gempy:
# These two lines are necessary only if gempy is not installed
import sys, os
sys.path.append("../")
# Importing gempy
import gempy as gp
# Embedding matplotlib figures into the notebooks
%matplotlib inline
# Aux imports
import numpy as np
sys.path.append("../gempy")
import Visualization as vv
import importlib
importlib.reload(vv)
<module 'Visualization' from '../gempy/Visualization.py'>
All data get stored in a python object InputData. Therefore we can use python serialization to save the input of the models. In the next chapter we will see different ways to create data but for this example we will use a stored one
geo_data = gp.read_pickle('../input_data/NoFault.pickle')
geo_data.n_faults = 0
print(geo_data)
<gempy.DataManagement.InputData object at 0x7f5081146898>
This geo_data object contains essential information that we can access through the correspondent getters. Such a the coordinates of the grid.
#w.set_foliations()
#w.set_interfaces()
#
interp_data = gp.InterpolatorInput(geo_data, u_grade = [3],
)
print(interp_data)
[2, 2] <gempy.DataManagement.InterpolatorInput object at 0x7f5021b6aa20>
sol = gp.compute_model(interp_data)
[3]
gp.get_surface(sol[1,:], interp_data, 2)
(array([[ 0. , 130.44838905, -40. ], [ 40. , 160. , -54.08843994], [ 40. , 133.32850456, -40. ], ..., [ 1960. , 1680.08468628, -120. ], [ 1960. , 1693.94668579, -80. ], [ 1960. , 1706.76925659, -40. ]]), array([[ 2, 1, 0], [ 0, 1, 3], [ 1, 4, 3], ..., [3869, 3870, 3784], [3785, 3871, 3787], [3870, 3871, 3785]], dtype=int32))
getattr(interp_data, 'potential_at_interfaces')
False
fo = interp_data.get_formation_number().keys()
f_l = list(fo)[:-1]
v_l = []
s_l = []
for i in range(1,5):
v, s = gp.get_surface(sol[1, :], interp_data, i)
v_l.append(v)
s_l.append(s)
importlib.reload(vv)
w = vv.vtkVisualization(geo_data)
w.set_surfaces(v_l, s_l, f_l)
w.set_interfaces()
w.set_foliations()
w.render_model()
v_l, s_l
([array([[ 2.91980886e+01, 0.00000000e+00, 7.80060000e+06], [ 4.00000000e+01, 3.40689063e+00, 7.80060000e+06], [ 4.00000000e+01, 0.00000000e+00, 7.80059108e+06], ..., [ 1.96000000e+03, 1.92000000e+03, 7.80103210e+06], [ 1.96000000e+03, 1.93012238e+03, 7.80104000e+06], [ 1.96000000e+03, 1.96000000e+03, 7.80106450e+06]]), array([[ 0.00000000e+00, 1.30448389e+02, 7.80196000e+06], [ 4.00000000e+01, 1.60000000e+02, 7.80194591e+06], [ 4.00000000e+01, 1.33328505e+02, 7.80196000e+06], ..., [ 1.96000000e+03, 1.68008469e+03, 7.80188000e+06], [ 1.96000000e+03, 1.69394669e+03, 7.80192000e+06], [ 1.96000000e+03, 1.70676926e+03, 7.80196000e+06]]), array([[ 1.73798883e+01, 0.00000000e+00, 7.80184000e+06], [ 4.00000000e+01, 4.00000000e+01, 7.80182455e+06], [ 4.00000000e+01, 6.35708451e+00, 7.80184000e+06], ..., [ 1.96000000e+03, 1.80416306e+03, 7.80188000e+06], [ 1.96000000e+03, 1.81917297e+03, 7.80192000e+06], [ 1.96000000e+03, 1.83327148e+03, 7.80196000e+06]]), array([[ 0.00000000e+00, 4.39861526e+02, 7.80000000e+06], [ 0.00000000e+00, 4.37541199e+02, 7.80004000e+06], [ 9.82853532e+00, 4.40000000e+02, 7.80004000e+06], ..., [ 1.96000000e+03, 1.92000000e+03, 7.80044187e+06], [ 1.95097244e+03, 1.96000000e+03, 7.80044000e+06], [ 1.96000000e+03, 1.96000000e+03, 7.80044593e+06]])], [array([[ 2, 1, 0], [ 0, 1, 3], [ 1, 4, 3], ..., [4410, 4489, 4411], [4412, 4411, 4489], [4412, 4489, 4490]], dtype=int32), array([[ 2, 1, 0], [ 0, 1, 3], [ 1, 4, 3], ..., [3869, 3870, 3784], [3785, 3871, 3787], [3870, 3871, 3785]], dtype=int32), array([[ 2, 1, 0], [ 1, 3, 0], [ 1, 4, 3], ..., [4489, 4490, 4394], [4395, 4491, 4397], [4490, 4491, 4395]], dtype=int32), array([[ 2, 1, 0], [ 2, 0, 3], [ 5, 4, 1], ..., [3608, 3555, 3553], [3608, 3609, 3611], [3608, 3611, 3610]], dtype=int32)])
interp_data.extent
array([0.2052282051282051, 0.7180487179487179, 0.2436897435897436, 0.7565102564102564, 0.2629205128205128, 0.7757410256410256], dtype=object)
w.create_surface(v,s,'Seal')
(vtkRenderingOpenGL2Python.vtkOpenGLActor)0x7fcfbdf02948
import vtk
Points = vtk.vtkPoints()
vtk.
Points.InsertPoints?
This solution can be plot with the correspondent plotting function. Blocks:
gp.plot_section(geo_data, sol[0], 25)
<gempy.Visualization.PlotData at 0x7f629ef6b9e8>
Potential field:
gp.plot_potential_field(geo_data, sol[1], 25)