In [1]:
from toolbox.processing import *
#%ls /home/stewart/su/2d_land_data/2D_Land_data_2ms/
file = "/home/stewart/su/2d_land_data/2D_Land_data_2ms/su/Line_001.su"
#file = "/home/sfletcher/Downloads/2d_land_data/2D_Land_data_2ms/Line_001.su"
#initialise file
data, params = initialise(file)
In [ ]:
#no coordinates in the headers, but we know energy point number and channel.
#%cat /home/stewart/su/2d_land_data/2D_Land_data_2ms/Line_001.TXT
#%cat /home/stewart/su/2d_land_data/2D_Land_data_2ms/Line_001.SPS
#%cat /home/stewart/su/2d_land_data/2D_Land_data_2ms/Line_001.RPS
#%cat /home/stewart/su/2d_land_data/2D_Land_data_2ms/header
In [ ]:
import numpy as np
from matplotlib import collections
dmap = np.memmap(file, dtype=toolbox.typeSU(1501), mode='r')
eps = np.unique(dmap['ep'])
for ep in eps[:1]:
panel = dmap[dmap['ep'] == ep].copy()
panel = toolbox.agc(panel, None, **params)
trace_centers = np.linspace(1,284, panel.size).reshape(-1,1)
trace_width = 284/(panel.size*0.5)
x = panel['trace'].copy()
x += trace_centers
y = np.meshgrid(np.arange(1501), np.arange(284))[0]
x = np.split(x.ravel(), 284)
y = np.split(y.ravel(), 284)
bits = [zip(x[a],y[a]) for a in range(len(x))]
fig = pylab.figure()
ax = fig.add_subplot(111)
col1 = collections.LineCollection(bits)
col1.set_color('k')
ax.add_collection(col1, autolim=True)
ax.autoscale_view()
pylab.xlim([0,284])
pylab.ylim([0,1500])
ax.set_ylim(ax.get_ylim()[::-1])
pylab.tight_layout()
pylab.show()
In [ ]:
import numpy as np
from matplotlib import collections
dmap = np.memmap(file, dtype=toolbox.typeSU(1501), mode='r')
eps = np.unique(dmap['ep'])
for ep in eps[:1]:
panel = dmap[dmap['ep'] == ep].copy()
panel = toolbox.agc(panel, None, **params)
trace_centers = np.linspace(1,284, panel.size).reshape(-1,1)
scalar = 284/(panel.size*0.5)
panel['trace'][:,-1] = np.nan
x = panel['trace'].ravel()
x[x < 0] = 0
y = np.meshgrid(np.arange(1501), np.arange(284))[0].ravel()
zero_crossings = np.where(x == 0)[0]+1
zero_crossings = zero_crossings[np.diff(zero_crossings) == 1]
#zero_crossings = np.where(np.diff(np.signbit(x)))[0]+1
x = ((panel['trace']*scalar)+trace_centers).ravel()
xverts = np.split(x, zero_crossings)
yverts = np.split(y, zero_crossings)
polygons = [zip(xverts[i], yverts[i]) for i in range(0, len(xverts)) if len(xverts[i]) > 2]
xlines = np.split(x, 284)
ylines = np.split(y, 284)
lines = [zip(xlines[a],ylines[a]) for a in range(len(xlines))]
fig = pylab.figure()
ax = fig.add_subplot(111)
col = collections.PolyCollection(polygons)
col.set_color('k')
ax.add_collection(col, autolim=True)
col1 = collections.LineCollection(lines)
col1.set_color('k')
ax.add_collection(col1, autolim=True)
ax.autoscale_view()
pylab.xlim([0,284])
pylab.ylim([0,1500])
ax.set_ylim(ax.get_ylim()[::-1])
pylab.tight_layout()
pylab.show()
In [ ]:
import numpy as np
from matplotlib import collections
%pylab tk
def polytrace(data, **kwargs):
segs = []
segl = []
nt = data.shape[-2]
for i in range(nt):
trace = data[i]
line = list(zip(trace, np.arange(1501)))
segl.append(line)
xx = trace #np.ma.array(trace, mask=(trace <= 0))
yy = np.arange(1501) #np.ma.array(np.arange(0,1501), mask=(trace <= 0))
curve = [(0, 0)]
curve.extend(list(zip(xx, yy)))
curve.extend([(0, 1501)])
segs.append(curve)
#print segs[0]
print ''
return segs, segl
#first lets do some checks. does of energy points should equal number of records?
print np.unique(data['ep']).size, np.unique(data['fldr']).size
#no duplicates - that makes it easier.
print 251*284
#284 traces per shot, 2 aux traces . lets have a look
dmap = np.memmap(file, dtype=toolbox.typeSU(1501), mode='r')
eps = np.unique(dmap['ep'])
for ep in eps[:1]:
panel = dmap[dmap['ep'] == ep].copy()
panel = toolbox.agc(panel, None, **params)
trace_centers = np.linspace(1,284, panel.size).reshape(-1,1)
trace_width = 284/(panel.size*0.5)
buf = panel['trace'].copy()
buf *= trace_width
segs, segl = polytrace(buf)
fig = pylab.figure()
ax = fig.add_subplot(111)
offs = (10.0, 0.0)
offs = list(zip(np.arange(238), np.zeros(238)))
col = collections.PolyCollection(segs, offsets=offs)
col.set_color('k')
ax.add_collection(col, autolim=True)
#col1 = collections.LineCollection(segl, offsets=offs)
#col1.set_color('k')
#ax.add_collection(col1, autolim=True)
ax.autoscale_view()
pylab.xlim([0,284])
pylab.ylim([0,1500])
ax.set_ylim(ax.get_ylim()[::-1])
pylab.tight_layout()
pylab.show()
In [ ]:
import numpy as np
from matplotlib import collections
import matplotlib.pyplot as pylab
#make some oscillating data
panel = np.meshgrid(np.arange(1501), np.arange(284))[0]
panel = np.sin(panel)
#generate coordinate vectors.
panel[:,-1] = np.nan #prevent wrapping when flatten 2d array
x = panel.flatten()
y = np.meshgrid(np.arange(1501), np.arange(284))[0].ravel()
#find indexes of each zero crossing
zero_crossings = np.where(np.diff(np.signbit(x)))[0]+1
#calculate scalar used to shift "traces" to plot corrdinates
trace_centers = np.linspace(1,284, panel.shape[-2]).reshape(-1,1)
gain = 0.5 #scale traces
#shift traces to plotting coordinate
x = ((panel*gain)+trace_centers).ravel()
#split each vector at each zero crossing
xverts = np.split(x, zero_crossings)
yverts = np.split(y, zero_crossings)
#we only want the vertices which outline positive values
if x[0] > 0:
steps = range(0, len(xverts),2)
else:
steps = range(1, len(xverts),2)
#turn vectors of coordinates into lists of coordinate pairs
polygons = [zip(xverts[i], yverts[i]) for i in steps if len(xverts[i]) > 2]
#this is so we can plot the lines as well
xlines = np.split(x, 284)
ylines = np.split(y, 284)
lines = [zip(xlines[a],ylines[a]) for a in range(len(xlines))]
#and plot
fig = pylab.figure()
ax = fig.add_subplot(111)
col = collections.PolyCollection(polygons)
col.set_color('k')
ax.add_collection(col, autolim=True)
col1 = collections.LineCollection(lines)
col1.set_color('k')
ax.add_collection(col1, autolim=True)
ax.autoscale_view()
pylab.xlim([0,284])
pylab.ylim([0,1500])
ax.set_ylim(ax.get_ylim()[::-1])
pylab.tight_layout()
pylab.show()
In [16]:
import numpy as np
from matplotlib import collections
dmap = np.memmap(file, dtype=toolbox.typeSU(1501), mode='r')
eps = np.unique(dmap['ep'])
for ep in eps[:1]:
panel = dmap[dmap['ep'] == ep].copy()
panel = toolbox.agc(panel, None, **params)
panel['trace'][:,-1] = np.nan
trace_centers = np.linspace(1,284, panel.size).reshape(-1,1)
scalar = 284/(panel.size*0.5)
y = np.meshgrid(np.arange(1501), np.arange(284))[0].ravel()
offsets = (np.meshgrid(np.arange(1501), np.arange(284))[1]+1).ravel()
x = ((panel['trace']*scalar)+trace_centers).ravel()
fig,ax = plt.subplots()
#or i in range(284):
#ax.plot(x[i],y[i],'k-')
ax.fill_betweenx(y,offsets,x,where=(x>offsets),color='k')
pylab.xlim([0,284])
pylab.ylim([0,1500])
ax.set_ylim(ax.get_ylim()[::-1])
pylab.tight_layout()
pylab.show()