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                 1.7e+03
$C_o$                 7.1e+04
drift equations  [3, 3, 3, 3]
nugget grad              0.01
nugget scalar           1e-06
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 0x7ff9aec11f10>
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                 1.7e+03
$C_o$                 7.1e+04
drift equations  [3, 3, 3, 3]
nugget grad              0.01
nugget scalar           1e-06
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 1.7e+03
$C_o$ 7.1e+04
drift equations [3, 3, 3, 3]
nugget grad 0.01
nugget scalar 1e-06
Rescaling rescaling factor 1.6e+03
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 0x7ff9f8dab9d0>
gp.plot_2d(geo_model_graben, direction='x')
Cell Number: mid Direction: x

Out:

<gempy.plot.visualization_2d.Plot2D object at 0x7ff9ba7c27d0>
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 0x7ff9fad72050>

sphinx_gallery_thumbnail_number = 5

ch1 5 fault relations ch1 5 fault relations

Out:

<gempy.plot.vista.GemPyToVista object at 0x7ff9bba32650>
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 0x7ff9bb940210>

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 0x7ff9bb903a10>
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 0x7ff9bba3b2d0>

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 0x7ff9ba647210>

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 36.765 seconds)

Gallery generated by Sphinx-Gallery