基于Python读取sgy文件并绘图

最新推荐文章于 2024-05-13 13:08:46 发布

MarsGril

最新推荐文章于 2024-05-13 13:08:46 发布

阅读量1.6k

点赞数 2

分类专栏： seismic processing Python 文章标签： python

本文链接：https://blog.csdn.net/wei_yanwen/article/details/132066067

版权

seismic processing 同时被 2 个专栏收录

5 篇文章 0 订阅

订阅专栏

Python

3 篇文章 0 订阅

订阅专栏



import os
import segyio
import random
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import matplotlib.colors as pcl
from matplotlib.colors import ListedColormap
plt.rcParams['font.sans-serif'] = ['Simhei']  #显示中文
plt.rcParams['axes.unicode_minus'] = False    #显示负号
from matplotlib.font_manager import FontProperties



class readsegy():
    def __init__(self, filename=None):
        self.filename = filename
    
    def readin_segy(self):
        with segyio.open(self.filename) as segyfile:
            # Memory map file for faster reading (especially if file is big...)
            segyfile.mmap()

            # Print binary header info
            print(f'====== Binary header ======\n{segyfile.bin}\n')
            # print(segyfile.bin[segyio.BinField.Traces])

            # Read headerword inline for trace 10
            print(f'====== Trace header for trace 10: ======\n{segyfile.header[10]}\n') #[segyio.TraceField.INLINE_3D])
            # print(segyfile.header[10][segyio.TraceField.INLINE_3D])

            # Print inline and crossline axis
            print(f'====== Inline ======\n {segyfile.xlines.min()} to {segyfile.xlines.max()} with step {segyfile.xlines[1] - segyfile.xlines[0]}')
            print(f'====== Crossline ======\n {segyfile.ilines.min()} to {segyfile.ilines.max()} with step {segyfile.ilines[1] - segyfile.ilines[0]}')
        
    def get_one_xline(self, lid=None):
        with segyio.open(self.filename) as segyfile:
            # Memory map file for faster reading (especially if file is big...)
            segyfile.mmap()
            # Read data along which (lid) xline
            data = segyfile.xline[segyfile.xlines[lid]]
            return data
    def get_one_iline(self, lid=None):
        with segyio.open(self.filename) as segyfile:
            # Memory map file for faster reading (especially if file is big...)
            segyfile.mmap()
            # Read data along which (lid) iline
            data = segyfile.iline[segyfile.ilines[lid]]
            return data 
        


# plotting

def get_plot_seismic(OUT=None, scale=0.1, title=None, ratio=16, nplot=None, saveflag=False):

    sc1 = OUT.mean() - scale*OUT.std()
    sc2 = OUT.mean() + scale*OUT.std()
    cN = pcl.Normalize(vmin=sc1, vmax=sc2)
    
    fig, ax = plt.subplots(figsize=(8,6), dpi=300)
    im = ax.imshow(OUT, cmap='seismic', norm=cN, interpolation='bicubic')
    ax.set_aspect(ratio)
    plt.ylabel(u'时间步', fontsize=10)#, fontproperties=SimHei) # y轴标题
    # plt.gca().invert_yaxis()
    plt.xlabel(u'水平距离', fontsize=10)#, fontproperties=SimHei) # x轴标题

    plt.title(title)
    nx = 512
    for ii in range(1, nplot):
        plt.axvline( int((ii)*nx), color='black') #, linewidth=2 )
    if saveflag:
        plt.savefig('currimage.pdf')
    return ax

依赖segyio和其他基本python库