Let's build a petrophysical composite display for a single well with:
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
Functions taken from Well-tie Calculus notebook,
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 [ ]: