In [1]:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
from itertools import product
% matplotlib inline
Get data
In [4]:
cd ../..
In [5]:
cd Onshore_Data_Cleanup/Data/seismic/SEGY/Windsor/2007-Elmworth\ Windsor
In [1]:
cd SEGY/PSTM/
In [47]:
from obspy.segy.segy import readSEGY
datadir = '.'
#segyfile = 'ELM07-KENN-01_PSTM.SGY'
segyfile = 'NS00_Wind-06_PSTM.sgy'
filename = datadir + "/" + segyfile
print "SEGY File: " + filename
# Read all traces
section = readSEGY(filename,unpack_headers=True)
print(section)
elevs = []
for i,trace in enumerate(section.traces):
elevs.append(trace.header.receiver_group_elevation)
nsamples = section.traces[-1].header.number_of_samples_in_this_trace
dt = (1/1000.0) * section.traces[-1].header.sample_interval_in_ms_for_this_trace
ntraces = len(section.traces)
tbase = np.arange(0,nsamples * dt, dt)
tstart = 0
tend = np.amax(tbase)
aspect = float(ntraces) / float(nsamples)
print ('ntraces', ntraces)
print ('nsamples', nsamples)
print ('dt', dt)
data = np.zeros((nsamples, ntraces))
for i, trace in enumerate(section.traces):
data[:,i] = trace.data
line_extents = {'first_trace': 1,
'last_trace': ntraces,
'start_time': tstart,
'end_time':tend
}
clip_val = np.percentile(data,99.0)
print "clip_val", clip_val
print "max_val", np.amax(data)
print "min_val", np.amin(data)
print "tstart", tstart
print "tend", tend
largest = max(np.amax(data), abs(np.amin(data)))
For our page, let's start with the width and heigh dimensions equal to the Golden ratio; 1.618
In [48]:
In [8]:
plt.imshow(data, cmap='Greys')
Out[8]:
In [9]:
def figsize(scale):
fig_width_pt = 469.755 # Get this from LaTeX using \the\textwidth
inches_per_pt = 1.0/72.27 # Convert pt to inch
golden_mean = (np.sqrt(5.0)-1.0)/2.0 # Aesthetic ratio (you could change this)
fig_width = fig_width_pt*inches_per_pt*scale # width in inches
fig_height = fig_width*golden_mean # height in inches
fig_size = (fig_width,fig_height)
return fig_size
In [10]:
sc = 2
page_dims = figsize(scale = sc)
In [11]:
page_dims
Out[11]:
In [12]:
# Allocate areas of the plot in square inches.
nrows = int(np.floor(page_dims[0])) # number of complete inches in x
ncols = int(np.floor(page_dims[-1])) # number of complete inches in y
print nrows, ncols
In [13]:
def make_grid_layout(nrows, ncols):
gs = gridspec.GridSpec(ncols, nrows)
gsb = gridspec.GridSpecBase(ncols, nrows)
fig = plt.figure(figsize = page_dims,
facecolor = 'w',
edgecolor='k',
dpi=None,
frameon=True)
x = np.arange(0, ncols, 1)
y = np.arange(0, nrows, 1)
xy = np.meshgrid(x,y,sparse=True)
index = [(a, b) for a in x for b in y]
for idx, g in zip(index,gs):
ax = fig.add_subplot(g)
ax.set_xticks([])
ax.set_yticks([])
ax.text(0.5,0.5, str(idx[0])+', '+str(idx[1]),ha='center',va='center',size=sc*5,alpha=.5)
fig.tight_layout()
fig.show()
return fig
tiles = make_grid_layout(nrows, ncols)
In [14]:
def make_grid_layout(nrows, ncols):
gs = gridspec.GridSpec(ncols, nrows)
gsb = gridspec.GridSpecBase(ncols, nrows)
fig = plt.figure(figsize = page_dims,
facecolor = 'w',
edgecolor='k',
dpi=None,
frameon=True)
x = np.arange(0, ncols, 1)
y = np.arange(0, nrows, 1)
xy = np.meshgrid(x,y,sparse=True)
index = [(a, b) for a in x for b in y]
# Make Background Grid
for idx, g in zip(index,gs):
ax = fig.add_subplot(g)
ax.set_xticks([])
ax.set_yticks([])
ax.patch.set_facecolor('grey')
ax.patch.set_alpha(0.25)
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.spines['left'].set_visible(False)
ax.spines['bottom'].set_visible(False)
ax.text(0.5,0.5, str(idx[0])+', '+str(idx[1]),ha='center',va='center',size=sc*5,alpha=.5)
# Make Overlay Grid
#'name', (start from top, span_y, start_from_left, end)
my_axes = {
'main': (2,7,2,-4),
'left':(2,7,0,2),
'right':(3,8,-4,0),
'top line plot':(0,2,2,9),
'base map':(0,3,9,13),
'logo':(0,2,0,2),
'bottom line plot':(-1,0,2,-4),
'lower left info':(-1,0,0,2),
'lower right info':(7,8,9,13),
#'inset a':(2,7,4,5),
#'scale':(7,9,7,8),
#'inset elevation':(1,2,2,10)
}
for key,val in my_axes.iteritems():
my_ax = fig.add_subplot(gs[val[0]:val[1],val[2]:val[3]])
my_ax.set_xticks([])
my_ax.set_yticks([])
my_ax.patch.set_facecolor('red')
my_ax.patch.set_alpha(0.2)
my_ax.text(0.5,0.5, key, ha='center',va='center',size=sc*8,alpha=1.0)
fig.tight_layout()
fig.show()
return fig
region_tiles = make_grid_layout(nrows, ncols)
In [15]:
def make_regions_layout(fig, nrows, ncols, tiles=False):
gs = gridspec.GridSpec(ncols, nrows)
gsb = gridspec.GridSpecBase(ncols, nrows)
fig = plt.figure(figsize = page_dims,
facecolor = 'w',
edgecolor='k',
dpi=None,
frameon=True)
x = np.arange(0, ncols, 1)
y = np.arange(0, nrows, 1)
xy = np.meshgrid(x,y,sparse=True)
index = [(a, b) for a in x for b in y]
# Make Background Grid
if tiles:
for idx, g in zip(index,gs):
ax = fig.add_subplot(g)
ax.set_xticks([])
ax.set_yticks([])
ax.patch.set_facecolor('grey')
ax.patch.set_alpha(0.25)
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.spines['left'].set_visible(False)
ax.spines['bottom'].set_visible(False)
ax.text(0.5,0.5, str(idx[0])+', '+str(idx[1]),ha='center',va='center',size=sc*5,alpha=.5)
# Make Overlay Grid
#'name', (start from top, span_y, start_from_left, end)
my_axes = {
'main': (2,7,1,-5),
'left':(2,7,0,1),
'right log':(2,7,-5,-4),
'top line plot':(0,2,1,-5),
'upper right':(0,2,8,13),
'logo':(0,2,0,1),
'bottom line plot':(-1,0,1,-5),
'lower left info':(-1,0,0,1),
'base map':(5,8,9,13),
#'inset a':(2,7,4,5),
#'scale':(7,9,7,8),
#'inset elevation':(1,2,2,10)
}
#gs.update(wspace=0.025, hspace = 0.025) # set the spacing between axes
for key,val in my_axes.iteritems():
my_ax = fig.add_subplot(gs[val[0]:val[1],val[2]:val[3]])
my_ax.set_xticks([])
my_ax.set_yticks([])
my_ax.patch.set_facecolor('red')
my_ax.patch.set_alpha(0.2)
my_ax.text(0.5,0.5, key, ha='center',va='center',size=sc*8,alpha=1.0)
gs.update(wspace=0.025, hspace = 0.025) # set the spacing between axes
#fig.tight_layout()
#fig.show()
return fig, gs, my_axes
canvas = plt.figure(figsize = (13,8), facecolor='#E8E8E8')
regions, gs, my_axes = make_regions_layout(canvas,nrows,ncols,tiles=False)
regions.show()
In [16]:
from mpl_toolkits.basemap import Basemap
In [25]:
from mpl_toolkits.axes_grid import make_axes_locatable
from mpl_toolkits.axes_grid1.axes_divider import HBoxDivider, VBoxDivider
import mpl_toolkits.axes_grid1.axes_size as Size
from matplotlib.collections import LineCollection
from matplotlib.colors import ListedColormap, BoundaryNorm
In [133]:
def make_basemap(fig, gs, map_axes):
map_ax = fig.add_subplot(gs[map_axes[0]:map_axes[1],map_axes[2]:map_axes[3]])
bmap = Basemap(projection='merc', lat_0=45, lon_0=-63,
resolution = 'h', area_thresh = 0.1,
llcrnrlon=-67.0, llcrnrlat=43.5,
urcrnrlon=-59.0, urcrnrlat=47.5)
bmap.drawcoastlines(color='#acacaf')
bmap.drawcountries()
bmap.fillcontinents(color='#d8d8db')
bmap.drawmapboundary(color = 'white')
bmap.ax = map_ax
return fig
def side_plot(fig, props, size, pad, ytext):
vplot = fig.add_subplot(gs[my_axes['main'][0]:my_axes['main'][1],
my_axes['main'][3]:my_axes['main'][3]+1],
)
vplot.plot(data[:,i], tbase, 'k', lw = 0.5, alpha = 0.5)
vplot.set_ylim(tbase[-1],tbase[0])
vplot.set_xlim(-3*clip_val,3*clip_val)
ztimes = np.array([ 0, 250, 495, 760, 980, 1170, 1360, 1555, 1755, 1950, 2150, 2350])
z = [0, 500, 1000, 1500, 2000, 2500, 3000, 3500, 4000, 4500, 5000, 5500, 6000, 6500]
vplot.set_yticks(ztimes)
vplot.yaxis.tick_right()
vplot.set_ylim(ztimes[-1], ztimes[0]) # turn axis upside down
vplot.yaxis.set_label_position("right")
vplot.tick_params(axis='y', labelsize=6)
vplot.set_xticks([])
vplot.set_frame_on(True)
for spine in vplot.spines.values():
spine.set_edgecolor('#F3F3F3')
vplot.set_ylabel("well log", fontsize = 6)
vplot.set_ylabel("Depth [m]", fontsize = 8)
return vplot
# Data manipulation:
def make_segments(x, y):
'''
Create list of line segments from x and y coordinates, in the correct format for LineCollection:
an array of the form numlines x (points per line) x 2 (x and y) array
'''
points = np.array([x, y]).T.reshape(-1, 1, 2)
segments = np.concatenate([points[:-1], points[1:]], axis=1)
return segments
# Interface to LineCollection:
def colorline(x, y, z=None, cmap=plt.get_cmap('copper'), norm=True, linewidth=5, alpha=1.0):
'''
Plot a colored line with coordinates x and y
Optionally specify colors in the array z
Optionally specify a colormap, a norm function and a line width
'''
# Default colors equally spaced on [0,1]:
if z is None:
z = y
# Special case if a single number:
if not hasattr(z, "__iter__"): # to check for numerical input -- this is a hack
z = np.array([z])
z = np.asarray(z)
if norm:
norm = plt.Normalize(np.amin(y), np.amin(y))
segments = make_segments(x, y)
lc = LineCollection(segments, array=z, cmap=cmap, norm=norm, linewidth=linewidth, alpha=alpha)
#ax = plt.gca()
#ax.add_collection(lc)
return lc
def elevation_plot(canvas, elevs, colored_line, ax2share, ytext, props):
# horizontal plot
hplot = canvas.add_subplot(gs[my_axes['top line plot'][0]:my_axes['top line plot'][1],
my_axes['top line plot'][2]:my_axes['top line plot'][3]],
sharex=ax2share
)
hplot.add_collection(colored_line)
#hplot.colorline(np.arange(len(elevs)),elevs)
hplot.set_ylim(0,6*np.amax(elevs))
hplot.set_xlim(0,ntraces)
hplot.set_frame_on(True)
hplot.set_yticks([])
hplot.set_xticks([])
hplot.patch.set_alpha(0.1)
for spine in hplot.spines.values():
spine.set_edgecolor('#F3F3F3')
hplot.grid(True)
hplot.text(0.0, 0.5*(np.ptp(elevs)), ytext, props, rotation=0)
return hplot
def plot_xsec(canvas, data, my_axes):
xsec = canvas.add_subplot(gs[my_axes['main'][0]:my_axes['main'][1],
my_axes['main'][2]:my_axes['main'][3]])
xsec.imshow(data, cmap='Greys')
xsec.imshow(data, interpolation="nearest", cmap = cm.Greys)
xsec.set_ylim(nsamples*dt, 0)
xsec.set_ylabel("Two-way-time [ms]",fontsize = 8)
#xsec.set_yticks(np.arange(0,nsamples*dt, 500))
#xsec.set_yticks(np.arange(0,nsamples*dt, 500))
xsec.tick_params(axis='both', which='major', labelsize=8)
xsec.tick_params(axis='both', which='minor', labelsize=8) #extent = [nsamples*dt, 0, 0, ntraces])
xsec.set_yticks(np.arange(0,nsamples*dt, 500))
xsec.set_yticklabels(np.arange(0,nsamples*dt, 500))
return xsec
def plot_log(md, log, xloc):
logplot = canvas.add_subplot(gs[my_axes['main'][0]:my_axes['main'][1],
my_axes['main'][2]:my_axes['main'][3]])
logplot.plot(log + xloc, 0.5*md, 'g', lw = 0.5, alpha = 0.5)
logplot.set_ylim(tbase[-1],tbase[0]) #turn axis upside down
#vplot.tick_params(axis='y', labelsize=6)
logplot.set_xticks([])
logplot.set_yticks([])
logplot.set_frame_on(False)
return
In [134]:
cd /Users/Evan/Dropbox/Agile/NS_DOE/Onshore_Data_Cleanup/Data/wells/per_well/wells_by_basin/Windsor/P-130/wireline_data/project
In [135]:
ls
In [136]:
w = 'P-130_out.LAS'
In [137]:
well = np.loadtxt(w, skiprows=112)
In [138]:
Well = {'z':well[:,0], 'DT':well[:,3], 'RHOB':well[:-1],'GR':well[:,-5]}
md = Well['z'][50:-400]
GR = Well['GR'][50:-400]
In [139]:
#% matplotlib
canvas = plt.figure(figsize = (20,10), facecolor='#F2F2F2')
#regions, gs, my_axes = make_regions_layout(canvas,nrows,ncols,tiles=False)
make_basemap(canvas, gs, my_axes['base map'])
xsec = plot_xsec(canvas, data, my_axes)
plot_log(md, GR, xloc=300)
elev_text = 'Elevation [m]'
bbox = {'fc':'w', 'pad':0, 'ec' :'none', 'alpha' : 0.5}
props = {'ha':'left', 'va':'center', 'fontsize': 6, 'bbox':bbox}
#side_plot = side_plot(canvas, props, size='10%', pad = 0.0, ytext = elev_text)
elc = colorline(np.arange(ntraces), s_elevs)
hplot = elevation_plot(canvas, s_elevs, elc, xsec, elev_text, props)
In [ ]:
pwd
In [112]:
plot(Well['z'][50:-400],Well['GR'][50:-400])
Out[112]:
In [104]:
In [ ]: