.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "Tutorials\02_POC_create-MC-ensemble.py"
.. LINE NUMBERS ARE GIVEN BELOW.
.. only:: html
.. note::
:class: sphx-glr-download-link-note
Click :ref:`here <sphx_glr_download_Tutorials_02_POC_create-MC-ensemble.py>`
to download the full example code
.. rst-class:: sphx-glr-example-title
.. _sphx_glr_Tutorials_02_POC_create-MC-ensemble.py:
02 - Monte Carlo simulation
======================
The following tutorial will lead you through an example workflow on how to create a Monte Carlo simulation of
geological models, meaning we will produce different geological geometries and also simulate their gravity
response.
.. GENERATED FROM PYTHON SOURCE LINES 10-13
Importing libraries
-------------------
First things first: let's import necessary libraries.
.. GENERATED FROM PYTHON SOURCE LINES 13-35
.. code-block:: default
import warnings
warnings.filterwarnings("ignore")
# Importing GemPy
import gempy as gp
from gempy.assets import topology as tp
from gempy.bayesian.fields import compute_prob, calculate_ie_masked
from gempy.assets.geophysics import GravityPreprocessing
# Importing auxilary libraries
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
plt.style.use('seaborn-talk')
#import sys
#sys.path.append('../../OpenWF/')
#from aux_functions import log_progress
# Check gempy version used for running the code
print(f"Code run with GemPy version: {gp.__version__}")
.. rst-class:: sphx-glr-script-out
Out:
.. code-block:: none
Code run with GemPy version: 2.2.9
.. GENERATED FROM PYTHON SOURCE LINES 36-41
Model Initialization
--------------------
First, we import the base Proof-of-Concept model (POC-model from here on), which was generated in the previous example. Using the loading method of GemPy `gp.load_model()` directly loads the model's input, already set with fault relations, surfaces assigned to a stack (series), etc.
Only thing left is to recompile and run the model.
.. GENERATED FROM PYTHON SOURCE LINES 41-50
.. code-block:: default
model_path = '../../models/2021-06-04_POC_base_model'
geo_model = gp.load_model('POC_PCT_model', path=model_path,
recompile=False)
# import DTM
dtm = np.load('../../models/Graben_base_model_topography.npy')
.. rst-class:: sphx-glr-script-out
Out:
.. code-block:: none
Active grids: ['regular']
Active grids: ['regular' 'topography']
.. GENERATED FROM PYTHON SOURCE LINES 51-52
Using the method `.get_additional_data()`, we can display a summary of model information and parameters, such as the kriging parameters.
.. GENERATED FROM PYTHON SOURCE LINES 52-55
.. code-block:: default
geo_model.get_additional_data()
.. raw:: html
<div class="output_subarea output_html rendered_html output_result">
<div>
<style scoped>
.dataframe tbody tr th:only-of-type {
vertical-align: middle;
}
.dataframe tbody tr th {
vertical-align: top;
}
.dataframe thead th {
text-align: right;
}
</style>
<table border="1" class="dataframe">
<thead>
<tr style="text-align: right;">
<th></th>
<th></th>
<th>values</th>
</tr>
</thead>
<tbody>
<tr>
<th rowspan="9" valign="top">Structure</th>
<th>isLith</th>
<td>True</td>
</tr>
<tr>
<th>isFault</th>
<td>True</td>
</tr>
<tr>
<th>number faults</th>
<td>5</td>
</tr>
<tr>
<th>number surfaces</th>
<td>11</td>
</tr>
<tr>
<th>number series</th>
<td>10</td>
</tr>
<tr>
<th>number surfaces per series</th>
<td>[1, 1, 1, 1, 1, 3, 1, 1, 1, 0]</td>
</tr>
<tr>
<th>len surfaces surface_points</th>
<td>[6, 8, 12, 8, 6, 8, 14, 14, 26, 24, 24]</td>
</tr>
<tr>
<th>len series surface_points</th>
<td>[6, 8, 12, 8, 6, 36, 26, 24, 24, 0]</td>
</tr>
<tr>
<th>len series orientations</th>
<td>[2, 2, 6, 4, 2, 12, 16, 12, 10, 0]</td>
</tr>
<tr>
<th rowspan="5" valign="top">Options</th>
<th>dtype</th>
<td>float64</td>
</tr>
<tr>
<th>output</th>
<td>geology</td>
</tr>
<tr>
<th>theano_optimizer</th>
<td>fast_compile</td>
</tr>
<tr>
<th>device</th>
<td>cpu</td>
</tr>
<tr>
<th>verbosity</th>
<td>None</td>
</tr>
<tr>
<th rowspan="3" valign="top">Kriging</th>
<th>range</th>
<td>32190.8</td>
</tr>
<tr>
<th>$C_o$</th>
<td>2.46726e+07</td>
</tr>
<tr>
<th>drift equations</th>
<td>[3, 3, 3, 3, 3, 3, 3, 3, 3, 3]</td>
</tr>
<tr>
<th rowspan="2" valign="top">Rescaling</th>
<th>rescaling factor</th>
<td>56916.7</td>
</tr>
<tr>
<th>centers</th>
<td>[14239.166495, 6900.0, -2581.92853]</td>
</tr>
</tbody>
</table>
</div>
</div>
<br />
<br />
.. GENERATED FROM PYTHON SOURCE LINES 56-57
Changing the kriging parameters affects the resulting models, e.g. the range represents the maximum correlation distance, or reducing the coefficient of correlation will yield a smoother, less "bumpy" model. For the POC-model, we set the `range` to 20000 and the correlation coefficient $C_o$ to 200000. Then we set up the interpolator, i.e. compile the functions which will calculate the scalar fields of our model surfaces.
.. GENERATED FROM PYTHON SOURCE LINES 57-75
.. code-block:: default
# adapt kriging to the parameters of previous example
# decrease the kriging range
geo_model.modify_kriging_parameters('range', 20000.)
geo_model.modify_kriging_parameters('$C_o$', 2e5)
# Set the interpolator function
# Create the theano model
gp.set_interpolator(geo_model,
compile_theano=True,
theano_optimizer='fast_compile',
verbose=[],
update_kriging=False);
# compute the model
sol = gp.compute_model(geo_model, compute_mesh=True)
.. rst-class:: sphx-glr-script-out
Out:
.. code-block:: none
Compiling theano function...
Level of Optimization: fast_compile
Device: cpu
Precision: float64
Number of faults: 5
Compilation Done!
Kriging values:
values
range 20000
$C_o$ 200000
drift equations [3, 3, 3, 3, 3, 3, 3, 3, 3, 3]
.. GENERATED FROM PYTHON SOURCE LINES 76-77
Now that the model is computed, lets have a look at a cross-section along the y-axis, so across the graben system:
.. GENERATED FROM PYTHON SOURCE LINES 77-81
.. code-block:: default
gp.plot_2d(geo_model, cell_number=25, direction='y', show_data=False, show_topography=False,
show_lith=True, show_results=True, show_boundaries=False);
.. image:: /Tutorials/images/sphx_glr_02_POC_create-MC-ensemble_001.png
:alt: Cell Number: 25 Direction: y
:class: sphx-glr-single-img
.. rst-class:: sphx-glr-script-out
Out:
.. code-block:: none
<gempy.plot.visualization_2d.Plot2D object at 0x0000017585940A60>
.. GENERATED FROM PYTHON SOURCE LINES 82-87
The two distinct domains in this model are directly visible: (i) the old graben system (extensional regime), covered by the (ii) thrusted, younger units.
Add Gravity grid
----------------
In the previous example, next to creating the model, we chose quasi-random locations for 15 gravity stations. The gravity signal of the base POC-model is simulated at these 15 stations. In the following workflows, we assume that these 15 stations were measured. So they serve as observed data for conditioning the MonteCarlo Ensemble of different geological geometries.
.. GENERATED FROM PYTHON SOURCE LINES 87-103
.. code-block:: default
# In[7]:
grav_stations = pd.read_csv('../../data/Data_for_MC/20210322_forw_grav_seed58.csv')
station_coordinates = np.stack((grav_stations.X.values,
grav_stations.Y.values,
grav_stations.Z.values), axis=1)
fig = plt.figure(figsize=[11,5])
cb = plt.scatter(grav_stations['X'], grav_stations['Y'], c=grav_stations['grav'],
marker='s', s=90, cmap='viridis')
plt.colorbar(cb, label='gravity')
plt.ylabel('y [m]')
plt.xlabel('x [m]');
.. image:: /Tutorials/images/sphx_glr_02_POC_create-MC-ensemble_002.png
:alt: 02 POC create MC ensemble
:class: sphx-glr-single-img
.. rst-class:: sphx-glr-script-out
Out:
.. code-block:: none
Text(0.5, 13.944444444444438, 'x [m]')
.. GENERATED FROM PYTHON SOURCE LINES 104-105
These stations are used for creating a centered grid around each station. The centered grid has an extent of 10 cells in x- and y-direction, and 15 cells in the z-direction.
.. GENERATED FROM PYTHON SOURCE LINES 105-112
.. code-block:: default
geo_model.set_centered_grid(station_coordinates, resolution = [10, 10, 15], radius=6000)
g = GravityPreprocessing(geo_model.grid.centered_grid)
tz = g.set_tz_kernel()
.. rst-class:: sphx-glr-script-out
Out:
.. code-block:: none
Active grids: ['regular' 'topography' 'centered']
.. GENERATED FROM PYTHON SOURCE LINES 113-118
We see that there are three active grids. On each, the gravity signal will be calculated. Of course, we can let it be calculated on each grid, but we may not need the information on e.g. the topography grid (which would for instance yield the geological map).
So we can set only the centered grid to active, which speeds up the simulation.
**Note** that you'll need to model also the regular grid, if you plan to export the `lith_block` geological voxel model later on!
As we want to also have the geometric changes in the lithological grid, we set `reset=False`. If we were to set it to `True`, only the 'centered' grid would be active.
.. GENERATED FROM PYTHON SOURCE LINES 118-121
.. code-block:: default
geo_model.set_active_grid('centered', reset=False)
.. rst-class:: sphx-glr-script-out
Out:
.. code-block:: none
Active grids: ['regular' 'topography' 'centered']
Grid Object. Values:
array([[ 140. , 140. , -6437.5 ],
[ 140. , 140. , -6312.5 ],
[ 140. , 140. , -6187.5 ],
...,
[30888.89 , 6285.71 , -3495.98176905],
[30888.89 , 6285.71 , -4948.49684561],
[30888.89 , 6285.71 , -6966.76 ]])
.. GENERATED FROM PYTHON SOURCE LINES 122-125
The centered grid will now be the only one where the model information is stored, meaning less computational time. Let's have a look how this comes in handy, when we start to modify the depth of units and calculate the gravity.
Before running the simulations, we need to assign densities to the rock units, otherwise it will raise an error.
.. GENERATED FROM PYTHON SOURCE LINES 125-131
.. code-block:: default
# add densities - from abdelfettah 2014 and SAPHYR
densities = [0, 0, 0, 0, 0, 2.466, 2.61, 2.53,
2.61, 2.47, 2.55, 2.67]
geo_model.add_surface_values(densities, ['density'])
.. raw:: html
<div class="output_subarea output_html rendered_html output_result">
<style type="text/css" >
#T_2f28591f_373c_11ec_9a60_00e04c6800carow0_col3 {
background-color: #5DA629;
} #T_2f28591f_373c_11ec_9a60_00e04c6800carow1_col3 {
background-color: #5DA629;
} #T_2f28591f_373c_11ec_9a60_00e04c6800carow2_col3 {
background-color: #015482;
} #T_2f28591f_373c_11ec_9a60_00e04c6800carow3_col3 {
background-color: #015482;
} #T_2f28591f_373c_11ec_9a60_00e04c6800carow4_col3 {
background-color: #015482;
} #T_2f28591f_373c_11ec_9a60_00e04c6800carow5_col3 {
background-color: #dbdbac;
} #T_2f28591f_373c_11ec_9a60_00e04c6800carow6_col3 {
background-color: #e588f3;
} #T_2f28591f_373c_11ec_9a60_00e04c6800carow7_col3 {
background-color: #ff792b;
} #T_2f28591f_373c_11ec_9a60_00e04c6800carow8_col3 {
background-color: #725c9a;
} #T_2f28591f_373c_11ec_9a60_00e04c6800carow9_col3 {
background-color: #cfc199;
} #T_2f28591f_373c_11ec_9a60_00e04c6800carow10_col3 {
background-color: #a5d490;
} #T_2f28591f_373c_11ec_9a60_00e04c6800carow11_col3 {
background-color: #c7848f;
}</style><table id="T_2f28591f_373c_11ec_9a60_00e04c6800ca" ><thead> <tr> <th class="blank level0" ></th> <th class="col_heading level0 col0" >surface</th> <th class="col_heading level0 col1" >series</th> <th class="col_heading level0 col2" >order_surfaces</th> <th class="col_heading level0 col3" >color</th> <th class="col_heading level0 col4" >id</th> <th class="col_heading level0 col5" >density</th> </tr></thead><tbody>
<tr>
<th id="T_2f28591f_373c_11ec_9a60_00e04c6800calevel0_row0" class="row_heading level0 row0" >9</th>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow0_col0" class="data row0 col0" >Thrust1_south</td>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow0_col1" class="data row0 col1" >Thrust1_series</td>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow0_col2" class="data row0 col2" >1</td>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow0_col3" class="data row0 col3" >#5DA629</td>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow0_col4" class="data row0 col4" >1</td>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow0_col5" class="data row0 col5" >0.000000</td>
</tr>
<tr>
<th id="T_2f28591f_373c_11ec_9a60_00e04c6800calevel0_row1" class="row_heading level0 row1" >10</th>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow1_col0" class="data row1 col0" >Thrust2_south</td>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow1_col1" class="data row1 col1" >Thrust2_series</td>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow1_col2" class="data row1 col2" >1</td>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow1_col3" class="data row1 col3" >#5DA629</td>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow1_col4" class="data row1 col4" >2</td>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow1_col5" class="data row1 col5" >0.000000</td>
</tr>
<tr>
<th id="T_2f28591f_373c_11ec_9a60_00e04c6800calevel0_row2" class="row_heading level0 row2" >0</th>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow2_col0" class="data row2 col0" >Fault2</td>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow2_col1" class="data row2 col1" >Fault2_series</td>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow2_col2" class="data row2 col2" >1</td>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow2_col3" class="data row2 col3" >#015482</td>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow2_col4" class="data row2 col4" >3</td>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow2_col5" class="data row2 col5" >0.000000</td>
</tr>
<tr>
<th id="T_2f28591f_373c_11ec_9a60_00e04c6800calevel0_row3" class="row_heading level0 row3" >1</th>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow3_col0" class="data row3 col0" >Fault5</td>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow3_col1" class="data row3 col1" >Fault5_series</td>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow3_col2" class="data row3 col2" >1</td>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow3_col3" class="data row3 col3" >#015482</td>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow3_col4" class="data row3 col4" >4</td>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow3_col5" class="data row3 col5" >0.000000</td>
</tr>
<tr>
<th id="T_2f28591f_373c_11ec_9a60_00e04c6800calevel0_row4" class="row_heading level0 row4" >2</th>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow4_col0" class="data row4 col0" >Fault6</td>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow4_col1" class="data row4 col1" >Fault6_series</td>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow4_col2" class="data row4 col2" >1</td>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow4_col3" class="data row4 col3" >#015482</td>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow4_col4" class="data row4 col4" >5</td>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow4_col5" class="data row4 col5" >0.000000</td>
</tr>
<tr>
<th id="T_2f28591f_373c_11ec_9a60_00e04c6800calevel0_row5" class="row_heading level0 row5" >6</th>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow5_col0" class="data row5 col0" >Tertiary</td>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow5_col1" class="data row5 col1" >Post_tectonic_series</td>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow5_col2" class="data row5 col2" >1</td>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow5_col3" class="data row5 col3" >#dbdbac</td>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow5_col4" class="data row5 col4" >6</td>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow5_col5" class="data row5 col5" >2.466000</td>
</tr>
<tr>
<th id="T_2f28591f_373c_11ec_9a60_00e04c6800calevel0_row6" class="row_heading level0 row6" >8</th>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow6_col0" class="data row6 col0" >Pink</td>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow6_col1" class="data row6 col1" >Post_tectonic_series</td>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow6_col2" class="data row6 col2" >2</td>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow6_col3" class="data row6 col3" >#e588f3</td>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow6_col4" class="data row6 col4" >7</td>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow6_col5" class="data row6 col5" >2.610000</td>
</tr>
<tr>
<th id="T_2f28591f_373c_11ec_9a60_00e04c6800calevel0_row7" class="row_heading level0 row7" >7</th>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow7_col0" class="data row7 col0" >Orange</td>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow7_col1" class="data row7 col1" >Post_tectonic_series</td>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow7_col2" class="data row7 col2" >3</td>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow7_col3" class="data row7 col3" >#ff792b</td>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow7_col4" class="data row7 col4" >8</td>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow7_col5" class="data row7 col5" >2.530000</td>
</tr>
<tr>
<th id="T_2f28591f_373c_11ec_9a60_00e04c6800calevel0_row8" class="row_heading level0 row8" >5</th>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow8_col0" class="data row8 col0" >Unconformity</td>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow8_col1" class="data row8 col1" >Detachement</td>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow8_col2" class="data row8 col2" >1</td>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow8_col3" class="data row8 col3" >#725c9a</td>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow8_col4" class="data row8 col4" >9</td>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow8_col5" class="data row8 col5" >2.610000</td>
</tr>
<tr>
<th id="T_2f28591f_373c_11ec_9a60_00e04c6800calevel0_row9" class="row_heading level0 row9" >4</th>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow9_col0" class="data row9 col0" >Upper-filling</td>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow9_col1" class="data row9 col1" >Syn_tectonic_series2</td>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow9_col2" class="data row9 col2" >1</td>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow9_col3" class="data row9 col3" >#cfc199</td>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow9_col4" class="data row9 col4" >10</td>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow9_col5" class="data row9 col5" >2.470000</td>
</tr>
<tr>
<th id="T_2f28591f_373c_11ec_9a60_00e04c6800calevel0_row10" class="row_heading level0 row10" >3</th>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow10_col0" class="data row10 col0" >Lower-filling</td>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow10_col1" class="data row10 col1" >Pre_tectonic_series</td>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow10_col2" class="data row10 col2" >1</td>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow10_col3" class="data row10 col3" >#a5d490</td>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow10_col4" class="data row10 col4" >11</td>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow10_col5" class="data row10 col5" >2.550000</td>
</tr>
<tr>
<th id="T_2f28591f_373c_11ec_9a60_00e04c6800calevel0_row11" class="row_heading level0 row11" >11</th>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow11_col0" class="data row11 col0" >basement</td>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow11_col1" class="data row11 col1" >Basement</td>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow11_col2" class="data row11 col2" >1</td>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow11_col3" class="data row11 col3" >#c7848f</td>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow11_col4" class="data row11 col4" >12</td>
<td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow11_col5" class="data row11 col5" >2.670000</td>
</tr>
</tbody></table>
</div>
<br />
<br />
.. GENERATED FROM PYTHON SOURCE LINES 132-135
MC Variation
------------
For varying the depth of units, we extract the indices of the units whose input points we want to modify. To guarantee that we always vary the original depth in each realization (and not the depth used in the previous realization), we first generate an initial-depth array, containing the original depth information of all input points:
.. GENERATED FROM PYTHON SOURCE LINES 135-138
.. code-block:: default
Z_init = geo_model.surface_points.df['Z'].copy()
.. GENERATED FROM PYTHON SOURCE LINES 139-140
Having all the undisturbed depth values, we extract all surface points belonging to the units whose inputs we want to vary:
.. GENERATED FROM PYTHON SOURCE LINES 140-146
.. code-block:: default
graben_lower = geo_model.surface_points.df.query("surface=='Lower-filling'")
graben_middle = geo_model.surface_points.df.query("surface=='Upper-filling'")
unconformity = geo_model.surface_points.df.query("surface=='Unconformity'")
.. GENERATED FROM PYTHON SOURCE LINES 147-148
Before running the Monte Carlo simulations, we set up the interpolator for a "fast-run", i.e. it optimizes runtime on cost of compilation time:
.. GENERATED FROM PYTHON SOURCE LINES 148-152
.. code-block:: default
gp.set_interpolator(geo_model, output=['gravity'],
theano_optimizer='fast_run',
update_kriging=True)
.. rst-class:: sphx-glr-script-out
Out:
.. code-block:: none
Setting kriging parameters to their default values.
Compiling theano function...
Level of Optimization: fast_run
Device: cpu
Precision: float64
Number of faults: 5
Compilation Done!
Kriging values:
values
range 32190.8
$C_o$ 2.46726e+07
drift equations [3, 3, 3, 3, 3, 3, 3, 3, 3, 3]
<gempy.core.interpolator.InterpolatorModel object at 0x0000017586B9C340>
.. GENERATED FROM PYTHON SOURCE LINES 153-156
Now we are good to go and run the Monte Carlo simulation. In the following, we fix a numpy random number seed so that this MC-simulation is reproducible
Then, we create empty arrays and dictionaries for the lithologies and gravity, respectively. In a `for` loop, we then vary depths of interface points and
compute a model.
.. GENERATED FROM PYTHON SOURCE LINES 156-188
.. code-block:: default
np.random.seed(1)
# allocate array for lithology blocks
lith_blocks = np.array([])
# create a dictionary to store gravity of simulations
grav = dict()
# get indices where the variable input points are
Lgraben = list(graben_lower.index)
Ugraben = list(graben_middle.index)
Uncon = list(unconformity.index)
Cindices = Lgraben + Ugraben + Uncon
# set number of realizations
n_iterations = 10
for i in range(n_iterations):
# vary surface points
Z_var = np.random.normal(0, 300, size=3)
Z_loc = np.hstack([Z_init[Lgraben] + Z_var[0],
Z_init[Ugraben] + Z_var[1],
Z_init[Uncon] + Z_var[2]])
# apply variation to model
geo_model.modify_surface_points(Cindices, Z=Z_loc)
# re-compute model
gp.compute_model(geo_model)
# store lithologies ONLY THERE IF REGULAR GRID IS ACTIVE
lith_blocks = np.append(lith_blocks, geo_model.solutions.lith_block)
# store gravity
grav[f"Real_{i}"] = geo_model.solutions.fw_gravity
lith_blocks = lith_blocks.reshape(n_iterations, -1)
.. GENERATED FROM PYTHON SOURCE LINES 189-192
Export models and gravity
-------------------------
For post-processing of use in different software (e.g. numerical simulators for heat- and mass-transport), knowing ways of exporting the MC-results, in this case the simulated gravity and the lithology-blocks, comes in handy. There are many different ways of saving stuff (e.g. pickle the simulation results), but here we present simple exports as `.csv` and `.npy` files.
.. GENERATED FROM PYTHON SOURCE LINES 192-203
.. code-block:: default
gravdf = pd.DataFrame.from_dict(grav)
# add station coordinates to the dataframe
gravdf["X"] = station_coordinates[:,0]
gravdf["Y"] = station_coordinates[:,1]
gravdf["Z"] =station_coordinates[:,2]
gravdf.head()
.. raw:: html
<div class="output_subarea output_html rendered_html output_result">
<div>
<style scoped>
.dataframe tbody tr th:only-of-type {
vertical-align: middle;
}
.dataframe tbody tr th {
vertical-align: top;
}
.dataframe thead th {
text-align: right;
}
</style>
<table border="1" class="dataframe">
<thead>
<tr style="text-align: right;">
<th></th>
<th>Real_0</th>
<th>Real_1</th>
<th>Real_2</th>
<th>Real_3</th>
<th>Real_4</th>
<th>Real_5</th>
<th>Real_6</th>
<th>Real_7</th>
<th>Real_8</th>
<th>Real_9</th>
<th>X</th>
<th>Y</th>
<th>Z</th>
</tr>
</thead>
<tbody>
<tr>
<th>0</th>
<td>-533.673251</td>
<td>-532.673751</td>
<td>-534.038468</td>
<td>-533.143688</td>
<td>-534.590898</td>
<td>-533.483330</td>
<td>-533.631910</td>
<td>-534.186293</td>
<td>-533.873075</td>
<td>-534.097548</td>
<td>21777.78</td>
<td>7142.86</td>
<td>405.17</td>
</tr>
<tr>
<th>1</th>
<td>-533.469712</td>
<td>-531.718800</td>
<td>-534.009147</td>
<td>-532.039845</td>
<td>-534.347122</td>
<td>-533.275549</td>
<td>-533.126723</td>
<td>-534.074861</td>
<td>-533.600943</td>
<td>-534.057343</td>
<td>22343.43</td>
<td>9142.86</td>
<td>455.98</td>
</tr>
<tr>
<th>2</th>
<td>-529.378333</td>
<td>-529.036879</td>
<td>-528.462048</td>
<td>-529.954048</td>
<td>-525.608291</td>
<td>-527.275759</td>
<td>-528.747860</td>
<td>-527.931698</td>
<td>-527.359286</td>
<td>-526.207417</td>
<td>16686.87</td>
<td>4285.71</td>
<td>389.34</td>
</tr>
<tr>
<th>3</th>
<td>-533.691969</td>
<td>-531.952383</td>
<td>-534.055335</td>
<td>-532.632034</td>
<td>-533.845509</td>
<td>-532.641484</td>
<td>-533.093253</td>
<td>-534.387168</td>
<td>-533.574950</td>
<td>-533.370863</td>
<td>21494.95</td>
<td>8000.00</td>
<td>424.80</td>
</tr>
<tr>
<th>4</th>
<td>-533.735260</td>
<td>-532.092981</td>
<td>-534.197462</td>
<td>-532.649674</td>
<td>-533.692935</td>
<td>-532.890857</td>
<td>-533.219227</td>
<td>-534.262520</td>
<td>-533.789781</td>
<td>-533.422763</td>
<td>21494.95</td>
<td>11428.57</td>
<td>469.85</td>
</tr>
</tbody>
</table>
</div>
</div>
<br />
<br />
.. GENERATED FROM PYTHON SOURCE LINES 204-205
This can be saved as usual with `df.to_csv('pathname')` using Pandas. For the lithological block model, one good option is to save it as a numpy array, using `numpy.save()`.
.. GENERATED FROM PYTHON SOURCE LINES 205-208
.. code-block:: default
np.save('../../data/outputs/MCexample_10realizations.npy', lith_blocks)
.. GENERATED FROM PYTHON SOURCE LINES 209-212
Quick model analysis
--------------------
Let's have a quick first look at the resulting gravity and lithological block models. From the gravity dictionary, we can quickly generate a dataframe, convenient for further model analysis.
.. GENERATED FROM PYTHON SOURCE LINES 212-217
.. code-block:: default
prob_block = gp.bayesian.fields.probability(lith_blocks)
ie_block = gp.bayesian.fields.information_entropy(prob_block)
.. GENERATED FROM PYTHON SOURCE LINES 218-219
The following plot shows the probability of unit 5 in the probability block. With faults not being excluded, and counting of units starting with 0, we can see that the index 5 relates to the `Lower-filling` surface. The plot shows where to expect the unit. Everywhere, this unit is present throughout the simulations, the probability plot shows a bright yellow (probability = 1). Where it is always absent, we see the dark violet (probability = 0). The blueish-greenish areas are in between, meaning that in some realizations, the `Lower-filling` unit is present there, in other realization it is not.
.. GENERATED FROM PYTHON SOURCE LINES 219-228
.. code-block:: default
layer = 5
gp.plot_2d(geo_model,
show_lith=False, show_boundaries=False, show_data=False,
regular_grid=prob_block[layer],
kwargs_regular_grid={'cmap': 'viridis',
'norm': None}
);
.. image:: /Tutorials/images/sphx_glr_02_POC_create-MC-ensemble_003.png
:alt: Cell Number: mid Direction: y
:class: sphx-glr-single-img
.. rst-class:: sphx-glr-script-out
Out:
.. code-block:: none
<gempy.plot.visualization_2d.Plot2D object at 0x0000017586724700>
.. GENERATED FROM PYTHON SOURCE LINES 229-232
In the for-loop above, we not only varied the bottom boundary of the `Lower-filling` unit, but also `Upper-filling` and `Unconformity`. Using the measure of information entropy, we can visualize the parts of the model, where the most change is happening, i.e. where entropy is largest. Black areas in the following plot have zero information entropy, as there is only one "microstate" for the system, i.e. the model ensemble.
This means, we'd always encounter the same unit at the same place in every ensemble member. The colored areas, however, are areas where we'd encounter different geological units between ensemble members.
.. GENERATED FROM PYTHON SOURCE LINES 232-240
.. code-block:: default
gp.plot_2d(geo_model,
show_lith=False, show_boundaries=False, show_data=False,
regular_grid=ie_block,
kwargs_regular_grid={'cmap': 'magma',
'norm': None}
);
.. image:: /Tutorials/images/sphx_glr_02_POC_create-MC-ensemble_004.png
:alt: Cell Number: mid Direction: y
:class: sphx-glr-single-img
.. rst-class:: sphx-glr-script-out
Out:
.. code-block:: none
<gempy.plot.visualization_2d.Plot2D object at 0x0000017581AF3F10>
.. GENERATED FROM PYTHON SOURCE LINES 241-242
Finally, let's have a look at the gravity. We'll simply have a look at mean and standard deviation of the simulated gravity of the ensemble:
.. GENERATED FROM PYTHON SOURCE LINES 242-269
.. code-block:: default
# make subplots with mean and std
gravdf_plt = pd.DataFrame.from_dict(grav)
gravdf_plt.to_csv('../../data/outputs/MCexample_10grav.csv', index=False)
fig, axs = plt.subplots(1,2, figsize=[15,5], sharey=True)
m_grav = np.mean(gravdf_plt, axis=1)
st_grav = np.std(gravdf_plt, axis=1)
m = axs[0].scatter(grav_stations['X'], grav_stations['Y'], c=m_grav,
marker='s', s=90, cmap='magma', zorder=2)
axs[0].contourf(dtm[:,:,0], dtm[:,:,1], dtm[:,:,2],20, cmap='gist_earth', zorder=0)
axs[0].contour(dtm[:,:,0], dtm[:,:,1], dtm[:,:,2],10, colors='gray', zorder=1)
s = axs[1].scatter(grav_stations['X'], grav_stations['Y'], c=st_grav,
marker='s', s=90, cmap='magma', zorder=2)
axs[1].contourf(dtm[:,:,0], dtm[:,:,1], dtm[:,:,2],20, cmap='gist_earth', zorder=0)
axs[1].contour(dtm[:,:,0], dtm[:,:,1], dtm[:,:,2],10, colors='gray', zorder=1)
fig.colorbar(m, ax=axs[0], label='gravity')
fig.colorbar(s, ax=axs[1], label='std of gravity')
axs[0].set_title('Ensemble mean')
axs[1].set_title('Ensemble standard deviation')
axs[0].set_ylabel('Y [m]')
axs[0].set_xlabel('X [m]')
axs[1].set_xlabel('X [m]')
fig.tight_layout()
.. image:: /Tutorials/images/sphx_glr_02_POC_create-MC-ensemble_005.png
:alt: Ensemble mean, Ensemble standard deviation
:class: sphx-glr-single-img
.. rst-class:: sphx-glr-timing
**Total running time of the script:** ( 4 minutes 54.026 seconds)
.. _sphx_glr_download_Tutorials_02_POC_create-MC-ensemble.py:
.. only :: html
.. container:: sphx-glr-footer
:class: sphx-glr-footer-example
.. container:: sphx-glr-download sphx-glr-download-python
:download:`Download Python source code: 02_POC_create-MC-ensemble.py <02_POC_create-MC-ensemble.py>`
.. container:: sphx-glr-download sphx-glr-download-jupyter
:download:`Download Jupyter notebook: 02_POC_create-MC-ensemble.ipynb <02_POC_create-MC-ensemble.ipynb>`
.. only:: html
.. rst-class:: sphx-glr-signature
`Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_