Modelamiento Geológico Estructural 3D en Python con Gempy - Tutorial

3DStructuralGeologicalModelinginPythonwithGempy.png

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>
output_7_1.png

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).values
array([[   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)
output_13_1.png
gp.plotting.plot_data_3D(geo_data)

<gempy.plotting.plot.vtkPlot at 0x10879630>
output_27_1.png

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()
output_22_0.png
plt.hist(fault_block[0],bins=10)
plt.show()
output_23_0.png

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)
output_25_0.png
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) )
descarga.png

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()
output_27_0.png
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.

Smiley face

Suscríbase a nuestro boletín gratuito para recibir noticias, datos interesantes y fechas de nuestros cursos en recursos hídricos.

Posted on April 11, 2019 and filed under TutorialModflow, TutorialPython.