Python 的obspy是一个非常成熟的地震处理的库,详细使用已经有很多博主的说明。
最近遇到了一个问题就是如果使用Madagascar里面的算法进行正演或者反演的话,需要segy文件或者rsf文件。需要使用python创建通用的segy文件(可以使用sfsegyread读取的文件),在这里记录一下。
import numpy as np
import segyio
from obspy import read, Trace, Stream, UTCDateTime
from obspy.core import AttribDict
from obspy.io.segy.segy import SEGYTraceHeader, SEGYBinaryFileHeader
from obspy.io.segy.core import _read_segy, _write_segy
import sys
stream = Stream()
for _i in range(10):
data = np.random.ranf(1000)
data = np.require(data, dtype=np.float32)
trace = Trace(data)
trace.stats.delta = 0.001
'''
SEGY does not support microsecond precision! Any microseconds will be discarded
'''
trace.stats.starttime = UTCDateTime(2023, 11, 2, 22, 15, 0)
'''
Add trace header
'''
if not hasattr(trace.stats, 'segy.trace_header'):
trace.stats.segy = {}
trace.stats.segy.trace_header = SEGYTraceHeader()
trace.stats.segy.trace_header.trace_sequence_number_within_line = _i + 1
trace.stats.segy.trace_header.receiver_group_elevation = 444
stream.append(trace)
stream.stats = AttribDict()
stream.stats.textual_file_header = 'Textual Header!'
stream.stats.binary_file_header = SEGYBinaryFileHeader()
stream.stats.binary_file_header.trace_sorting_code = 5
_write_segy(stream, 'test1.segy', data_encoding=1, byteorder='>', textual_header_encoding=None) # (official function)
st1 = _read_segy("test1.segy") # (official function)
'''
default write and read in system —— the second methods to read and write segy file
'''
# stream.write("test.segy", format="SEGY", data_encoding=1, byteorder=sys.byteorder)
# st2 = read("TEST.segy")