Petrophysical Composite for Interpretation Work

Let's build a petrophysical composite display for a single well with:

  • Logs (GammaRay Track, Resistivity Track, Neutron-Density Track, Sonic, Synthetic)
  • Stiplog Track
  • Textual cutting descriptions
  • Markers or tops annotations
  • Can switch between time and depth

Let's use the P-129 well from the demo data


In [ ]:
%cd ~/Dropbox/geotransect_test/scripts

In [ ]:
import matplotlib.pyplot as plt
import numpy as np
% matplotlib

In [ ]:
from las import LASReader

In [ ]:
from lithology import intervals_from_las3_string

In [ ]:
from utils import summarize_rock
from legends import helen, legend

In [ ]:
r = {'code': 4,
     'colour':'Red',
     'grainsize': None,
     'lithology': 'Sandstone',
     'map colour': 'F2FF42' ,
     'summary': 'Sandstone, F-M, Grey' ,
     'width': 4}

In [ ]:
# well_dir = 'data/wells/P-129/'

In [ ]:
cd ../data/wells/

In [ ]:
lithwell = LASReader('P-129/lithology_log/P-129_striplog.las', unknown_as_other=True)

In [ ]:
striplog = intervals_from_las3_string(lithwell.other)

In [ ]:
striplog.keys()

In [ ]:


In [ ]:
logs = LASReader('P-129/wireline_log/P-129_out.LAS')

In [ ]:
Z = logs.data['DEPT']
GR = logs.data['GR']
DT = logs.data['DT'] 
DPHISS = logs.data['DPHI_SAN']
NPHISS = logs.data['NPHI_SAN']
DTS = logs.data['DTS']
RT_HRLT = logs.data['RT_HRLT']
RHOB = logs.data['RHOB']

In [ ]:
log_dict = {
            'GR':GR, 
            'DT':DT, 
            'DPHI_SAN':DPHISS,
            'NPHI_SAN':NPHISS,
            'DTS':DTS, 
            'RT_HRLT':RT_HRLT,
            'RHOB':RHOB
            }

we only want to plot a segment of the well between a top and base


In [ ]:
from matplotlib.patches import Rectangle
from matplotlib.ticker import Formatter, FixedLocator

In [ ]:
import csv
from StringIO import StringIO
import re

In [ ]:
top = 500
bottom = 1000
zrange = (top, bottom)
strip_width = 100

In [ ]:
def plot_striplog(ax, striplog, legend, width=1, ladder=False, summaries=False, minthick=1, alpha=0.75):
    
    for t, b, l in zip(striplog['tops'],striplog['bases'], striplog['liths']):
        origin = (0, t)
        colour = '#' + l['map colour'].strip()
        thick = b - t
        
        if ladder:
            w = legend[colour[1:]]['width']
        else:
            w = width
            
        rect = Rectangle(origin, w, thick, color=colour, edgecolor='k', 
                         linewidth=1.0, alpha=alpha)
        ax.add_patch(rect)

        if summaries:
            if thick >= minthick:
                ax.text(6, t+thick/2, summarize(l), verticalalignment='center')

    return ax

In [ ]:
def get_curve_params(abbrev, fname):
    """
    returns a dictionary of petrophysical parameters for 
    plotting purposes
    """
    params = {}
    with open(fname, 'rU') as csvfile:
        reader = csv.DictReader(csvfile) 
        for row in reader:
            if row['acronymn'] == abbrev:
                params['acronymn'] = row['acronymn']
                params['track'] = int(row['track'])
                params['units'] = row['units']
                params['xleft'] = float(row['xleft'])
                params['xright'] = float(row['xright'])
                params['logarithmic'] = row['logarithmic']
                params['hexcolor']= row['hexcolor']
                params['fill_left_cond']= bool(row['fill_left_cond'])
                params['fill_left']= row['fill_left']
                params['fill_right_cond']= bool(row['fill_right_cond'])
                params['fill_right']= row['fill_right']
                params['xticks'] = row['xticks'].split(',')
    return params

In [ ]:
pwd

In [ ]:
fname = '../../scripts/Petrophysics_display_template.csv'
params = get_curve_params('RT_HRLT', fname)

In [ ]:
import re

In [ ]:
test = ['0.2', '2.0', '20.0', '200']
string = [re.sub("'","'", s) for s in test]
print string

In [ ]:
print " * curves we have are: *"
for curve in log_dict:
    print "       ", curve

In [ ]:
# Subplots adjustment params

left  = 0.125  # the left side of the subplots of the figure
right = 0.9    # the right side of the subplots of the figure
bottom = 0.1   # the bottom of the subplots of the figure
top = 0.9      # the top of the subplots of the figure
wspace = 0.1   # the amount of width reserved for blank space between subplots
hspace = 0.5   # the amount of height reserved for white space between subplots

Fixing editing / despiking the logs


In [ ]:
def rolling_window(a, window):
        shape = a.shape[:-1] + (a.shape[-1] - window + 1, window)
        strides = a.strides + (a.strides[-1],)
        rolled = np.lib.stride_tricks.as_strided(a, 
                                                 shape=shape, 
                                                 strides=strides)
        rolled = np.median(rolled, -1)
        rolled = np.pad(rolled, window/2, mode='edge')
        return rolled
    
def despike(curve, curve_sm, max_clip): 
    spikes = np.where(curve - curve_sm > max_clip)[0]
    spukes = np.where(curve_sm - curve > max_clip)[0]
    out = np.copy(curve)
    out[spikes] = curve_sm[spikes] + max_clip  # Clip at the max allowed diff
    out[spukes] = curve_sm[spukes] - max_clip  # Clip at the min allowed diff
    return out

In [ ]:
w = 31   # window length must be an even number
frac = 0.05 # fraction of extreme values to clip (small == big clip)

cmin = 1.0
cmax = 3.5

smRHOB = rolling_window(RHOB,w)#, cmin, cmax
dsRHOB = despike(RHOB, smRHOB, frac*np.mean(RHOB))

values = smRHOB

plt.figure(figsize = (3,10))
#plt.plot(np.clip(RHOB, cmin, cmax), Z, c = 'b', alpha = 0.25)
#plt.plot(dsRHOB, Z, c = 'r', alpha = 0.25)
plt.plot(values, Z, c = 'g', alpha = 0.55)

plt.xlim([-1.0, 3.5])
plt.ylim(Z[-1], Z[0])

In [ ]:
import matplotlib.ticker as ticker
import matplotlib as mpl
import matplotlib.transforms as transforms

In [ ]:
% matplotlib

In [ ]:
# make plot
window = 51 # window length for smoothing must be an odd integer
frac = 0.05
p = 4
lw = 1
smooth = True
has_striplog = True
height = 3*p  # in inches
width = 2*p # in inches
fs = 12  #font size for curve labels

ntracks = 4
naxes = 0
ncurv_per_track = np.zeros(4)

if has_striplog:
    ncurv_per_track[0]=1

for curve, values in log_dict.iteritems(): 
    naxes += 1
    params =  get_curve_params(curve, fname)
    ncurv_per_track[params['track']] += 1

fig, axs = plt.subplots(1, p, figsize=(width, height))
fig.subplots_adjust(left, bottom, right, top, wspace, hspace)
fig.set_facecolor('w')

axss = axs[0]
axs0 = [axss, axs[0].twiny(), axs[0].twiny()]
axs1 = [axs[1]]
axs2 = [axs[2]]
axs3 = [axs[3]]

axs = [axs0, axs1, axs2, axs3]

plot_striplog(axs0[0], striplog, legend, width = 5, alpha = 0.75, ladder=True)

axs0[0].set_ylim([striplog['bases'][-1], 0])

# Plot each curve with a white fill to fake the curve fill.

label_shift = np.zeros(len(axs))

for curve, values in log_dict.iteritems():
    
    params = get_curve_params(curve, fname)
    
    i = params['track']
    
    if i == 0:
        j = 1
    
    j = 0  # default number of tracks to index into
    
    ncurves = ncurv_per_track[i]
    
    label_shift[i] += 1
    
    units = '$%s$' % params['units']
        
    #crvtxt = params['acronymn']+ '\n' + '\n' + units

    linOrlog = params['logarithmic']
    
    sxticks = np.array(params['xticks'])
    xticks = np.array(sxticks, dtype = float)
    
    if linOrlog == 'log':
        midline = np.log(np.mean(xticks))
        xpos = midline 
    else:
        midline = np.mean(xticks)
        xpos = midline 
    
    if smooth == True:
        values = rolling_window(values,window)
        
    if params['acronymn']=='GR':
        j = 1  # second axis in first track
        label_shift[i] = 1
    if params['acronymn']=='RHOB':
        j = 2
        label_shift[i] = 1
        
    #fill_left
    if params['fill_left_cond']:
        print i,j
        axs[i][j].fill_betweenx(Z, params['xleft'], values,
                            facecolor='white', alpha = 1.0, zorder = 11)
    # fill right        
    if params['fill_right_cond']:
        axs[i][j].fill_betweenx(Z, values, params['xright'],
                            facecolor='white', alpha = 1.0, zorder = 11)
    # plot curve        
    axs[i][j].plot(values, Z, color = params['hexcolor'], lw = lw)
        
    # set scale of curve
    axs[i][j].set_xlim([params['xleft'],params['xright']])
    axs[i][j].xaxis.set_ticks([xticks])
    axs[i][j].set_xticklabels([sxticks])
    
    for label in axs[i][j].axes.xaxis.get_ticklabels():
        label.set_rotation(90)
            
    axs[i][j].tick_params(axis='x', direction='out')
        
    # set curve label and units beneath it
    # if this is the first curve on the axis
    
    # curve label
    
    trans = transforms.blended_transform_factory(axs[i][j].transData, axs[i][j].transData)
    axs[i][j].text(xpos, -140 - 40*(label_shift[i] - 1), params['acronymn'],  
                   horizontalalignment='center',                    
                   verticalalignment='bottom',
                   fontsize=12, color=params['hexcolor'],
                   transform=trans)
    # curve units
    if label_shift[i] <= 1:
        axs[i][j].text(xpos, -120, units,
                   horizontalalignment='center',
                   verticalalignment='top',
                   fontsize=12, color='k',
                   transform=trans)
        # set ticks and tick labels
        #axs[i][j].set_xticks(xticks)
        #axs[i][j].set_xticklabels(xticks)
    
    axs[i][j].set_xscale(linOrlog)
    axs[i][j].set_ylim([striplog['bases'][-1], 0]) 
    axs[i][j].xaxis.tick_top()
    axs[i][j].xaxis.set_label_position('top')
    
    if i != 0:
        axs[i][j].set_yticks([])
    print '...............'
fig.show()

In [ ]:
% matplotlib
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot([1,3,4,5])
alpha = ['a','b','c','d','e','f','g']
ax.xaxis.set_ticklabels([])
ax.xaxis.set_ticks([])
ax.yaxis.set_ticks([1.1,1.2,2.1,4.5])

In [ ]:
plt.show()

In [ ]:


In [ ]:


In [ ]:


In [ ]: