KRIGING SIMPLE 2D CON GSLIB REIMPLEMENTADA EN PYTHON

HEBER HERNÁNDEZ G. – NUBE MINERA (MAYO 2019)

GSLIB es el acrónimo de “Geoestatistical Software Library”, el cual es uno de los primeros directorios contenedores de software geoestadístico en el mundo, desarrollado en el departamento de ingeniería del Petroleo de la Universidad de Stanford.

Su primera publicación fue en el año 1992, con una amplia colección de algoritmos para realizar análisis estructural, estimación y simulación geoestadística en dominios espaciales bidimensionales y tridimensionales.

Algunos investigadores lo consideran como el principal desarrollo informático en la historia de la Geoestadística (Remy, 2009).

Debido a su diseño modular y su gran eficiencia ha sido empleada como herramienta de enseñanza e investigación en el desarrollo de nuevos algoritmos y metodologías (Deutsch y Journel, 1998; Remy, 2009).

Su principal desventaja son las limitadas capacidades de visualización de resultados, puesto que solo soporta la visualización mediante archivos Post-Script, además de la inexistencia de una interfaz gráfica freeware de usuario (Wingle et al., 1999).

De lo anterior, es que dada la re-implementación de GSLIB en Python, específicamente en la biblioteca GeostatsPy, por el Dr. Michael J. Pyrcz (https://github.com/GeostatsGuy/GeostatsPy), las limitaciones existentes en cuanto a visualización y tiempos de rutinas, que hace 20 años atrás eran un verdadero problema, ya han sido resueltos.


En esta publicación de Nube Minera, se realiza una demostración de la biblioteca GeostatsPy con una estimación a nivel 2D con Kriging Simple.


La demostración requiere tener algunas bibliotecas adicionales a GeostatsPy:

import geostatspy.GSLIB as GSLIB # Utilidades GSLIB y funciones envueltas.
import geostatspy.geostats as geostats # GSLIB convertido a Python.
import os # Para configurar el directorio de trabajo.
import numpy as np # Biblioteca de funciones matemáticas.
import pandas as pd # Biblioteca para manipulación de marcos de datos.
import matplotlib.pyplot as plt # Biblioteca para salida gráfica.

Con Pandas se procede a llamar al archivo (descargar aquí) que contiene las muestras 2D:

ruta ="C:/Users/Heber/Desktop/Demostracion/" #Carpeta en la que se encuentra el archivo de muestras.
archivo ="muestras_ks.csv" #Nombre del archivo de muestras.

Con OS, se debe establecer la dirección de los ejecutables de GSLIB, los cuales se descargan desde la pagina oficial del proyecto (http://www.statios.com/Quick/gslib.html):

os.chdir("C:/Users/Heber/Desktop/Demostracion/GSLIB_FOLDER") # ruta de ejecutables de GSLIB.

Con la siguiente función de Pandas, se lee el archivo de muestras, y ya con esto se pueden generar múltiples análisis con los datos:

df = pd.read_csv(ruta + archivo, sep=";") #lectura de las muestras, las cuales se encuentran separadas por punto y coma.

El primer método que utilizaremos de GeostatsPy es GSLIB.locmap(), con el cual se desplegara un mapa de ubicación para las muestras.

Si le interesa conocer mas acerca de este ejecutable original de GSLIB puede consultar en: http://www.gslib.com/gslib_help/locmap.html

Dado que no es el propósito realizar un análisis estadístico exhaustivo sobre los datos, solo se resumirá que la variable de interés presenta las siguientes características:

  • N° de muestras: 144
  • Media: 4.475617
  • Desviación estándar: 1.988807
  • Mínimo: 0.203012
  • Máximo: 8.980486

Esto ultimo se puede obtener utilizando la función «describe» de Pandas sobre los datos:

df.describe()

A continuación se presenta  un histograma ploteado con Matplotlib que denota una leve asimetría positiva con un factor de 0.36.

plt.hist(df['muestra'], bins="sturges",color='skyblue', histtype='stepfilled')
plt.grid(False)
plt.xlabel("Variable")
plt.ylabel("Frecuencia")
plt.title("Histograma muestras")
plt.grid(True)

Para el análisis de continuidad espacial, se ha utilizado el semivariograma experimental para dos direcciones: 0° y 90° azimut. Luego también se ha incluido un semivariograma omnidireccional en el mismo gráfico.

Desde GeostatPy se utiliza GSLIB.gamv_2d() en este propósito. Mas información sobre este ejecutable en: http://www.statios.com/Help/gamv.html 

Luego con GSLIB.make_variogram()  y sobre el semivariograma omnidireccional, se ajusta un único modelo elemental del tipo esférico con alcance de 300 metros y meseta 1.2.

Finalmente con el modelo estructural y una media conocida asumida de 4.45 para la variable, es que se utiliza el estimador de Kriging Simple re implementado en Python como geostats.kb2d(). 

Puede consultar el ejecutable original de la GSLIB en: http://www.statios.com/Help/kb2d.html

Para la visualización ocupamos GSLIB.pixelplt_st() de GeostatsPy:

En conclusión, la biblioteca GeostatsPy entre muchas otras funcionalidades, puede generar modelos 2D de estimación con Kriging Simple de manera rápida y con buenos resultados exportables a archivos CSV. 

Si es de su interés aprender un poco de Python enfocado a Geoestadística Lineal y familiarizarse con los códigos no expuestos en esta publicación que hacen relación a la creación del semivariograma y la estimación, lo invito a visitar el siguiente ENLACE. 

junio 16, 2019
top
Nube Minera © 2024