python成图wigb

 注意事项：1.需要安装obspy,numpy,matplotlib.

                   2.运行代码直接选择数据即可

                   3.可以读取obspy支持的所有主流地震数据

本文参考引用了matlab的wigb.m以及下面这个博主的代码(31条消息) 鱼鱼太好了吧的博客_CSDN博客-Madagascar,C语言,数据结构领域博主https://blog.csdn.net/qq_45317164?type=blog

import numpy as np
from tkinter import filedialog
from obspy import read
import matplotlib.pyplot as plt

# 读取数据存为list
st=read(filedialog.askopenfilename())
data_list = []
trace=len(st) # 道数
npts=len(st[0]) # 每道采样点数

for i in range(trace): # 读取stream为list数据
    data_list.append(st[i].data)
data_list=np.array(data_list)

print(f'数据共 {trace} 道,每道共 {npts} 个数据')
#data_list=np.array(data_list[0:1:1][:])

def wigb(data, scal:int=0, x=0, z:list=[], a_max:float=0):
    nx=len(data)    # x轴
    nz=len(data[0]) # y轴
    #print(f'nx{nx},nz{nz}')
    trace_max = [np.max(np.abs(data[:,i])) for i in range(nx)]
    #print(f'trace_max {trace_max}')
    if a_max==0:
        a_max=np.mean(trace_max)
    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(f"错误：道数必须大于 1 道！请重新选择数据！")
        return

    # 计算向量x中相邻元素的平均间距，并将其存储在变量"dx"中
    dx1 = np.abs(x[1:nx] - x[0:nx - 1])
    dx = np.median(dx1)

    # 计算向量z中相邻元素的差值（dz）
    dz = z[1] - z[0]

    # 存储矩阵"data"的最大值和最小值分别在变量x_max和x_min中
    x_max = np.max(data)
    x_min = np.min(data)

    # 将矩阵"data"乘以dx/a_max来缩放它，然后再乘以缩放因子"scal"进行调整
    k = dx / a_max
    data = data * k * scal

    # 打印有关数据范围和绘制的最大值的信息
    print(f'数据范围为：{x_min}到{x_max},画图最大值为：{a_max}\n')

    # 设置绘图的显示范围和属性
    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()

    z_start = z[0]
    z_end = z[-1]

    # 定义填充颜色、线条颜色和线宽
    fillcolor = [0, 0, 0]
    linecolor = [0, 0, 0]
    linewidth = 1.0

    # 遍历数据矩阵 a 的每一列，并绘制对应的 wiggle 形式
    for i in range(0,nx):
        # 如果当前追踪的最大振幅不为零，则为该追踪生成地震剖面图形
        if trace_max[i] != 0:
            tr=data[i][:]
            s=np.sign(tr)
            #print(f'tr={len(tr)},len(s)={len(s)}')
            # 寻找零点
            i1=[]
            for j in range(0,nz):
                if j==nz-1:
                    continue
                if s[j] != s[j+1]:
                    i1.append(j)
            # 根据之前找到的过零点位置列表i1来计算附加的零振幅位置
            i1=np.array(i1)
            # 列表i1为空，表示没有过零点，将变量zadd设置为空的NumPy数组np.array(())
            if len(i1)==0:
                zadd = np.array(())
            else:
                # i1非空，计算出附加的零振幅位置
                zadd = i1 + tr[i1] / (tr[i1] - tr[i1 + 1])
            # 创建一个形状与zadd相同的全零数组的零振幅位置的振幅值
            aadd = np.zeros(zadd.shape)

            # 处理正振幅和零振幅的位置信息
            zpos=np.where(tr>0)
            tmp=np.append(zpos,zadd)
            zz=np.sort(tmp)
            iz = np.argsort(tmp)
            aa = np.append(tr[zpos], aadd)
            iz = iz.astype(int)
            aa = aa[iz]

            # 假设存在正振幅，将变量a0设置为0，z0设置为1.00
            if tr[0] > 0:
                a0=0
                z0=1.00
            # 假设不存在正振幅，将变量a0设置为零,z0振幅位置列表zadd的第一个元素
            else:
                a0=0
                z0=zadd[0]

            # 假设存在正振幅，将变量a1设置为0，z1设置为地震数据的最后一个位置（nz）
            if tr[nz - 1] > 0:
                a1=0
                z1=nz
            # 假设不存在正振幅，将变量a1设置为0，z1设置为零振幅位置列表zadd中的最大值
            else:
                a1=0
                z1=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 = z_start + zz * dz - dz
            plt.fill(aa + x[i], zzz + 1, color=fillcolor)
            plt.plot(x[i] + [0, 0], [z_start, z_end], color=[1, 1, 1])
            plt.plot(x[i] + tr, z, color=linecolor, linewidth=linewidth)
        else:
            plt.plot([x[i], x[i]], [z_start, z_end], color=linecolor, linewidth=linewidth)
    plt.show()

wigb(data_list)