In [1]:
import numpy as np
import matplotlib.pyplot as plt
% matplotlib inline
In [2]:
# Make dummy seismic data
ntr = 300
nsamps = 300
dz = 15.0 # sample interval in metres
seismic = np.random.rand(nsamps, ntr)
srd = 0.0 #seismic reference datum
max_z = 4500.0 # max depth of cross-section
z_dom = np.arange(0, max_z, dz) # vertical axes domain extents
start_trace = 50 #left most position of first trace along x-section
In [3]:
# Make cross-section plane
dx = 100 # offset sample interval
xmin = 0 # xmin
xmax = 600 * dx # xmax
zmin = -500 # zmin
zmax = 5000 # zmax
nx = xmax - xmin
nz = zmax - zmin
ve = 5.0 # vertical exaggeration
xsec = np.zeros((nz, nx))
In [4]:
# Make dummy elevation profile
from matplotlib.collections import LineCollection
offset = np.arange(xmin, xmax+dx, dx)
c = 0.0001
k = -0.9
maxelev = 700.0
elevs = maxelev*np.exp(k*offset/30000)*(np.sin(c*offset))**2
# make_colors
colors = np.zeros(len(offset))
for i,val in enumerate(offset):
if i < 100:
colors[i] = 40.0
if 100 < i < 210:
colors[i] = 20.0
if 210 < i < 237:
colors[i] = 60.0
if 237 < i < 415:
colors[i] = 80.0
if 415 < i < 550:
colors[i] = 100.0
#print color
# 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:
scatter(offset, elevs, s= 50, c=colors, edgecolor = 'none', cmap='jet')
np.ptp(colors)
In [ ]:
# Make horizon(s)
seisx = np.arange(0, (ntr+1)*dx, dx)
hrz1 = 2750 - 0.15*seisx
hrz2 = 2500 - 0.17*seisx
plot(seisx + start_trace*dx, hrz1, color = 'orange')
plot(seisx + start_trace*dx, hrz2, color='y')
In [ ]:
# Make dummy strip log(s)
xpos = 17500 # well penetrates the xsex at this x pos.
kb = -400.0 # KB elevation
td = 1800.0 # Total Depth
width = 2000 #width of log in plot units (probably metres) for plotting
intervals = np.arange(kb, td, 100)
tops = intervals[:-1]
bases = intervals[1:]
liths = []
for i in range(len(intervals)-1):
liths.append(int(width*np.random.rand(1)))
In [ ]:
from matplotlib.patches import Rectangle
lith_fig = plt.figure(figsize=(3,10))
well_ax = lith_fig.add_subplot(111)
for t, b, w in zip(tops, bases, liths):
#print w
color = 'yellow'
if w > width/2:
color = 'brown'
rect1 = Rectangle((0,t), w, b - t, color=color, alpha = 0.5)
well_ax.add_patch(rect1)
well_ax.set_xlim([0,1.25*width])
well_ax.set_ylim([np.amax(bases), 0])
well_ax.set_xticks([])
well_ax.set_yticklabels([])
well_ax.set_yticks([])
# Hide the right and top spines
well_ax.spines['right'].set_visible(False)
well_ax.spines['top'].set_visible(False)
well_ax.spines['bottom'].set_visible(False)
# Turn on the right spine
well_ax.spines['left'].set_position(('data', 0))
#well_ax.patch.set_facecolor('white')
well_ax.patch.set_alpha(0.5)
In [ ]:
In [5]:
fsize = (17,22) # figure size in inches
fbkgd = '#F7F7F7' # light grey so we can see the "canvas"
xsec_extents = [0.1,0.1,0.8,0.8]
In [5]:
In [6]:
# make figure instance
fig = plt.figure(figsize = fsize, facecolor = fbkgd) # create figure, set size, set facecolor
# make cross-section plane
xsec_ax = fig.add_axes(xsec_extents)
xsec_ax.imshow(xsec, extent = [xmin, xmax, zmax, zmin], aspect = ve, alpha = 0.00)
fig.canvas.draw()
xlabels=np.array(xsec_ax.get_xticks()/1000.0)
ylabels=np.array(xsec_ax.get_yticks()/1000.0)
xsec_ax.set_xticklabels(xlabels.astype(int))
#xsec_ax.set_yticklabels(ylabels.astype(int))
#xsec_ax.set_ylabel('Depth [km]',labelpad=20)
new_xlabels = [r'%s m' % l.get_text() for l in
xsec_ax.xaxis.get_ticklabels()]
new_ylabels = [r'%s m' % l.get_text() for l in
xsec_ax.yaxis.get_ticklabels()]
xsec_ax.xaxis.set_ticklabels(new_xlabels)
xsec_ax.yaxis.set_ticklabels(new_ylabels)
xsec_ax.set_xlabel('Distance')
xsec_ax.text(xmin, zmax,' transect name ',horizontalalignment='left',
verticalalignment='bottom')
#xsec_ax.grid()
# make elevation plot:
elev_ax = fig.add_axes(xsec_extents)
# elev_ax.plot(offset, -elevs, 'k', lw = 2)
elev_ax.scatter(offset, -elevs, s= 50, c=colors, edgecolor = 'none', cmap='jet', zorder = 1)
# make seismic axes (or more than one)
seis_ax = fig.add_axes(xsec_extents)
seis_ax.imshow(seismic, extent = [start_trace*dx, start_trace*dx + ntr*dx, max_z, srd], cmap = 'Greys', aspect = ve, alpha = 0.2)
seis_ax.plot(seisx + start_trace*dx, hrz1, lw = 3, c = 'orange')
seis_ax.plot(seisx + start_trace*dx, hrz2, lw = 3, c = 'y')
seis_ax.set_xlim(left = xmin, right = xmax)
seis_ax.set_ylim(bottom = zmax, top = zmin)
seis_ax.text(start_trace*dx + ntr*dx, max_z,' seismic line name ',horizontalalignment='right',
verticalalignment='bottom')
# make strip log of well (try to do more than one)
well_ax = fig.add_axes(xsec_extents)
well_ax.text(xpos, td*1.05,'Well\nname',horizontalalignment='center',
verticalalignment='top')
for t, b, w in zip(tops, bases, liths):
color = 'yellow'
if w > width/2:
color = 'brown'
rect1 = Rectangle((xpos, t), w, b - t, color=color, alpha = 0.5)
well_ax.add_patch(rect1)
well_ax.vlines(x=xpos, ymin=kb, ymax=td, color='k', zorder=2)
tops = {'Alpha': td*0.2, 'Beta' : td*0.5, 'Epsilon' : td*0.75}
kw = dict(rotation=0, verticalalignment='center', horizontalalignment='right')
tops_bbox=dict(facecolor='white', edgecolor='black',
boxstyle='round', alpha = 0.5)
for name, val in tops.iteritems():
well_ax.text(xpos - 500.0, val, name, fontsize = 8,
bbox = tops_bbox, **kw)
well_ax.scatter(xpos, val, s=50, c='k')
# insert map area
"""
from mpl_toolkits.axes_grid1.inset_locator import inset_axes, zoomed_inset_axes
map_ax = fig.add_axes([0.65, 0.375, 0.25, 0.10])
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
"""
In [ ]:
In [ ]:
In [ ]: