1.5: Fault relations

Importing gempy

import gempy as gp

# Aux imports
import numpy as np
import pandas as pd
import os

# Importing the function to find the interface
from gempy.utils.input_manipulation import find_interfaces_from_block_bottoms
import matplotlib.pyplot as plt

np.random.seed(1515)
pd.set_option('precision', 2)

We import a model from an existing folder.

cwd = os.getcwd()
if 'examples' not in cwd:
    data_path = os.getcwd() + '/examples/'
else:
    data_path = cwd + '/../../'

geo_model = gp.load_model(r'Tutorial_ch1-9a_Fault_relations',
                          path=data_path + 'data/gempy_models/Tutorial_ch1-9a_Fault_relations',
                          recompile=True)

Out:

Active grids: ['regular']
Setting kriging parameters to their default values.
Compiling theano function...
Level of Optimization:  fast_compile
Device:  cpu
Precision:  float64
Number of faults:  1
Compilation Done!
Kriging values:
                        values
range                 1732.05
$C_o$                71428.57
drift equations  [3, 3, 3, 3]
nugget grad              0.01
nugget scalar             0.0
fault_series1 series2 series1 basement_series
fault_series1 False False True True
series2 False False False False
series1 False False False False
basement_series False False False False


order_series BottomRelation isFault isFinite isActive
fault_series1 1 Fault True False True
series2 2 Erosion False False True
series1 3 Erosion False False True
basement_series 4 Erosion False False False


surface series order_surfaces color id
5 fault1 fault_series1 1 #015482 1
7 rock5 series2 1 #9f0052 2
6 rock4 series2 2 #ffbe00 3
0 rock3 series1 1 #728f02 4
3 rock2 series1 2 #443988 5
1 rock1 series1 3 #ff3f20 6
4 basement basement_series 1 #325916 7


gp.compute_model(geo_model, compute_mesh=False)

Out:

Lithology ids
  [7. 7. 7. ... 2. 2. 2.]

Out:

array([7., 7., 7., ..., 2., 2., 2.])

Out:

array([[2., 2., 2., ..., 1., 1., 1.]])
gp.plot_2d(geo_model, cell_number=[25], show_data=True)

# Graben example
# --------------
Cell Number: 25 Direction: y

Out:

<gempy.plot.visualization_2d.Plot2D object at 0x7f3194e63280>
geo_model_graben = gp.load_model(r'Tutorial_ch1-9b_Fault_relations',
                                 path=data_path + 'data/gempy_models/Tutorial_ch1-9b_Fault_relations', recompile=True)

geo_model.meta.project_name = "Faults_relations"

Out:

Active grids: ['regular']
Setting kriging parameters to their default values.
Compiling theano function...
Level of Optimization:  fast_compile
Device:  cpu
Precision:  float64
Number of faults:  2
Compilation Done!
Kriging values:
                        values
range                 1732.05
$C_o$                71428.57
drift equations  [3, 3, 3, 3]
nugget grad              0.01
nugget scalar             0.0
surface series order_surfaces color id
7 fault2 fault_series2 1 #015482 1
5 fault1 fault_series1 1 #9f0052 2
6 rock4 series1 1 #ffbe00 3
0 rock3 series1 2 #728f02 4
3 rock2 series1 3 #443988 5
1 rock1 series1 4 #ff3f20 6
4 basement basement_series 1 #325916 7


values
Structure isLith True
isFault True
number faults 2
number surfaces 6
number series 4
number surfaces per series [1, 1, 4, 0]
len surfaces surface_points [9, 9, 15, 15, 12, 12]
len series surface_points [9, 9, 54, 0]
len series orientations [1, 1, 10, 0]
Options dtype float64
output geology
theano_optimizer fast_compile
device cpu
verbosity [0]
Kriging range 1732.05
$C_o$ 71428.57
drift equations [3, 3, 3, 3]
nugget grad 0.01
nugget scalar 0.0
Rescaling rescaling factor 1600.0
centers [500.0, 500.0, -650.0]


Displaying the input data:

gp.plot_2d(geo_model_graben, direction='y')
Cell Number: mid Direction: y

Out:

<gempy.plot.visualization_2d.Plot2D object at 0x7f31e2428f10>
gp.plot_2d(geo_model_graben, direction='x')
Cell Number: mid Direction: x

Out:

<gempy.plot.visualization_2d.Plot2D object at 0x7f31e24346a0>
order_series BottomRelation isFault isFinite isActive
fault_series2 1 Fault True False True
fault_series1 2 Fault True False True
series1 3 Erosion False False True
basement_series 4 Erosion False False False


order_series BottomRelation isFault isFinite isActive
fault_series2 1 Fault True False True
fault_series1 2 Fault True False True
series1 3 Erosion False False True
basement_series 4 Erosion False False False


fault_series2 fault_series1 series1 basement_series
fault_series2 False True True True
fault_series1 False False True True
series1 False False False False
basement_series False False False False


Out:

Lithology ids
  [7. 7. 7. ... 3. 3. 3.]
gp.plot_2d(geo_model_graben, cell_number=[25], show_data=True)
Cell Number: 25 Direction: y

Out:

<gempy.plot.visualization_2d.Plot2D object at 0x7f31fea41e50>

sphinx_gallery_thumbnail_number = 5

ch1 5 fault relations ch1 5 fault relations

Out:

<gempy.plot.vista.GemPyToVista object at 0x7f3233dbed00>
gp.plot_2d(geo_model_graben, cell_number=[25], show_scalar=True, series_n=0)

gp.plot_2d(geo_model_graben, cell_number=[25], show_scalar=True, series_n=1)
  • Cell Number: 25 Direction: y
  • Cell Number: 25 Direction: y

Out:

<gempy.plot.visualization_2d.Plot2D object at 0x7f321a978910>

Offset parameter (Experimental)

geo_model_graben._interpolator.theano_graph.offset.set_value(1)
gp.compute_model(geo_model_graben, compute_mesh=False)

Out:

Lithology ids
  [7. 7. 7. ... 3. 3. 3.]
Cell Number: mid Direction: y

Out:

<gempy.plot.visualization_2d.Plot2D object at 0x7f321a612670>
gp.plot_2d(geo_model_graben, series_n=2, show_scalar=True)
Cell Number: mid Direction: y

Out:

<gempy.plot.visualization_2d.Plot2D object at 0x7f3194be6d30>

Out:

array([0.17363373, 0.1848186 , 0.19600352, ..., 1.36450254, 1.37569005,
       1.38687763])

Finding the faults intersection:

Sometimes we need to find the voxels that contain the each fault. To do so we can use gempy’s functionality to find interfaces as follows. Lets use the first fault as an example:

gp.plot_2d(geo_model_graben,
           regular_grid=geo_model_graben.solutions.block_matrix[0, 0, :125000],
           show_data=True)
Cell Number: mid Direction: y

Out:

<gempy.plot.visualization_2d.Plot2D object at 0x7f3166f84340>

Remember the fault block is stored on:

Out:

array([1., 1., 1., ..., 2., 2., 2.])

Now we can find where is the intersection of the values 1 by calling the following function. This will return Trues on those voxels on the intersection

intersection = find_interfaces_from_block_bottoms(
    geo_model_graben.solutions.block_matrix[0, 0, :125000].reshape(50, 50, 50), 1, shift=1)

We can manually plotting together to see exactly what we have done

ax = gp.plot_2d(geo_model_graben,
                block=geo_model_graben.solutions.block_matrix[0, 0, :125000],
                show_data=True, show=False)

plt.imshow(intersection[:, 25, :].T, origin='lower', extent=(0, 1000, -1000, 0), alpha=.5)
plt.show()

gp.save_model(geo_model)
Cell Number: mid Direction: y

Out:

True

Total running time of the script: ( 0 minutes 15.158 seconds)

Gallery generated by Sphinx-Gallery