最近在设计开发地震数据服务,需要对SEGY格式的地震数据进行读取,并解析SEGY卷头和道头数据,以及实际数据,之前一直采取C++、C#或java解析处理,本文主要利用python3读取SEGY格式,并对其进行解析及统计量分析,体会到了python在数据读取、处理和分析方面的强大之处和便利性。代码送上。
1.解析卷头的400字节二进制数据都有啥。经过大量分析,这部分数据主要以2字节或4字节整形存放,代码如下:
#coding=utf-8
import sys
import math
from array import array
import numpy as np
startPos=1472974
arr = array('L') # 长整型 4字节
# 读取卷头400字节的二进制卷头数据,并输出打印
def read400(fileName,startPos,len):
try:
# with open(sys.argv[1], 'rb') as f:
with open(fileName, 'rb') as f:
f.seek(startPos)
arr.fromfile(f, len)
arr.byteswap() # 字节反序一下
f.close()
except (IOError,Exception):
print('usage: scriptname segyfilename')
sys.exit(1)
# 输出道console
for i in arr:
print(i)
if __name__ == '__main__':
filename = "D:\project\PythonProject\aaa.sgy"
read400(filename,3200,400)
2.解析SEGY中指定任意道的实际数据,并计算其统计量,包括min/max/mean/rms/sum等。一般数据体中的样本数据是单精度浮点类型,占用4字节,高低位存储方式解析出来查看数据即可,如果存在异常,就采用arr2.byteswap() ,将字节反序即可,代码如下:
startPos=1472974
bufferLen=42345
samples=3000
max=min=mean=rms=0.0
sum=0.0
sumAvg=0.0
# 读取segy数据中指定任意道的实际数据,并计算其统计量,包括min/max/mean/rms等。
def readSegyBytes(fileName, traceid):
global max, min, sum,sumAvg
print('readSegyBytes--'+str(traceid))
arr2 = array('f') # 浮点类型 4字节
startPos=3600+(traceid-1)*samples*4+traceid*240 #第i道,注意是3000个样本点,每个占4字节,共计12000字节。
bufferLen=samples #3000个样本点
try:
with open(fileName, 'rb') as f:
f.seek(startPos)
arr2.fromfile(f, bufferLen)
arr2.byteswap() #字节反序一下
f.close()
except (IOError, Exception):
print('usage: scriptname segyfilename')
sys.exit(1)
sampleData = np.array(arr2);
print('mean='+str(sampleData.mean()))
print('min=' + str(sampleData.min()))
print('max=' + str(sampleData.max()))
print('rms=' + str(sampleData.std()))
print('var=' + str(sampleData.var()))
print('sum=' + str(sampleData.sum()))
if __name__ == '__main__':
traceid=4
readSegyBytes(filename,traceid)