代码
- 注意参数
- nTrace = 311
- nSample = 783
# 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 = ".\data\LX_SEGY005.segy"
nTrace = 311
nSample = 783
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()
