Chapter 1: GemPy Basic

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

All data get stored in a python object InputData. This object can be easily stored in a Python pickle. However, these files have the limitation that all dependecies must have the same versions as those when the pickle were created. For these reason to have more stable tutorials we will generate the InputData from raw data—i.e. csv files exported from Geomodeller.

These csv files can be found in the input_data folder in the root folder of GemPy. These tables contains uniquely the XYZ (and poles, azimuth and polarity in the foliation case) as well as their respective formation name (but not necessary the formation order).

# Importing the data from csv files and settign extent and resolution
geo_data = gp.create_data([0,2000,0,2000,-2000,0],[ 50,50,50],
                         path_f = os.pardir+"/input_data/FabLessPoints_Foliations.csv",
                         path_i = os.pardir+"/input_data/FabLessPoints_Points.csv")
import pandas as pn
gp.set_interfaces(geo_data,
                  pn.DataFrame(data = {"X" : (1200, 1800),
                                       "Y" : (1000,  1000),
                                       "Z" : (- 1600, -1800),
                                       "formation" : np.tile("Layer 3", 2),
                                                      }), append=True)



gp.set_foliations(geo_data,
                   pn.DataFrame(data = {"X" : (1500),
                      "Y" : (1000 ),
                      "Z" : ( -1400),
                      "azimuth" : ( 90),
                      "dip" : (20),
                      "polarity" : (1),
                      "formation" : ['Layer 3'],
                                       }),append=True)


gp.set_interfaces(geo_data,
                  pn.DataFrame(data = {"X" : (500, 600),
                                       "Y" : (1000,  1000),
                                       "Z" : (- 150, -1700),
                                       "formation" : np.tile("fault2", 2),
                                                      }), append=True)



gp.set_foliations(geo_data,
                   pn.DataFrame(data = {"X" : (400),
                      "Y" : (1000),
                      "Z" : ( -1300),
                      "azimuth" : ( 90),
                      "dip" : (90),
                      "polarity" : (1),
                      "formation" : ['fault2'],
                                       }),append=True)

With the command get data is possible to see all the input data. However, at the moment the (depositional) order of the formation and the separation of the series (More explanation about this in the next notebook) is totally arbitrary.

gp.get_data(geo_data)
X Y Z azimuth dip formation polarity series
interfaces 0 800 1000 -1600 NaN NaN MainFault NaN fault
1 1200 1000 -400 NaN NaN MainFault NaN fault
2 1100 1000 -700 NaN NaN MainFault NaN fault
3 900 1000 -1300 NaN NaN MainFault NaN fault
4 1000 1000 -1000 NaN NaN MainFault NaN fault
5 800 200 -1400 NaN NaN Reservoir NaN Rest
6 800 1800 -1400 NaN NaN Reservoir NaN Rest
7 600 1000 -1050 NaN NaN Reservoir NaN Rest
8 300 1000 -950 NaN NaN Reservoir NaN Rest
9 2000 1000 -1275 NaN NaN Reservoir NaN Rest
10 1900 1000 -1300 NaN NaN Reservoir NaN Rest
11 1300 1000 -1100 NaN NaN Reservoir NaN Rest
12 1600 1000 -1200 NaN NaN Reservoir NaN Rest
13 1750 1000 -1250 NaN NaN Reservoir NaN Rest
14 1000 100 -1100 NaN NaN Reservoir NaN Rest
15 1000 1975 -1050 NaN NaN Reservoir NaN Rest
16 1000 1900 -1100 NaN NaN Reservoir NaN Rest
17 1150 1000 -1050 NaN NaN Reservoir NaN Rest
18 1000 25 -1050 NaN NaN Reservoir NaN Rest
19 1300 1000 -600 NaN NaN SecondaryReservoir NaN Rest
20 600 1000 -550 NaN NaN SecondaryReservoir NaN Rest
21 900 1000 -650 NaN NaN SecondaryReservoir NaN Rest
22 1600 1000 -700 NaN NaN SecondaryReservoir NaN Rest
23 1900 1000 -800 NaN NaN SecondaryReservoir NaN Rest
24 2000 1000 -775 NaN NaN SecondaryReservoir NaN Rest
25 2000 1000 -875 NaN NaN Seal NaN Rest
26 600 1000 -650 NaN NaN Seal NaN Rest
27 1300 1000 -700 NaN NaN Seal NaN Rest
28 1900 1000 -900 NaN NaN Seal NaN Rest
29 900 1000 -750 NaN NaN Seal NaN Rest
30 1600 1000 -800 NaN NaN Seal NaN Rest
31 600 1000 -1350 NaN NaN Overlying NaN Rest
32 300 1000 -1250 NaN NaN Overlying NaN Rest
33 1000 1000 -1300 NaN NaN Overlying NaN Rest
34 2000 1000 -1575 NaN NaN Overlying NaN Rest
35 1300 1000 -1400 NaN NaN Overlying NaN Rest
36 1600 1000 -1500 NaN NaN Overlying NaN Rest
37 1750 1000 -1550 NaN NaN Overlying NaN Rest
38 1900 1000 -1600 NaN NaN Overlying NaN Rest
39 1800 1000 200 NaN NaN Layer 3 NaN unco
40 1200 1000 200 NaN NaN Layer 3 NaN unco
foliations 0 917.45 1000 -1135.4 270 71.565 MainFault 1 fault
1 1450 1000 -1150 90 18.435 Reservoir 1 Rest
2 1500 1000 400 90 0 Layer 3 1 unco

To set the number of (depositional) series and which formation belongs to which, it is possible to use the function set_series. Here, there are two important things to notice:

  • set_series requires a dictionary. In Python dictionaries are not order (keys dictionaries since Python 3.6 are) and therefore there is not guarantee of having the right order.
    • The order of the series are vital for the method (from younger to older).
    • The order of the formations (as long they belong to the correct series) are only important for the color code. If the order of the pile differs from the final result the color of the interfaces and input data will be different
  • Every fault is treated as an independent series and have to be at the top of the pile. The relative order between the distinct faults represent the tectonic relation between them (from younger to older as well).

The order of the series (for Python < 3.6, otherwise passing the correct order of keys in the dictionary is enough) can be given as attribute as in the cell below. For the order of formations can be passed as attribute as well or using the specific function set_order_formations.

# Assigning series to formations as well as their order (timewise)
gp.set_series(geo_data, {"fault":'MainFault',
                      "Rest": ('SecondaryReservoir','Seal', 'Reservoir', 'Overlying'),
                        'unco' : 'Layer 3', 'Fault_serie2':'fault2'},
                       order_series = ["fault", 'Fault_serie2', 'Rest', 'unco'],
                       order_formations=['MainFault', 'fault2',
                                         'SecondaryReservoir', 'Seal','Reservoir', 'Overlying', 'Layer 3',
                                         ])
<gempy.strat_pile.StratigraphicPile at 0x7fdf34fda898>
../_images/ch1-Copy1_8_1.png

As an alternative the stratigraphic pile is interactive given the right backend (try %matplotlib notebook or %matplotlib qt5). These backends sometimes give some trouble though. Try to execute the cell twice:

%matplotlib notebook
gp.get_stratigraphic_pile(geo_data)
<IPython.core.display.Javascript object>
<gempy.strat_pile.StratigraphicPile at 0x7f6e21c1d470>

Notice that the colors depends on the order and therefore every time the cell is executed the colors are always in the same position. Be aware of the legend to be sure that the pile is as you wish!! (In the future every color will have the annotation within the rectangles to avoid confusion)

This geo_data object contains essential information that we can access through the correspondent getters. Such a the coordinates of the grid.

print(gp.get_grid(geo_data))
[[    0.             0.         -2000.        ]
 [    0.             0.         -1959.18371582]
 [    0.             0.         -1918.36730957]
 ...,
 [ 2000.          2000.           -81.63265228]
 [ 2000.          2000.           -40.81632614]
 [ 2000.          2000.             0.        ]]

The main input the potential field method is the coordinates of interfaces points as well as the orientations. These pandas dataframes can we access by the following methods:

Interfaces Dataframe

gp.get_data(geo_data, 'interfaces').head()
X Y Z formation series annotations
0 1000 1000 -1000 MainFault fault ${\bf{x}}_{\alpha \,{\bf{1}},0}$
1 800 1000 -1600 MainFault fault ${\bf{x}}_{\alpha \,{\bf{1}},1}$
2 1200 1000 -400 MainFault fault ${\bf{x}}_{\alpha \,{\bf{1}},2}$
3 1100 1000 -700 MainFault fault ${\bf{x}}_{\alpha \,{\bf{1}},3}$
4 900 1000 -1300 MainFault fault ${\bf{x}}_{\alpha \,{\bf{1}},4}$

Foliations Dataframe

Now the formations and the series are correctly set.

gp.get_data(geo_data, 'foliations').head()
X Y Z dip azimuth polarity formation series annotations
0 917.45 1000.0 -1135.398 71.565 270.0 1 MainFault fault ${\bf{x}}_{\beta \,{\bf{1}},0}$
1 1450.00 1000.0 -1150.000 18.435 90.0 1 Reservoir Rest ${\bf{x}}_{\beta \,{\bf{4}},0}$

It is important to notice the columns of each data frame. These not only contains the geometrical properties of the data but also the formation and series at which they belong. This division is fundamental in order to preserve the depositional ages of the setting to model.

A projection of the aforementioned data can be visualized in to 2D by the following function. It is possible to choose the direction of visualization as well as the series:

%matplotlib inline
gp.plot_data(geo_data, direction='y')
../_images/ch1-Copy1_20_0.png

GemPy supports visualization in 3D as well trough vtk. These plots are interactive. Try to drag and drop a point or interface! In the perpendicular views only 2D movements are possible to help to place the data where is required.

gp.plot_data_3D(geo_data)

The ins and outs of Input data objects

As we have seen objects DataManagement.InputData (usually called geo_data in the tutorials) aim to have all the original geological properties, measurements and geological relations stored.

Once we have the data ready to generate a model, we will need to create the next object type towards the final geological model:

interp_data = gp.InterpolatorInput(geo_data, u_grade=[3,3, 3, 3])#, verbose=['faults_matrix', 'covariance_matrix'])
print(interp_data)
Level of Optimization:  fast_run
Device:  cpu
Precision:  float32
<gempy.DataManagement.InterpolatorInput object at 0x7fdf34ec0320>
len(interp_data.interpolator.tg.len_series_f.get_value())-1
3
interp_data.get_formation_number()
{'DefaultBasement': 0,
 'Layer 3': 6,
 'MainFault': 1,
 'Overlying': 5,
 'Reservoir': 4,
 'Seal': 3,
 'SecondaryReservoir': 2}

By default (there is a flag in case you do not need) when we create a interp_data object we also compile the theano function that compute the model. That is the reason why takes long.

gempy.DataManagement.InterpolatorInput (usually called interp_data in the tutorials) prepares the original data to the interpolation algorithm by scaling the coordinates for better and adding all the mathematical parametrization needed.

gp.get_kriging_parameters(interp_data)
range 0.8882311582565308 3464.1015172
Number of drift equations [2 2]
Covariance at 0 0.01878463476896286
Foliations nugget effect 0.009999999776482582

These later parameters have a default value computed from the original data or can be changed by the user (be careful of changing any of these if you do not fully understand their meaning).

At this point, we have all what we need to compute our model. By default everytime we compute a model we obtain:

  • Lithology block model
    • with the lithology values in 0
    • with the potential field values in 1
  • Fault block model
    • with the faults zones values (i.e. every divided region by each fault has one number) in 0
    • with the potential field values in 1
interp_data.geo_data_res.interfaces
X Y Z annotations formation formation number isFault order_series series
0 0.410356 0.5001 0.371895 ${\bf{x}}_{\alpha \,{\bf{1}},0}$ MainFault 1 True 1 fault
1 0.512921 0.5001 0.679587 ${\bf{x}}_{\alpha \,{\bf{1}},1}$ MainFault 1 True 1 fault
2 0.487279 0.5001 0.602664 ${\bf{x}}_{\alpha \,{\bf{1}},2}$ MainFault 1 True 1 fault
3 0.435997 0.5001 0.448818 ${\bf{x}}_{\alpha \,{\bf{1}},3}$ MainFault 1 True 1 fault
4 0.461638 0.5001 0.525741 ${\bf{x}}_{\alpha \,{\bf{1}},4}$ MainFault 1 True 1 fault
5 0.538562 0.5001 0.628305 ${\bf{x}}_{\alpha \,{\bf{2}},0}$ SecondaryReservoir 2 False 2 Rest
6 0.359074 0.5001 0.641126 ${\bf{x}}_{\alpha \,{\bf{2}},1}$ SecondaryReservoir 2 False 2 Rest
7 0.435997 0.5001 0.615485 ${\bf{x}}_{\alpha \,{\bf{2}},2}$ SecondaryReservoir 2 False 2 Rest
8 0.615485 0.5001 0.602664 ${\bf{x}}_{\alpha \,{\bf{2}},3}$ SecondaryReservoir 2 False 2 Rest
9 0.692408 0.5001 0.577023 ${\bf{x}}_{\alpha \,{\bf{2}},4}$ SecondaryReservoir 2 False 2 Rest
10 0.718049 0.5001 0.583433 ${\bf{x}}_{\alpha \,{\bf{2}},5}$ SecondaryReservoir 2 False 2 Rest
11 0.718049 0.5001 0.557792 ${\bf{x}}_{\alpha \,{\bf{3}},0}$ Seal 3 False 2 Rest
12 0.359074 0.5001 0.615485 ${\bf{x}}_{\alpha \,{\bf{3}},1}$ Seal 3 False 2 Rest
13 0.538562 0.5001 0.602664 ${\bf{x}}_{\alpha \,{\bf{3}},2}$ Seal 3 False 2 Rest
14 0.692408 0.5001 0.551382 ${\bf{x}}_{\alpha \,{\bf{3}},3}$ Seal 3 False 2 Rest
15 0.435997 0.5001 0.589844 ${\bf{x}}_{\alpha \,{\bf{3}},4}$ Seal 3 False 2 Rest
16 0.615485 0.5001 0.577023 ${\bf{x}}_{\alpha \,{\bf{3}},5}$ Seal 3 False 2 Rest
17 0.410356 0.294972 0.423177 ${\bf{x}}_{\alpha \,{\bf{4}},0}$ Reservoir 4 False 2 Rest
18 0.410356 0.705228 0.423177 ${\bf{x}}_{\alpha \,{\bf{4}},1}$ Reservoir 4 False 2 Rest
19 0.359074 0.5001 0.512921 ${\bf{x}}_{\alpha \,{\bf{4}},2}$ Reservoir 4 False 2 Rest
20 0.282151 0.5001 0.538562 ${\bf{x}}_{\alpha \,{\bf{4}},3}$ Reservoir 4 False 2 Rest
21 0.718049 0.5001 0.455228 ${\bf{x}}_{\alpha \,{\bf{4}},4}$ Reservoir 4 False 2 Rest
22 0.692408 0.5001 0.448818 ${\bf{x}}_{\alpha \,{\bf{4}},5}$ Reservoir 4 False 2 Rest
23 0.538562 0.5001 0.5001 ${\bf{x}}_{\alpha \,{\bf{4}},6}$ Reservoir 4 False 2 Rest
24 0.615485 0.5001 0.474459 ${\bf{x}}_{\alpha \,{\bf{4}},7}$ Reservoir 4 False 2 Rest
25 0.653946 0.5001 0.461638 ${\bf{x}}_{\alpha \,{\bf{4}},8}$ Reservoir 4 False 2 Rest
26 0.461638 0.269331 0.5001 ${\bf{x}}_{\alpha \,{\bf{4}},9}$ Reservoir 4 False 2 Rest
27 0.461638 0.7501 0.512921 ${\bf{x}}_{\alpha \,{\bf{4}},10}$ Reservoir 4 False 2 Rest
28 0.461638 0.730869 0.5001 ${\bf{x}}_{\alpha \,{\bf{4}},11}$ Reservoir 4 False 2 Rest
29 0.5001 0.5001 0.512921 ${\bf{x}}_{\alpha \,{\bf{4}},12}$ Reservoir 4 False 2 Rest
30 0.461638 0.2501 0.512921 ${\bf{x}}_{\alpha \,{\bf{4}},13}$ Reservoir 4 False 2 Rest
31 0.359074 0.5001 0.435997 ${\bf{x}}_{\alpha \,{\bf{5}},0}$ Overlying 5 False 2 Rest
32 0.282151 0.5001 0.461638 ${\bf{x}}_{\alpha \,{\bf{5}},1}$ Overlying 5 False 2 Rest
33 0.461638 0.5001 0.448818 ${\bf{x}}_{\alpha \,{\bf{5}},2}$ Overlying 5 False 2 Rest
34 0.718049 0.5001 0.378305 ${\bf{x}}_{\alpha \,{\bf{5}},3}$ Overlying 5 False 2 Rest
35 0.538562 0.5001 0.423177 ${\bf{x}}_{\alpha \,{\bf{5}},4}$ Overlying 5 False 2 Rest
36 0.615485 0.5001 0.397536 ${\bf{x}}_{\alpha \,{\bf{5}},5}$ Overlying 5 False 2 Rest
37 0.653946 0.5001 0.384715 ${\bf{x}}_{\alpha \,{\bf{5}},6}$ Overlying 5 False 2 Rest
38 0.692408 0.5001 0.371895 ${\bf{x}}_{\alpha \,{\bf{5}},7}$ Overlying 5 False 2 Rest
39 0.666767 0.5001 0.320613 ${\bf{x}}_{\alpha \,{\bf{6}},0}$ Layer 3 6 False 3 unco
40 0.512921 0.5001 0.320613 ${\bf{x}}_{\alpha \,{\bf{6}},1}$ Layer 3 6 False 3 unco
lith_block, fault_block = gp.compute_model(interp_data)

This solution can be plot with the correspondent plotting function. Blocks:

%matplotlib inline
gp.plot_section(geo_data, lith_block[0], 25, plot_data=True)
../_images/ch1-Copy1_33_0.png
%matplotlib inline
gp.plot_section(geo_data, fault_block[0], 25, plot_data=True)
../_images/ch1-Copy1_34_0.png

Potential field:

gp.plot_potential_field(geo_data, lith_block[1], 25, N=500, cmap='viridis')
../_images/ch1-Copy1_36_0.png

From the potential fields (of lithologies and faults) it is possible to extract vertices and simpleces to create the 3D triangles for a vtk visualization.

ver, sim = gp.get_surfaces(interp_data,lith_block[1], fault_block[1], original_scale=True)
---------------------------------------------------------------------------

ValueError                                Traceback (most recent call last)

<ipython-input-49-b70d47564eaf> in <module>()
----> 1 ver, sim = gp.get_surfaces(interp_data,lith_block[1], fault_block[1], original_scale=True)


~/PycharmProjects/gempy/gempy/GemPy_f.py in get_surfaces(interp_data, potential_lith, potential_fault, n_formation, step_size, original_scale)
    580             else:
    581                 v, s = get_surface(potential_fault, interp_data, pot_int, n,
--> 582                                    step_size=step_size, original_scale=original_scale)
    583                 vertices.append(v)
    584                 simplices.append(s)


~/PycharmProjects/gempy/gempy/GemPy_f.py in get_surface(potential_block, interp_data, pot_int, n_formation, step_size, original_scale)
    556             spacing=((interp_data.geo_data_res.extent[1] - interp_data.geo_data_res.extent[0]) / interp_data.geo_data_res.resolution[0],
    557                      (interp_data.geo_data_res.extent[3] - interp_data.geo_data_res.extent[2]) / interp_data.geo_data_res.resolution[1],
--> 558                      (interp_data.geo_data_res.extent[5] - interp_data.geo_data_res.extent[4]) / interp_data.geo_data_res.resolution[2]))
    559
    560


~/anaconda3/lib/python3.6/site-packages/skimage/measure/_marching_cubes_lewiner.py in marching_cubes_lewiner(volume, level, spacing, gradient_direction, step_size, allow_degenerate, use_classic)
    197         level = float(level)
    198         if level < volume.min() or level > volume.max():
--> 199             raise ValueError("Surface level must be within volume data range.")
    200     # spacing
    201     if len(spacing) != 3:


ValueError: Surface level must be within volume data range.
%debug
> /home/miguel/anaconda3/lib/python3.6/site-packages/skimage/measure/_marching_cubes_lewiner.py(199)marching_cubes_lewiner()
    197         level = float(level)
    198         if level < volume.min() or level > volume.max():
--> 199             raise ValueError("Surface level must be within volume data range.")
    200     # spacing
    201     if len(spacing) != 3:

ipdb> up
> /home/miguel/PycharmProjects/gempy/gempy/GemPy_f.py(558)get_surface()
    556             spacing=((interp_data.geo_data_res.extent[1] - interp_data.geo_data_res.extent[0]) / interp_data.geo_data_res.resolution[0],
    557                      (interp_data.geo_data_res.extent[3] - interp_data.geo_data_res.extent[2]) / interp_data.geo_data_res.resolution[1],
--> 558                      (interp_data.geo_data_res.extent[5] - interp_data.geo_data_res.extent[4]) / interp_data.geo_data_res.resolution[2]))
    559 
    560 

ipdb> pot_int[n_formation-1]
90.703789
ipdb> potential_block.max()
95.060684
ipdb> potential_block.min()
90.841385
ipdb> exit
gp.plot_surfaces_3D(geo_data, ver, sim, alpha=1)

Additionally is possible to update the model and recompute the surfaces in real time. To do so, we need to pass the data rescaled. To get an smooth response is important to have the theano optimizer flag in fast_run and run theano in the gpu. This can speed up the modeling time in a factor of 20.

gp.get_data(interp_data.geo_data_res, verbosity=1)
G_x G_y G_z X Y Z annotations azimuth dip formation formation number isFault order_series polarity series
interfaces 0 NaN NaN NaN 0.474459 0.5001 0.49369 ${\bf{x}}_{\alpha \,{\bf{1}},0}$ NaN NaN MainFault 1 True 1 NaN fault
1 NaN NaN NaN 0.423177 0.5001 0.339844 ${\bf{x}}_{\alpha \,{\bf{1}},1}$ NaN NaN MainFault 1 True 1 NaN fault
2 NaN NaN NaN 0.525741 0.5001 0.647536 ${\bf{x}}_{\alpha \,{\bf{1}},2}$ NaN NaN MainFault 1 True 1 NaN fault
3 NaN NaN NaN 0.5001 0.5001 0.570613 ${\bf{x}}_{\alpha \,{\bf{1}},3}$ NaN NaN MainFault 1 True 1 NaN fault
4 NaN NaN NaN 0.448818 0.5001 0.416767 ${\bf{x}}_{\alpha \,{\bf{1}},4}$ NaN NaN MainFault 1 True 1 NaN fault
5 NaN NaN NaN 0.679587 0.5001 0.314203 ${\bf{x}}_{\alpha \,{\bf{2}},0}$ NaN NaN fault2 2 True 2 NaN Fault_serie2
6 NaN NaN NaN 0.269331 0.5001 0.711638 ${\bf{x}}_{\alpha \,{\bf{2}},1}$ NaN NaN fault2 2 True 2 NaN Fault_serie2
7 NaN NaN NaN 0.705228 0.5001 0.544972 ${\bf{x}}_{\alpha \,{\bf{3}},0}$ NaN NaN SecondaryReservoir 3 False 3 NaN Rest
8 NaN NaN NaN 0.628305 0.5001 0.570613 ${\bf{x}}_{\alpha \,{\bf{3}},1}$ NaN NaN SecondaryReservoir 3 False 3 NaN Rest
9 NaN NaN NaN 0.371895 0.5001 0.609074 ${\bf{x}}_{\alpha \,{\bf{3}},2}$ NaN NaN SecondaryReservoir 3 False 3 NaN Rest
10 NaN NaN NaN 0.551382 0.5001 0.596254 ${\bf{x}}_{\alpha \,{\bf{3}},3}$ NaN NaN SecondaryReservoir 3 False 3 NaN Rest
11 NaN NaN NaN 0.448818 0.5001 0.583433 ${\bf{x}}_{\alpha \,{\bf{3}},4}$ NaN NaN SecondaryReservoir 3 False 3 NaN Rest
12 NaN NaN NaN 0.730869 0.5001 0.551382 ${\bf{x}}_{\alpha \,{\bf{3}},5}$ NaN NaN SecondaryReservoir 3 False 3 NaN Rest
13 NaN NaN NaN 0.448818 0.5001 0.557792 ${\bf{x}}_{\alpha \,{\bf{4}},0}$ NaN NaN Seal 4 False 3 NaN Rest
14 NaN NaN NaN 0.730869 0.5001 0.525741 ${\bf{x}}_{\alpha \,{\bf{4}},1}$ NaN NaN Seal 4 False 3 NaN Rest
15 NaN NaN NaN 0.371895 0.5001 0.583433 ${\bf{x}}_{\alpha \,{\bf{4}},2}$ NaN NaN Seal 4 False 3 NaN Rest
16 NaN NaN NaN 0.551382 0.5001 0.570613 ${\bf{x}}_{\alpha \,{\bf{4}},3}$ NaN NaN Seal 4 False 3 NaN Rest
17 NaN NaN NaN 0.628305 0.5001 0.544972 ${\bf{x}}_{\alpha \,{\bf{4}},4}$ NaN NaN Seal 4 False 3 NaN Rest
18 NaN NaN NaN 0.705228 0.5001 0.519331 ${\bf{x}}_{\alpha \,{\bf{4}},5}$ NaN NaN Seal 4 False 3 NaN Rest
19 NaN NaN NaN 0.423177 0.294972 0.391126 ${\bf{x}}_{\alpha \,{\bf{5}},0}$ NaN NaN Reservoir 5 False 3 NaN Rest
20 NaN NaN NaN 0.423177 0.705228 0.391126 ${\bf{x}}_{\alpha \,{\bf{5}},1}$ NaN NaN Reservoir 5 False 3 NaN Rest
21 NaN NaN NaN 0.364376 0.504601 0.481594 ${\bf{x}}_{\alpha \,{\bf{5}},2}$ NaN NaN Reservoir 5 False 3 NaN Rest
22 NaN NaN NaN 0.294972 0.5001 0.50651 ${\bf{x}}_{\alpha \,{\bf{5}},3}$ NaN NaN Reservoir 5 False 3 NaN Rest
23 NaN NaN NaN 0.730869 0.5001 0.423177 ${\bf{x}}_{\alpha \,{\bf{5}},4}$ NaN NaN Reservoir 5 False 3 NaN Rest
24 NaN NaN NaN 0.705228 0.5001 0.416767 ${\bf{x}}_{\alpha \,{\bf{5}},5}$ NaN NaN Reservoir 5 False 3 NaN Rest
25 NaN NaN NaN 0.551382 0.5001 0.468049 ${\bf{x}}_{\alpha \,{\bf{5}},6}$ NaN NaN Reservoir 5 False 3 NaN Rest
26 NaN NaN NaN 0.628305 0.5001 0.442408 ${\bf{x}}_{\alpha \,{\bf{5}},7}$ NaN NaN Reservoir 5 False 3 NaN Rest
27 NaN NaN NaN 0.666767 0.5001 0.429587 ${\bf{x}}_{\alpha \,{\bf{5}},8}$ NaN NaN Reservoir 5 False 3 NaN Rest
28 NaN NaN NaN 0.474459 0.730869 0.468049 ${\bf{x}}_{\alpha \,{\bf{5}},9}$ NaN NaN Reservoir 5 False 3 NaN Rest
29 NaN NaN NaN 0.512921 0.5001 0.480869 ${\bf{x}}_{\alpha \,{\bf{5}},10}$ NaN NaN Reservoir 5 False 3 NaN Rest
30 NaN NaN NaN 0.474459 0.269331 0.468049 ${\bf{x}}_{\alpha \,{\bf{5}},11}$ NaN NaN Reservoir 5 False 3 NaN Rest
31 NaN NaN NaN 0.474459 0.7501 0.480869 ${\bf{x}}_{\alpha \,{\bf{5}},12}$ NaN NaN Reservoir 5 False 3 NaN Rest
32 NaN NaN NaN 0.474459 0.2501 0.480869 ${\bf{x}}_{\alpha \,{\bf{5}},13}$ NaN NaN Reservoir 5 False 3 NaN Rest
33 NaN NaN NaN 0.474459 0.5001 0.416767 ${\bf{x}}_{\alpha \,{\bf{6}},0}$ NaN NaN Overlying 6 False 3 NaN Rest
34 NaN NaN NaN 0.730869 0.5001 0.346254 ${\bf{x}}_{\alpha \,{\bf{6}},1}$ NaN NaN Overlying 6 False 3 NaN Rest
35 NaN NaN NaN 0.628305 0.5001 0.365485 ${\bf{x}}_{\alpha \,{\bf{6}},2}$ NaN NaN Overlying 6 False 3 NaN Rest
36 NaN NaN NaN 0.666767 0.5001 0.352664 ${\bf{x}}_{\alpha \,{\bf{6}},3}$ NaN NaN Overlying 6 False 3 NaN Rest
37 NaN NaN NaN 0.705228 0.5001 0.339844 ${\bf{x}}_{\alpha \,{\bf{6}},4}$ NaN NaN Overlying 6 False 3 NaN Rest
38 NaN NaN NaN 0.551382 0.5001 0.391126 ${\bf{x}}_{\alpha \,{\bf{6}},5}$ NaN NaN Overlying 6 False 3 NaN Rest
39 NaN NaN NaN 0.371895 0.5001 0.403946 ${\bf{x}}_{\alpha \,{\bf{6}},6}$ NaN NaN Overlying 6 False 3 NaN Rest
40 NaN NaN NaN 0.294972 0.5001 0.429587 ${\bf{x}}_{\alpha \,{\bf{6}},7}$ NaN NaN Overlying 6 False 3 NaN Rest
41 NaN NaN NaN 0.679587 0.5001 0.288562 ${\bf{x}}_{\alpha \,{\bf{7}},0}$ NaN NaN Layer 3 7 False 4 NaN unco
42 NaN NaN NaN 0.525741 0.5001 0.339844 ${\bf{x}}_{\alpha \,{\bf{7}},1}$ NaN NaN Layer 3 7 False 4 NaN unco
foliations 0 -0.948683 -1.7427e-16 0.316229 0.453292 0.5001 0.458972 ${\bf{x}}_{\beta \,{\bf{1}},0}$ 270 71.565 MainFault 1 True 1 1 fault
1 0.642788 3.93594e-17 0.766044 0.577023 0.5001 0.416767 ${\bf{x}}_{\beta \,{\bf{2}},0}$ 90 40 fault2 2 True 2 1 Fault_serie2
2 0.316229 1.93634e-17 0.948683 0.589844 0.5001 0.455228 ${\bf{x}}_{\beta \,{\bf{5}},0}$ 90 18.435 Reservoir 5 False 3 1 Rest
3 0.34202 2.09427e-17 0.939693 0.602664 0.5001 0.391126 ${\bf{x}}_{\beta \,{\bf{7}},0}$ 90 20 Layer 3 7 False 4 1 unco
ver_s, sim_s = gp.get_surfaces(interp_data,lith_block[1],
                               fault_block[1],
                               original_scale=False)
---------------------------------------------------------------------------

ValueError                                Traceback (most recent call last)

<ipython-input-50-88ef9c3214df> in <module>()
      1 ver_s, sim_s = gp.get_surfaces(interp_data,lith_block[1],
      2                                fault_block[1],
----> 3                                original_scale=False)


~/PycharmProjects/gempy/gempy/GemPy_f.py in get_surfaces(interp_data, potential_lith, potential_fault, n_formation, step_size, original_scale)
    580             else:
    581                 v, s = get_surface(potential_fault, interp_data, pot_int, n,
--> 582                                    step_size=step_size, original_scale=original_scale)
    583                 vertices.append(v)
    584                 simplices.append(s)


~/PycharmProjects/gempy/gempy/GemPy_f.py in get_surface(potential_block, interp_data, pot_int, n_formation, step_size, original_scale)
    556             spacing=((interp_data.geo_data_res.extent[1] - interp_data.geo_data_res.extent[0]) / interp_data.geo_data_res.resolution[0],
    557                      (interp_data.geo_data_res.extent[3] - interp_data.geo_data_res.extent[2]) / interp_data.geo_data_res.resolution[1],
--> 558                      (interp_data.geo_data_res.extent[5] - interp_data.geo_data_res.extent[4]) / interp_data.geo_data_res.resolution[2]))
    559
    560


~/anaconda3/lib/python3.6/site-packages/skimage/measure/_marching_cubes_lewiner.py in marching_cubes_lewiner(volume, level, spacing, gradient_direction, step_size, allow_degenerate, use_classic)
    197         level = float(level)
    198         if level < volume.min() or level > volume.max():
--> 199             raise ValueError("Surface level must be within volume data range.")
    200     # spacing
    201     if len(spacing) != 3:


ValueError: Surface level must be within volume data range.
gp.plot_surfaces_3D_real_time(interp_data, ver_s, sim_s)

In the same manner we can visualize the fault block:

gp.plot_section(geo_data, fault_block[0], 25)
../_images/ch1-Copy1_46_0.png
gp.plot_potential_field(geo_data, fault_block[1], 25)
../_images/ch1-Copy1_47_0.png