【Python 绘制地震接收道集 Wigb函数】Plotting Seismic Wiggle Traces

使用Python Matplotlib绘制seismic wiggle traces

根据matlab版本的wigb.m改写。
Author(s): Xingong Li (Intseis, Integrated Seismic Solutions)
Copyright 1998-2003 Xingong
Revision: 1.2 Date: Dec/2002

代码

调用格式为wigb(data,scal,x,z,amx)
data为地震数据，数据结构为ndarray，shape：(nt,nc) (道长，道数)
scal 放缩系数
x,z 坐标轴
amx 最大值
后面四个参数是可选的
直接上代码：

import numpy as np
import matplotlib.pyplot as plt

def wigb(a=[],scal=0,x=0,z=0,amx=0):
#  WIGB: Plot seismic data using wiggles.
#
#  WIGB(a,scal,x,z,amx) 
#
#  IN  a:     地震数据 (a ndarray, traces are columns)
#      scale: multiple data by scal
#      x:     x轴(often offset)
#      z:     y轴 (often time)
#
#  Note
#
#    If only 'a' is enter, 'scal,x,z,amn,amx' are set automatically; 
#    otherwise, 'scal' is a scalar; 'x, z' are vectors for annotation in 
#    offset and time, amx are the amplitude range.
#
    if a==[]:
        nx, nz = 10, 10
        a = np.random.random((nz,nx))
        # print(a)
    nz, nx = a.shape

    trmx = np.max(np.abs(a),axis=0)
    if amx==0:
        amx = np.mean(trmx)
    if x==0:
        x = np.arange(1,nx+1,1)
        z = np.arange(1,nz+1,1)
    if scal==0:
        scal = 1
    
    if nx <=1:
        print('ERR:PlotWig: nx has to be more than 1')
        return 
    
# take the average as dx
    dx1 = np.abs(x[1:nx]-x[0:nx-1])
    dx = np.median(dx1)


    dz = z[1]-z[0]
    xmax, xmin = a.max(), a.min()

    a = a * dx / amx
    a = a * scal

    print(' PlotWig: data range [%f, %f], plotted max %f \n'%(xmin,xmax,amx))

    # set display range

    x1 = min(x) - 2.0*dx
    x2 = max(x) - 2.0*dx
    z1 = min(z) - dz
    z2 = max(z) - dz

    fig = plt.figure()
    plt.xlim([x1,x2])
    plt.ylim([z1,z2])
    plt.gca().invert_yaxis()


    zstart, zend = z[0], z[-1]

    fillcolor = [0, 0, 0]
    linecolor = [0, 0, 0]
    linewidth = 1.


    for i in range(0,nx):
        if not trmx[i]==0:
            tr = a[:,i]
            s = np.sign(tr)
            i1 = []
            for j in range(0,nz):
                if j==nz-1:
                    continue
                if not s[j]==s[j+1]:
                    i1.append(j)
            npos = len(i1)

            i1 = np.array(i1)
            if len(i1)==0:
                zadd=np.array(())
            else:
                zadd = i1 + tr[i1] / (tr[i1]-tr[i1+1])
            aadd = np.zeros(zadd.shape)

            zpos = np.where(tr>0)
            tmp = np.append(zpos,zadd)
            zz = np.sort(tmp)
            # iz =np.array(())
            # for j in range(0,len(tmp)):
            #     iz=np.append(iz,np.where(zz==tmp[j]))
            iz = np.argsort(tmp)
            aa = np.append(tr[zpos],aadd)
            iz=iz.astype(int)
            aa = aa[iz]


            if tr[0]>0:
                a0,z0=0,1.00
            else:
                a0,z0=0,zadd[0]
            if tr[nz-1]>0:
                a1,z1=0,nz
            else:
                a1,z1=0,zadd.max()
            
            zz = np.append(np.append(np.append(z0,zz),z1),z0)
            aa = np.append(np.append(np.append(a0,aa),a1),a0)

            zzz = zstart + zz * dz - dz
            plt.fill(aa+x[i], zzz+1, color=fillcolor)
            plt.plot(x[i]+[0,0],[zstart,zend],color=[1,1,1])

            plt.plot(x[i]+tr,z,color=linecolor,linewidth=linewidth)
        else:
            plt.plot([x[i],x[i]],[zstart,zend],color=linecolor,linewidth=linewidth)
    plt.show()


        

if __name__=='__main__':
    wigb()

样例