Cómo rellenar el área entre 2 curvas con un color o una paleta de colores en matplotlib

El análisis de perfiles de pozo se apoya mucho en gráficos para poder identificar características petrofísicas generales de los reservorios que atraviesa el pozo tan sólo observando los logs. Para ello, se utilizan varias herramientas visuales, de las cuales el relleno de áreas entre curvas es la herramienta que trataremos en este artículo.

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 matplotlib_cmap_fillbetweenx.ipynb.

Como comentamos al comienzo del artículo, el análisis de perfiles de pozo utiliza mucho los logs para obtener una interpretación petrofísica inicial del pozo de manera visual. El relleno de áreas se usa cuando queremos resaltar cruces de ciertas curvas (densidad-neutron, caliper-bitsize), o cambios en los valores de una curva (gamma ray).

En esta oportunidad, discutiremos la metodología para rellenar el área entre 2 curvas, o entre un valor del track y una curva; con un color, o una paleta de colores utilizando matplotlib.

Si les parece bien, empecemos:

IMPORTAR LIBRERIAS

import numpy as np
import pandas as pd
import lasio

import matplotlib.pyplot as plt
from matplotlib import cm
%matplotlib inline

Como siempre, lo primero que haremos es importar todas las librerías que utilizaremos en el notebook.


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()
pd.set_option('display.max_columns', None)           # Muestra todas las columnas
well = well.reindex(sorted(well.columns), axis=1)    # Ordena las columnas alfabeticamente
well

fill_betweenx() en matplotlib

Para rellenar el área entre curvas utilizando matplotlib se utilizan las funciones fill_between(), (si queremos rellenar entre valores distintos del eje y) y fill_betweenx() (si el relleno es con valores diferentes ubicados en el eje x).

Nosotros, en los gráficos símil logs, utilizaremos fill_betweenx() .

Algo importante a tener en cuenta, es que fill_betweenx() necesita que las 2 curvas estén en la misma escala para poder rellenar correctamente el área sin tener que hacer transformaciones matemáticas.

En este caso en particular, vamos a rellenar el cruce entre las curvas de densidad y neutrón.

Entonces, el siguiente paso es obtener la porosidad de la densidad a partir de la curva de densidad bulk (RHOB), ya que la curva de porosidad de la densidad no está en el archivo .las.

ℹ️ Calcularemos 2 porosidades de densidad, una para martriz arenisca (DPHISS) y otra para matriz calcárea (DPHILS). Pero utilizaremos solamente una en el ejemplo (DPHILS). 
def den_por_sand(bulk):
    return (2.65-bulk)/(2.65-1)    

def den_por_lime(bulk):
    return (2.71-bulk)/(2.71-1)    

well['DPHISS']= well[['RHOB']].apply(den_por_sand)
well['DPHILS']= well[['RHOB']].apply(den_por_lime)

Seguidamente, ubicaremos las curvas de porosidad neutrón y porosidad densidad en un mismo track y pintaremos las áreas bajo las curvas de la siguiente manera:

  • Relleno golor gris: cuando la por_neu > por_den (generalmente una alta porosidad neutrónica con respecto a la porosidad de la densidad indica presencia de arcillas)
  • Relleno color amarillo: cuando la por_den > por_neu (generalmente asociado a reservorios limpios o también a presencia de Hidrocarburos, etc)
ax3.fill_betweenx(well.index, well[crv_neu_por], well[crv_den_por], where=well[crv_neu_por] > well[crv_den_por], interpolate = True, color="grey",alpha=0.3)

ax3.fill_betweenx(well.index, well[crv_neu_por], well[crv_den_por], where=well[crv_neu_por] < well[crv_den_por], interpolate = True, color="yellow",alpha=0.9)

Donde:

  • well.index Curva del eje y
  • well[crv_neu_por], well[crv_den_por] Las 2 curvas en el eje x (si en lugar de 2 curvas, utilizamos una y un valor numérico, rellena el área entre ese valor y la curva).
  • where lo utilizamos, si es necesario, para declarar donde queremos que se rellene (si las curvas se cruzan por ejemplo, que es nuestro caso).
  • interpolate útil cuando las curvas se cruzan y utilizamos where (nuestro caso).
  • color color del relleno.
  • alpha claridad del color, (entre 0 y 1), mientras más cerca de 0, el color es más claro.

ℹ️ Si desean conocer más sobre esta función, aquí está el link https://matplotlib.org/stable/api/_as_gen/matplotlib.axes.Axes.fill_betweenx.html .

El código completo del gráfico es el siguiente:

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

crv_den_por = 'DPHILS'
crv_neu_por = 'NPHI'

por_izq = .45
por_der = -.15

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

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

ax3 = plt.subplot2grid((1,5), (0,0), rowspan=1, colspan = 1) 
ax31 = ax3.twiny()


###############################################################################
### DEN POROSITY ###

ax3.set_ylim(bot_plot, top_plot)
ax3.plot(well[crv_den_por], well.index, color = "red", linewidth = .2)
ax3.set_xlabel(f"Density Porosity ({well[crv_den_por].name})")
ax3.set_xlim(por_izq,por_der)
ax3.xaxis.label.set_color("red")
ax3.xaxis.set_ticks_position("top")
ax3.xaxis.set_label_position("top")
ax3.spines["top"].set_position(("axes", 1.))
ax3.tick_params(axis='x', colors="r")
ax3.set_xticks([0.45,0.35,0.25,0.15,0.05,-0.05,-0.15])
ax3.spines["top"].set_edgecolor("r")
################################################################################
### NEU POROSITY ###

ax31.plot(well[crv_neu_por], well.index, color = "blue", linewidth =.2)
ax31.set_xlabel(f"Neutron Porosity ({well[crv_neu_por].name})")
ax31.xaxis.label.set_color("b")
ax31.set_xlim(por_izq,por_der)
ax31.tick_params(axis='x', colors="blue")
ax31.spines["top"].set_position(("axes", 1.07))
ax31.spines["top"].set_visible(True)
ax31.set_xticks([0.45,0.35,0.25,0.15,0.05,-0.05,-0.15])
ax31.spines["top"].set_edgecolor("blue")
################################################################################

ax3.fill_betweenx(well.index, well[crv_neu_por], well[crv_den_por], where=well[crv_neu_por] > well[crv_den_por], interpolate = True, color="grey",alpha=0.3)
ax3.fill_betweenx(well.index, well[crv_neu_por], well[crv_den_por], where=well[crv_neu_por] < well[crv_den_por], interpolate = True, color="yellow",alpha=0.9)

plt.tight_layout()
plt.show()

colormaps (cmaps) en matplotlib

Un mapa de colores (a menudo denominado tabla de colores o paleta) es una matriz de colores que se utiliza para asignar datos de píxeles (representados como índices en la tabla de colores) a los valores de color reales.

La idea de elegir un buen mapa de colores, es encontrar una buena representación de los datos con los que estemos trabajando mediante una gama de colores 3D.

Los colormaps se clasifican en secuenciales, divergentes, cíclicos, cualitativos, etc, según ciertas características (entre ellas, cambios en la luminosidad, saturación, etc).

Los mapas de colores se utilizan mucho en scatter plots para pintar de diferentes colores los puntos según su valor (Ejemplo). Sin embargo, en nuestro caso los vamos a usar junto con fill_betweenx() para que la paleta de colores se aplique al área entre una curva y el comienzo del track.

matplotlib cuenta con varios mapas de colores para su utilización, de los cuales aquí les dejaremos estos ejemplos.

ℹ️ Si desean más información sobre colormaps, por favor visiten: https://matplotlib.org/stable/tutorials/colors/colormaps.html .

Ahora les mostraremos una de las maneras de obtener un mapa de colores, utilizar la funcion get_cmap:

plt.get_cmap("nombre del colormap", "cantidad de elementos del cmap")

Seguidamente, veremos 3 ejemplos para ver cómo trabaja esta función:


Ejemplo 1:

plt.get_cmap('plasma', 5)

Ejemplo 2:

plt.get_cmap('plasma', 50)

En este ejemplo aumentamos la cantidad de elementos del cmap con respecto al ejemplo 1.


Ejemplo 3:

plt.get_cmap('plasma_r', 50)

Aquí le agregamos al final del nombre del cmap el término _r, lo cual invierte la paleta de colores.


Volviendo a nuestro ejemplo, y como ya explicamos anteriormente, los cambios en los valores de una sola curva pueden ser indicios de cambios litológicos en los reservorios que atraviesa el pozo (por ejemplo, muchas veces los valores de la curva gamma ray (GR) van a estar relacionados a la presencia o no de arcillosidad, mientras más arcillosidad, más altos valores de gamma ray).

Si podemos definir, o al menos aproximar, un valor de GR de 0 arcillosidad y de 100% de arcillosidad, entonces los valores intermedios van a tener una proporción de arcilla y otra de roca/matriz (arenisca, caliza, etc). Si esto lo podemos distinguir visualmente mediante el uso de diferentes colores, entonces nos facilitará la tarea.

Para tratar de lograr esto, vamos a utilizar conjuntamente cmap y fill_betweenx() en el código siguiente:

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


### CURVA: nombre, escala y color
crv_1 = 'GR'             # nombre de la curva
crv_izq = 0              # Limite izquierdo curva
crv_der = 300            # Limite derecho curva
crv_color = '#145A32'    # Color de la curva en el plot


#### CMAP (mapa de colores) ###
crv_cmap_min = 0      # Valor mínimo color de la curva para el cmap (arena/caliza/clean)
crv_cmap_max = 120    # Valor máximo color de la curva para el cmap (arcilla/shale)

cmap_len = 100        # cantidad de elementos del cmap y cantidad de elementos 
                        # entre los valores mín-máx de la curva a la que le aplico el cmap.

cmap_name = 'YlOrBr'  # Nombre del cmap que voy a utilizar

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

cmap = plt.get_cmap(cmap_name, cmap_len)

cmap_crv = np.linspace(crv_cmap_min, crv_cmap_max, cmap_len)


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

ax1 = plt.subplot2grid((1,4), (0,0), rowspan=1, colspan = 1) 

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

## CMAP y FILL_BETWEENX
for item in cmap_crv:
    color = cmap(item/crv_cmap_max) 
    ax1.fill_betweenx(well.index, crv_izq , well[crv_1], where = well[crv_1] >= item,  color = color)

plt.show()

Listo, ya hemos visto cómo utilizar estas herramientas de relleno de áreas. Para finalizar, vamos a hacer un poco de log analysis poniendo ambos tracks uno al lado del otro. De esta manera vamos a tener, posiblemente, una visión más completa de los reservorios que atraviesan nuestro pozo:

Código:

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


########## CURVAS ##########
crv_1 = 'GR'             # nombre de la curva
crv_izq = 0              # Limite izquierdo curva
crv_der = 300            # Limite derecho curva
crv_color = '#145A32'    # Color de la curva en el plot

crv_den_por = 'DPHILS'
crv_neu_por = 'NPHI'
por_izq = .45
por_der = -.15


########## CMAP ##########
crv_cmap_min = 0         # Valor mínimo color de la curva para el cmap (arena/clean)
crv_cmap_max = 120       # Valor máximo color de la curva para el cmap (arcilla/shale)

cmap_len = 100           # cantidad de elementos del cmap y cantidad de elementos 
                        # entre los valores mín-máx de la curva a la que le aplico el cmap.

cmap_name = 'YlOrBr'     # Nombre del cmap que voy a utilizar


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

cmap = plt.get_cmap(cmap_name, cmap_len)
cmap_crv = np.linspace(crv_cmap_min, crv_cmap_max, cmap_len)

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

ax1 = plt.subplot2grid((1,6), (0,0), rowspan=1, colspan = 1) 
ax3 = plt.subplot2grid((1,6), (0,1), rowspan=1, colspan = 1, sharey = ax1) 
ax31 = ax3.twiny()

############################################################
### GAMMA RAY ###

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

plt.setp(ax3.get_yticklabels(), visible = False)
ax3.plot(well[crv_den_por], well.index, color = "red", linewidth = .2)
ax3.set_xlabel(f"Density Porosity ({well[crv_den_por].name})")
ax3.set_xlim(por_izq,por_der)
ax3.xaxis.label.set_color("red")
ax3.xaxis.set_ticks_position("top")
ax3.xaxis.set_label_position("top")
ax3.spines["top"].set_position(("axes", 1.))
ax3.tick_params(axis='x', colors="r")
ax3.set_xticks([0.45,0.35,0.25,0.15,0.05,-0.05,-0.15])
ax3.spines["top"].set_edgecolor("r")

############################################################
### NEU POR ###

ax31.plot(well[crv_neu_por], well.index, color = "blue", linewidth =.2)
ax31.set_xlabel(f"Neutron Porosity ({well[crv_neu_por].name})")
ax31.xaxis.label.set_color("b")
ax31.set_xlim(por_izq,por_der)
ax31.tick_params(axis='x', colors="blue")
ax31.spines["top"].set_position(("axes", 1.07))
ax31.spines["top"].set_visible(True)
ax31.set_xticks([0.45,0.35,0.25,0.15,0.05,-0.05,-0.15])
ax31.spines["top"].set_edgecolor("blue")

############################################################
### FILL_BETWEENX y CMAP ###

ax3.fill_betweenx(well.index, well[crv_neu_por], well[crv_den_por], where=well[crv_neu_por] > well[crv_den_por], color="grey",alpha=0.4)
ax3.fill_betweenx(well.index, well[crv_neu_por], well[crv_den_por], where=well[crv_neu_por] < well[crv_den_por], color="yellow",alpha=0.4)

for item in cmap_crv:
    color = cmap(item/crv_cmap_max) 
    ax1.fill_betweenx(well.index, crv_izq , well[crv_1], where = well[crv_1] >= item,  color = color)
############################################################  

fig.subplots_adjust(wspace = 0.1)
plt.tight_layout()
plt.show()

En resumen, hemos cargado los datos de un archivo .las y hemos hecho un gráfico utilizando la función fill_betweenx() de matplotlib para rellenar con color el área entre 2 curvas.

Seguidamente, graficamos otro log utilizando fill_betweenx() y un mapa de colores o cmap conjuntamente, para rellenar el área entre una curva y un valor del track con una paleta de colores.

Para finalizar, juntamos ambos gráficos en uno solo, lo cual, en teoría, nos va a ayudar con la interpretación petrofísica inicial que hacemos al mirar un log.

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