Cómo filtrar una curva utilizando scipy

La resolución vertical de las curvas de las diferentes herramientas de perfilaje suelen estar pensadas para brindarnos un óptimo análisis del pozo de unos pocos metros, incluso centímetros. Cuando estamos haciendo trabajos donde debemos examinar cientos o miles de metros del pozo, posiblemente una diferente resolución vertical de las curvas nos ayude a diferenciar eventos macros o de mayor escala, o simplemente nos elimine un poco de ruido.

Bienvenidas/os!!, muchas gracias por visitarnos, si están comenzando en Python y este es el primer artículo al que entran en nuestro sitio, les recomendamos visitar Primeros Pasos.

ℹ️ Recuerden que todos los ejemplos de códigos en nuestro sitio web van a encontrarlos en nuestro repositorio GitHub. Los códigos de este artículo los encontrarán en curve_smooth.ipynb.

En este artículo veremos los pasos para filtrar una sola curva con el filtro savgol (Savitzky-Golay) en scipy, y plotearla junto a la curva filtrada utilizando matplotlib.

Comencemos!!

IMPORTAR LIBRERIAS

import lasio
import numpy as np
import pandas as pd

from scipy.signal import savgol_filter

import matplotlib.pyplot as plt
%matplotlib inline

CARGA DE LOS DATOS

La carga de los datos la haremos utilizando la librería lasio en primera instancia, y luego la convertiremos a un dataframe pandas.

las = lasio.read(r"D:\Pablo\Python\working_files\Equinor\15_9-23.las")
well = las.df()
well = well.reindex(sorted(well.columns), axis=1)  # Ordenamos las columnas alfabéticamente
pd.set_option('display.max_columns', 200)          # Muestra las primeras 200 columnas
well

FILTRADO Y PLOTEO

Ahora nos toca elegir la curva que queremos filtrar, aplicarle el filtro y ver ambas en un plot.

Para simplificar las cosas, vamos a utilizar variables en lugar de valores fijos para ciertas cosas (por ejemplo, el nombre de la curva a filtrar será crv_original, las escalas del plot serán crv_izq y crv_der, etc). Esto lo hacemos para no tener que estar modificando TODO el código cada vez que elijamos parámetros y/o curvas diferentes.

En este caso, debido a la longitud del código, vamos a tratar que no sea tan abrumador explicándolo paso a paso.


Elegir la curva original y darle un nombre a la filtrada

### CURVAS
crv_original = 'GR'
crv_filtered = 'GR_F'

Acá vamos a elegir la curva del dataframe well que queremos filtrar y pondremos su nombre a la variable crv_original. También elegiremos el nombre de la curva resultante de la aplicación del filtro. (crv_filtered)

⚠️ Las comillas deben estar encerrando el nombre de la curva para que funcione el código de abajo. También recuerden que python diferencia entre mayúsculas y minúsculas.

Escalas para graficar las curvas en el plot

# Escala curvas
crv_izq = 0     
crv_der = 300

Seguidamente elegimos las escalas en las cuales queremos que se grafiquen las curvas

ℹ️ Notemos que en lugar de caracterizar las escalas utilizando min y max, utilizamos izq y der. Esto se debe a que, a veces, ploteamos curvas en sentido inverso, con el valor mínimo a la derecha, y el máximo a la izquierda (por ejemplo: porosidades, Delta-T, etc), por ello creemos que es más conveniente definir las escalas utilizando izq-der.

Parámetros del filtro Savitzky-Golay

### FILTRO
savgol_wl = 355  # Cantidad de muestras para filtrar (debe ser impar)
savgol_pol = 3   # orden del polinomio

En este caso sólo utilizamos 2 parámetros, windows length (savgol_wl) y polyorder (savgol_pol). Si desean experimentar las opciones de este filtro, visiten https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.savgol_filter.html.

⚠️ savgol_wl debe ser un número impar.
⚠️ savgol_pol debe ser siempre menor que savgol_wl. 
#### Top-Bottom del plot

top_plot = las.well.STRT.value  # Valor default TOP del plot (modificar si es necesario)
bot_plot = las.well.STOP.value  # Valor default BOTTOM del plot (modificar si es necesario)

Por último, elegiremos el intervalo que queremos graficar. Por default, están los valores top y bot del las que importamos con lasio.

En nuestra opinión, el resto del código no necesita cambios substanciales para que podamos trabajar correctamente.

ℹ️ Pueden definir más variables para hacer más dinámico el código si lo desean (color de las curvas por ejemplo, etc).

⚠️ Para que todo funcione correctamente, el objeto lasio debe llamarse las, y el dataframe well. Si los modificaron en los pasos anteriores, deben hacer las modificaciones pertinentes en este código.

⚠️ El código de esta celda va a crear un dataframe auxiliar llamado well2 cuando trabaje con el filtro savgol. Por favor no utilicen este nombre para nombrar al dataframe original, porque sobrescribirá sus datos. 
############################ VALORES DE ENTRADA ############################
### CURVAS
crv_original = 'GR'
crv_filtered = 'GR_F'

# Escala curvas
crv_izq = 0     # a veces ploteamos en sentido inverso (por ej: porosidades)
crv_der = 300   # por eso es preferible definir izq-der en lugar de min-max 

### FILTRO
savgol_wl = 355  # Cantidad de muestras para filtrar (debe ser impar)
savgol_pol = 3   # polinomio


## TOP-BOT plot
top_plot = las.well.STRT.value  # Valor default TOP del plot (modificar si es necesario)
bot_plot = las.well.STOP.value  # Valor default BOTTOM del plot (modificar si es necesario)


############################ CODIGO ######################################################################

#### Filtrar curva

well2=pd.DataFrame(well[well[crv_original].notnull()][crv_original])             # Tuve que crear otro df porque hay Nan en GR
well2[crv_filtered] = savgol_filter(well2[crv_original], savgol_wl, savgol_pol)  # (Para aplicar el filtro tengo que sacar 
well[crv_filtered]=well2[crv_filtered]                                           # los valores nulos del calculo)



### PLOT #####

fig, ax = plt.subplots(figsize=(50,10))#Set up the plot axes

ax1 = plt.subplot2grid((1,6), (0,0), rowspan=1, colspan = 1) 
ax2 = ax1.twiny()

ax1.plot(well[crv_original], well.index, color = "#53B923", linewidth = 1)
ax1.set_ylim(bot_plot, top_plot)
ax1.set_xlabel(f"Curva original ({crv_original})")
ax1.xaxis.label.set_color("green")
ax1.set_xlim(crv_izq, crv_der)
ax1.tick_params(axis='x', colors="green")
ax1.spines["top"].set_position(("axes", 1.0))
ax1.spines["top"].set_edgecolor("green")
ax1.title.set_color('green')
ax1.xaxis.set_ticks_position("top")
ax1.xaxis.set_label_position("top")


ax2.plot(well[crv_filtered], well.index, color = "black", linewidth = 2)
ax2.set_xlabel(f"Curva filtrada ({crv_filtered})")
ax2.xaxis.label.set_color("black")
ax2.set_xlim(crv_izq, crv_der)
ax2.tick_params(axis='x', colors="black")
ax2.spines["top"].set_position(("axes", 1.08))
ax2.spines["top"].set_edgecolor("black")
ax2.title.set_color('black')
ax2.xaxis.set_ticks_position("top")
ax2.xaxis.set_label_position("top")


plt.tight_layout()
plt.show()

En resumen, cargamos los datos de un archivo .las y elegimos una de sus curvas. Seguidamente la filtramos utilizando el filtro Savitzky-Golay y graficamos ambas en un plot.

Les agradecemos su tiempo y esperamos fervientemente que hayan disfrutado este artículo. Si tienen alguna consulta, desean hacer algún comentario o sugerencia para mejorar el contenido, o simplemente indicarles qué les pareció este artículo, debajo pueden hacerlo.

Esperamos reencontrarlos en algún otro artículo del sitio. Hasta luego!


comments powered by Disqus