繪製地震擺動軌跡使用matplotlib

Question

我試圖重新使用matplotlib繪製的上述風格。

原始數據存儲在2D numpy數組中，其中快軸是時間。

繪製線條很簡單。我試圖有效地獲得陰影區域。

我當前的嘗試看起來像：

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 #lazy prevents polygon wrapping 
x = panel.ravel() 
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 scalars used to shift "traces" to plotting corrdinates 
trace_centers = np.linspace(1,284, panel.shape[-2]).reshape(-1,1) 
gain = 0.5 #scale traces 

#shift traces to plotting coordinates 
x = ((panel*gain)+trace_centers).ravel() 

#split coordinate vectors at each zero crossing 
xpoly = np.split(x, zero_crossings) 
ypoly = np.split(y, zero_crossings) 

#we only want the polygons which outline positive values 
if x[0] > 0: 
    steps = range(0, len(xpoly),2) 
else: 
    steps = range(1, len(xpoly),2) 

#turn vectors of polygon coordinates into lists of coordinate pairs 
polygons = [zip(xpoly[i], ypoly[i]) for i in steps if len(xpoly[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() 

，其結果是

這裏有兩個問題：

1. 它不填寫完全，因爲我對分裂數組索引最接近零交叉點，而不是確切的零交叉點。我假設計算每個過零點將是一個大計算命中。

2. 表現。考慮到問題的嚴重程度，我的筆記本電腦一秒鐘左右就能呈現出來，但這並不是那麼糟糕，但我想把它降低到100ms-200ms。

由於使用情況，我僅限於Python與numpy/scipy/matplotlib。有什麼建議麼？

跟帖：

原來線性內插過零點可以用很少的計算量來完成。通過在數據中插入插值，將負值設置爲nans，並使用一次調用pyplot.fill，可以在300ms左右繪製500,000個奇數採樣。

作爲參考，湯姆的方法在相同的數據下面花費了大約8秒。

下面的代碼假定輸入一個numpy recarray，其中dtype模仿一個地震unix頭/軌跡定義。

def wiggle(frame, scale=1.0): 
     fig = pylab.figure() 
     ax = fig.add_subplot(111)   
     ns = frame['ns'][0] 
     nt = frame.size 
     scalar = scale*frame.size/(frame.size*0.2) #scales the trace amplitudes relative to the number of traces 
     frame['trace'][:,-1] = np.nan #set the very last value to nan. this is a lazy way to prevent wrapping 
     vals = frame['trace'].ravel() #flat view of the 2d array. 
     vect = np.arange(vals.size).astype(np.float) #flat index array, for correctly locating zero crossings in the flat view 
     crossing = np.where(np.diff(np.signbit(vals)))[0] #index before zero crossing 
     #use linear interpolation to find the zero crossing, i.e. y = mx + c. 
     x1= vals[crossing] 
     x2 = vals[crossing+1] 
     y1 = vect[crossing] 
     y2 = vect[crossing+1] 
     m = (y2 - y1)/(x2-x1) 
     c = y1 - m*x1  
     #tack these values onto the end of the existing data 
     x = np.hstack([vals, np.zeros_like(c)]) 
     y = np.hstack([vect, c]) 
     #resort the data 
     order = np.argsort(y) 
     #shift from amplitudes to plotting coordinates 
     x_shift, y = y[order].__divmod__(ns) 
     ax.plot(x[order] *scalar + x_shift + 1, y, 'k') 
     x[x<0] = np.nan 
     x = x[order] *scalar + x_shift + 1 
     ax.fill(x,y, 'k', aa=True) 
     ax.set_xlim([0,nt]) 
     ax.set_ylim([ns,0]) 
     pylab.tight_layout() 
     pylab.show() 

完整的代碼在https://github.com/stuliveshere/PySeis

Answer 1

發表你可以很容易地fill_betweenx做到這一點。從文檔：

在兩條水平曲線之間製作填充多邊形。

呼叫簽名：

fill_betweenx（Y，X1，X2 = 0，其中=無，** kwargs）創建一個 PolyCollection填充x1和x2之間的區域，其中在那裏== TRUE

這裏的重要部分是where的論點。

所以，你想擁有x2 = offset，然後有where = x>offset

例如：

import numpy as np 
import matplotlib.pyplot as plt 

fig,ax = plt.subplots() 

# Some example data 
y = np.linspace(700.,900.,401) 
offset = 94. 
x = offset+10*(np.sin(y/2.)* 
     1/(10. * np.sqrt(2 * np.pi)) * 
     np.exp(- (y - 800)**2/(2 * 10.**2)) 
     ) # This function just gives a wave that looks something like a seismic arrival 

ax.plot(x,y,'k-') 
ax.fill_betweenx(y,offset,x,where=(x>offset),color='k') 

ax.set_xlim(93,95) 

plt.show() 

你需要做的fill_betweenx爲每個偏移。例如：

import numpy as np 
import matplotlib.pyplot as plt 

fig,ax = plt.subplots() 

# Some example data 
y = np.linspace(700.,900.,401) 
offsets = [94., 95., 96., 97.] 
times = [800., 790., 780., 770.] 

for offset, time in zip(offsets,times): 
    x = offset+10*(np.sin(y/2.)* 
     1/(10. * np.sqrt(2 * np.pi)) * 
     np.exp(- (y - time)**2/(2 * 10.**2)) 
     ) 

    ax.plot(x,y,'k-') 
    ax.fill_betweenx(y,offset,x,where=(x>offset),color='k') 

ax.set_xlim(93,98) 

plt.show() 

Answer 2

這是相當容易的，如果你有SEGY格式和/或txt格式的地震道做（你需要讓他們在.txt格式最終）。花了很長時間才找到最好的方法。對python和編程也是新的，所以請溫和。

將SEGY文件轉換爲.txt文件我使用了SeiSee（http://dmng.ru/en/freeware.html;不介意俄羅斯網站，它是一個合法的程序）。爲了加載和顯示你需要numpy和matplotlib。

以下代碼將加載地震痕跡，轉置它們並繪製它們。很明顯，你需要加載你自己的文件，改變垂直和水平範圍，並用vmin和vmax來玩一下。它也使用灰色的顏色表。該代碼會產生這樣的圖像：http://goo.gl/0meLyz

import numpy as np 
import matplotlib.pyplot as plt 

traces = np.loadtxt('yourtracestxtfile.txt') 
traces = np.transpose(traces) 

seismicplot = plt.imshow(traces[3500:4500,500:900], cmap = 'Greys',vmin = 0,vmax = 1,aspect = 'auto') #Tip: traces[vertical range,horizontal range]