Python读写SEGY数据的简单实现

目前只实现了二维规则数组的读写，未考虑大小端问题。

SEGY数据格式

3600字节文件头+240字节道头+道数据+240字节道头+道数据+......+240字节道头+道数据

本文只适用所有道数据采样点数、采样间隔一样的情况，此时只需要在文件头中写一些信息。文件头需要写的内容如下：

字节	格式	描述
3213-3214	short	每个道集的数据道数。
3217-3218	unsigned short	微秒（us）形式的采样间隔。
3221-3222	short	数据道采样点数。
3225-3226	unsigned short	数据采样格式编码。5=4字节IEEE浮点数

Python代码

使用seek移动文件读写指针，使用标准库struct实现数据类型转换，使用numpy构造数组。

写数据代码：

def save_segy(path, data, dt=1000):
    short2byte = lambda x: struct.pack('h', x)
    ushot2byte = lambda x: struct.pack('H', x)
    nx, nt = data.shape
    with open(path, 'wb') as f:
        f.seek(3212)
        f.write(short2byte(nx))
        f.seek(2, 1)
        f.write(ushot2byte(dt))
        f.seek(2, 1)
        f.write(ushot2byte(nt))
        f.seek(2, 1)
        f.write(short2byte(5))
        f.seek(3600, 0)
        for i in range(nx):
            f.seek(240, 1)
            f.write(struct.pack('f' * nt, *data[i]))

 读数据代码：

def read_segy(path):
    byte2short = lambda x: struct.unpack('h', x)
    byte2ushot = lambda x: struct.unpack('H', x)
    with open(path, 'rb') as f:
        f.seek(3212)
        nx = byte2short(f.read(2))[0]
        f.seek(6, 1)
        nt = byte2ushot(f.read(2))[0]
        f.seek(3600, 0)
        data = np.zeros((nx, nt), dtype='float32')
        for i in range(nx):
            f.seek(240, 1)
            dat = f.read(4 * nt)
            data[i] = np.array(struct.unpack('f' * nt, dat))
    return data