Gempy es una librería de código abierto de Python para generar modelos geológicos estructurales 3D completos. La librería es un desarrollo completo para crear modelos geológicos de interfases, fallas y orientaciones de capas, también relaciona la secuencia de las capas geológicas para representar intrusiones de rocas y órdenes de fallas.
El algoritmo para modelamiento geológico basado en la interpolación universal “cokriging” con el soporte de librerías matemáticas avanzadas de Python como Numpy, PyMC3 y Theano. Gempy crea un modelo grillado que puede ser visualizado en secciones 2D con Matplotlib o como geometrías 3d en VTK que permitirán la representación del modelo geológico en Paraview para divisiones, filtrado, transparencias y estilo.
Este tutorial es un ejemplo básico de una configuración geológica estratificada con 5 capas y una falla. Con finalidad de hacer que el tutorial sea accesible a la mayoría de usuarios, hemos creado un tutorial complementario acerca de Gempy en Windows con una distribución de repositorio de Anaconda.
Tutorial
Parte 1: Instalación de Gempy en Windows
Part 2: Ejemplo de modelo geológico con Gempy
Código en Python
Este es el código en Pytho de este tutorial:
Configurando el entorno de Python
En esta parte importamos las librerías requeridas para este tutorial. El código requiere Gempy tanto como Numpy y Matplotlib. Configuramos la opción de Jupyter para representación interactiva de gráficos de Matplotlib después de las celdas de código (%matplotlib inline).
Fijarse que las advertencias son únicamente mensajes que el usuario tiene que tomar en cuenta cuando corre el código, no significan un fallo en el código. Como este tutorial es Windows algunas librerías complementarias no pudieron ser instaladas, pero el rendimiento total del código de modelamiento geológico es completo.
# Embedding matplotlib figures in the notebooks
%matplotlib inline
# Importing GemPy
import gempy as gp
# Importing auxiliary libraries
import numpy as np
import matplotlib.pyplot as plt
WARNING (theano.configdefaults): g++ not available, if using conda: `conda install m2w64-toolchain`
E:\Software\Anaconda36\lib\site-packages\theano\configdefaults.py:560: UserWarning: DeprecationWarning: there is no c++ compiler.This is deprecated and with Theano 0.11 a c++ compiler will be mandatory
warnings.warn("DeprecationWarning: there is no c++ compiler."
WARNING (theano.configdefaults): g++ not detected ! Theano will be unable to execute optimized C-implementations (for both CPU and GPU) and will default to Python implementations. Performance will be severely degraded. To remove this warning, set Theano flags cxx to an empty string.
WARNING (theano.tensor.blas): Using NumPy C-API based implementation for BLAS functions.
E:\Software\Anaconda36\lib\site-packages\gempy\posterior_analysis.py:31: UserWarning: tqdm package not installed. No support for dynamic progress bars.
warnings.warn("tqdm package not installed. No support for dynamic progress bars.")Creación del objeto del modelo geológico y definición estratigráfica
El tutorial crea una grilla de 100 columnas x 100 filas x 100 capas sobre una extensión de 2km x 2km x 2km. Resoluciones mayores son posibles pero el tiempo computacional es más alto. El sistema de referencia son coordenadas locales, siguientes tutoriales evaluarán el rendimiento de Gempy con coordenadas UTM.
La orientación y contactos geológicos son importados de archivos CSV y convertidos en DataFrames de Pandas. Luego las series geológicas (fallas / formaciones) son definidas tanto como las secuencias de formación geológica.
Es importante mencionar que las fallas tienen que ser insertadas independientemente, donde el más joven es la primera entrada.
# Importing the data from CSV-files and setting extent and resolution
geo_data = gp.create_data([0,2000,0,2000,0,2000],[100,100,100],
path_o = "../Txt/simple_fault_model_orientations.csv", # importing orientation (foliation) data
path_i = "../Txt/simple_fault_model_points.csv") # importing point-positional interface data
gp.get_data(geo_data).loc[:,['X','Y','Z','formation']].head()| X | Y | Z | formation | ||
|---|---|---|---|---|---|
| interfaces | 52 | 700.0 | 1000.0 | 900.0 | Main_Fault |
| 53 | 600.0 | 1000.0 | 600.0 | Main_Fault | |
| 54 | 500.0 | 1000.0 | 300.0 | Main_Fault | |
| 55 | 800.0 | 1000.0 | 1200.0 | Main_Fault | |
| 56 | 900.0 | 1000.0 | 1500.0 | Main_Fault |
# Assigning series to formations as well as their order (timewise)
gp.set_series(geo_data, {"Fault_Series":'Main_Fault',
"Strat_Series": ('Sandstone_2','Siltstone', 'Shale', 'Sandstone_1')},
order_series = ["Fault_Series", 'Strat_Series'],
order_formations=['Main_Fault',
'Sandstone_2','Siltstone', 'Shale', 'Sandstone_1',
], verbose=0)Gráfico de secuencias geológicas
Gempy tiene algunas cualidades útiles para representar las series geológicas definidas y las secuencias de formaciones.
gp.get_sequential_pile(geo_data)
<gempy.plotting.sequential_pile.StratigraphicPile at 0x5a58400>
Revisión de la información de entrada
Los diferentes conjuntos de datos para la construcción del modelo geológico pueden ser accedidas a través de las funciones de Gempy ".get_"
# Review of the centroid coordinates from the model grid
gp.get_grid(geo_data).valuesarray([[ 10., 10., 10.],
[ 10., 10., 30.],
[ 10., 10., 50.],
...,
[ 1990., 1990., 1950.],
[ 1990., 1990., 1970.],
[ 1990., 1990., 1990.]], dtype=float32)# Defined interfases from the input CSV data
gp.get_data(geo_data, 'interfaces').loc[:,['X','Y','Z','formation']].head()| X | Y | Z | formation | |
|---|---|---|---|---|
| 52 | 700.0 | 1000.0 | 900.0 | Main_Fault |
| 53 | 600.0 | 1000.0 | 600.0 | Main_Fault |
| 54 | 500.0 | 1000.0 | 300.0 | Main_Fault |
| 55 | 800.0 | 1000.0 | 1200.0 | Main_Fault |
| 56 | 900.0 | 1000.0 | 1500.0 | Main_Fault |
# Defined layer orientations from the input CSV data
gp.get_data(geo_data, 'orientations').loc[:,['X','Y','Z','formation','azimuth']]| X | Y | Z | formation | azimuth | |
|---|---|---|---|---|---|
| 2 | 500 | 1000 | 864.602 | Main_Fault | 270 |
| 1 | 400 | 1000 | 1400.000 | Sandstone_2 | 90 |
| 0 | 1000 | 1000 | 950.000 | Shale | 90 |
Representación gráfica de la información de entrada
En esta parte las representaciones 2D y 3D fueron realizadas para presentar las interfaces y orientaciones.
gp.plot_data(geo_data, direction='y')E:\Software\Anaconda36\lib\site-packages\gempy\gempy_front.py:927: FutureWarning: gempy plotting functionality will be moved in version 1.2, use gempy.plotting module instead
warnings.warn("gempy plotting functionality will be moved in version 1.2, use gempy.plotting module instead", FutureWarning)
gp.plotting.plot_data_3D(geo_data)
<gempy.plotting.plot.vtkPlot at 0x10879630>
Interpolación Geológica
Una vez la información de entrada está lista, podemos definir la información y parámetros para la interpolación con el método “interpolationData” de la librería de Gempy El modelo geológico es calculado bajo el método "compute_model". Resultados del proceso del modelo son la litología y fallas en la misma dimensión del array como el “geo_data” The geological model is calculated under the "compute_model" method. Results from the model process are the lithology and faults on the same array dimensions as the geo_data.
interp_data = gp.InterpolatorData(geo_data, u_grade=[1,1], output='geology', compile_theano=True, theano_optimizer='fast_compile')
Compiling theano function...
Compilation Done!
Level of Optimization: fast_compile
Device: cpu
Precision: float32
Number of faults: 1
interp_data.geo_data_res.formations.as_matrix
<bound method NDFrame.as_matrix of value formation_number
Main_Fault 1 1
Sandstone_2 2 2
Siltstone 3 3
Shale 4 4
Sandstone_1 5 5
basement 6 6>
interp_data.geo_data_res.get_formations()
[Main_Fault, Sandstone_2, Siltstone, Shale, Sandstone_1]
Categories (5, object): [Main_Fault, Sandstone_2, Siltstone, Shale, Sandstone_1]
lith_block, fault_block = gp.compute_model(interp_data)Exploración del modelo litológico
El bloque litológico tiene 2 partes, la primera tiene información acerca de la formación litológica, mientras que el segundo representa la orientación. En esta parte la distribución de la litología y la información de fallas separadas son representados como histogramas.
lith_block[0]
array([ 6.3131361 , 6.24877167, 6.19397354, ..., 2.00398016,
2.00626612, 2.00983 ], dtype=float32)
plt.hist(lith_block[0],bins=100)
plt.show()
plt.hist(fault_block[0],bins=10)
plt.show()
Geological model representation
Como cualquier array de Numpy, el resultado de los bloques litológicos puede ser representado en Matplotlib. Sin embargo, Gempy tiene métodos espaciales para la representación de la sección de corte. Con el uso de los complementos de Jupyter una representación de las secciones geológicas en dirección al eje Y es desarrollado con un certificación con una manilla para cambiar filas.
gp.plotting.plot_section(geo_data, lith_block[0], cell_number=50, direction='y', plot_data=False)
import ipywidgets as widgets
def plotCrossSection(cell):
gp.plotting.plot_section(geo_data, lith_block[0], cell_number=cell, direction='y', plot_data=False)
widgets.interact(plotCrossSection, cell=widgets.IntSlider(min=0,max=99,step=1,value=50) )
If you're reading this message in the Jupyter Notebook or JupyterLab Notebook, it may mean that the widgets JavaScript is still loading. If this message persists, it likely means that the widgets JavaScript library is either not installed or not enabled. See the Jupyter Widgets Documentation for setup instructions.
If you're reading this message in another frontend (for example, a static rendering on GitHub or NBViewer), it may mean that your frontend doesn't currently support widgets.
<function __main__.plotCrossSection>
gp.plotting.plot_scalar_field(geo_data, lith_block[1], cell_number=50, N=6,
direction='y', plot_data=False)
plt.colorbar()
plt.show()
ver_s, sim_s = gp.get_surfaces(interp_data,lith_block[1],
fault_block[1],
original_scale=True)
gp.plotting.plot_surfaces_3D(geo_data, ver_s, sim_s)
<gempy.plotting.plot.vtkPlot at 0x1f327550>
Datos de entrada
Los datos de entrada para este tutorial se pueden descargar de este enlace.
