使用python3读IBM格式segy文件并绘制变密度剖面图

# code for read IBM segy format seismic data files 
import matplotlib.pyplot as mp
import numpy as np 
import sys
import struct
import binascii

def ibm2ieee(ibm):
    if ibm == 0:
        return 0.0
    sign = ibm >> 31 & 0x01
    exponent = ibm >> 24 & 0x7f
    mantissa = (ibm & 0x00ffffff) / float(pow(2, 24))
    return (1 - 2 * sign) * mantissa * pow(16, exponent - 64)

fileName = "vp_marmousi-ii.segy"
nTrace = 13601
nSample = 2801
fSegy = open(fileName,"rb")
data = np.zeros((nSample,nTrace))
fSegy.seek(3600,0)    # 跳过3600字节卷头
for iTrace in range(nTrace):
    fSegy.seek(240,1) # 每一道都跳过240字节道头
    for iSample in range(nSample):
        tempValue = fSegy.read(4)
        hexValue = tempValue.hex()
        decValue = int(hexValue,16)
        unintValue = struct.unpack('>L',struct.pack('>I', decValue))[0] # < > 用于控制大小端转换
        floatValue = ibm2ieee(unintValue)
        data[iSample][iTrace] = floatValue
fSegy.close()
print(data)
mp.matshow(data, cmap=mp.cm.gray)
mp.show()```