# -*- coding: utf-8 -*- show_zgy_section.py 读取 ZGY(Schlumberger/OpenZGY)格式三维地震数据体并显示 inline / crossline / 时间(深度)切片剖面。 依赖安装 pip install openzgy numpy matplotlib 用法示例 python show_zgy_section.py data.zgy # 打印数据体信息 python show_zgy_section.py data.zgy --inline 500 # 显示某条 inline 剖面 python show_zgy_section.py data.zgy --crossline 300 # 显示某条 crossline 剖面 python show_zgy_section.py data.zgy --timeslice 2000 # 显示某时间切片 python show_zgy_section.py data.zgy --inline 500 --cmap seismic --save inline500.png import argparse import sys import numpy as np import matplotlib.pyplot as plt def open_reader(zgy_file): 打开 ZGY 文件返回 reader 对象失败时给出中文提示。 try: from openzgy.api import ZgyReader except ImportError: sys.exit(缺少依赖库 openzgy请先执行pip install openzgy) try: return ZgyReader(zgy_file) except Exception as e: sys.exit(f无法打开 ZGY 文件 {zgy_file}: {e}) def _geometry(reader): 返回 (il_start, il_step, xl_start, xl_step, z_start, z_step, z_unit)。 il0, xl0 reader.annotstart # 起始真实坐标 (inline, crossline) dil, dxl reader.annotinc # 步长 z0, dz, zunit reader.zstart, reader.zinc, reader.zunitname or ms return il0, dil, xl0, dxl, z0, dz, zunit def print_info(reader): 打印数据体基本 geometry 信息。 il0, dil, xl0, dxl, z0, dz, zunit _geometry(reader) ni, nx, nz reader.size print( * 60) print(f数据体大小 (inline, crossline, samples): {ni} x {nx} x {nz}) print(fInline 范围: {il0} ~ {il0 (ni - 1) * dil} (step {dil})) print(fCrossline 范围: {xl0} ~ {xl0 (nx - 1) * dxl} (step {dxl})) print(fZ 范围: {z0} ~ {z0 (nz - 1) * dz} (step {dz}, 单位 {zunit})) print( * 60) def _real_to_index(real, start, step): 把真实道号/采样值换算成数组下标。 idx int(round((real - start) / step)) return idx def read_section(reader, kind, value): 读取一个剖面返回 (data, x轴标签, y轴标签)。 kind: inline / crossline / timeslice il0, dil, xl0, dxl, z0, dz, zunit _geometry(reader) ni, nx, nz reader.size if kind inline: idx _real_to_index(value, il0, dil) if not (0 idx ni): sys.exit(fInline {value} 超出范围 (数组下标 {idx}共 {ni} 条)) # read(start, data): start 为起点下标data 为预先分配的 float32 缓冲 buf np.zeros((1, nx, nz), dtypenp.float32) reader.read((idx, 0, 0), buf) data buf[0, :, :].T # (nz, nx) x xl0 np.arange(nx) * dxl y z0 np.arange(nz) * dz return data, x, y, fInline {value} if kind crossline: idx _real_to_index(value, xl0, dxl) if not (0 idx nx): sys.exit(fCrossline {value} 超出范围 (数组下标 {idx}共 {nx} 条)) buf np.zeros((ni, 1, nz), dtypenp.float32) reader.read((0, idx, 0), buf) data buf[:, 0, :].T # (nz, ni) x il0 np.arange(ni) * dil y z0 np.arange(nz) * dz return data, x, y, fCrossline {value} if kind timeslice: idx _real_to_index(value, z0, dz) if not (0 idx nz): sys.exit(f时间切片 {value} 超出范围 (数组下标 {idx}共 {nz} 个采样)) buf np.zeros((ni, nx, 1), dtypenp.float32) reader.read((0, 0, idx), buf) data buf[:, :, 0] # (ni, nx) x xl0 np.arange(nx) * dxl y il0 np.arange(ni) * dil return data, x, y, fTime slice {value} {zunit} sys.exit(f未知剖面类型: {kind}) def plot_section(data, x, y, title, xlabel, ylabel, cmap, saveNone, clip99): 绘制二维剖面imshow 色标clip 为振幅截断百分位压制异常值。 vmax np.percentile(np.abs(data[np.isfinite(data)]), clip) if np.isfinite(data).any() else 1 vmax vmax if vmax 0 else 1.0 fig, ax plt.subplots(figsize(12, 6)) extent [x[0], x[-1], y[-1], y[0]] # y 轴默认上小下大深度向下 im ax.imshow(data, cmapcmap, aspectauto, extentextent, vmin-vmax, vmaxvmax, interpolationbilinear) ax.set_title(title) ax.set_xlabel(xlabel) ax.set_ylabel(ylabel) fig.colorbar(im, axax, labelAmplitude) fig.tight_layout() if save: fig.savefig(save, dpi200) print(f剖面已保存到: {save}) else: plt.show() def main(): p argparse.ArgumentParser(description显示 ZGY 地震体剖面) p.add_argument(zgy_file, helpZGY 文件路径) p.add_argument(--inline, typefloat, defaultNone, help要显示的 inline 号) p.add_argument(--crossline, typefloat, defaultNone, help要显示的 crossline 号) p.add_argument(--timeslice, typefloat, defaultNone, help要显示的时间/深度切片值) p.add_argument(--cmap, defaultseismic, helpcolormap默认 seismic) p.add_argument(--save, defaultNone, help保存为图片文件(如 out.png)不填则弹出窗口) args p.parse_args() reader open_reader(args.zgy_file) print_info(reader) _, _, _, _, _, _, zunit _geometry(reader) sections [ (inline, args.inline, Crossline, zunit), (crossline, args.crossline, Inline, zunit), (timeslice, args.timeslice, Crossline, Inline), ] shown False n_sections sum(1 for _, v, _, _ in sections if v is not None) for kind, value, xlab, ylab in sections: if value is None: continue data, x, y, title read_section(reader, kind, value) # 多个剖面同时保存时自动在文件名中加剖面标识避免互相覆盖 save args.save if save and n_sections 1: base, dot, ext save.rpartition(.) base base or save ext dot ext if dot else save f{base}_{kind}{int(value)}{ext} plot_section(data, x, y, title, xlab, ylab, args.cmap, save, clip99) shown True if not shown: print(未指定剖面使用 --inline / --crossline / --timeslice 之一来显示剖面。) print(提示可先只传文件名查看数据体范围再选取合适的剖面号。) if __name__ __main__: main() zgy地震数据体的显示
# -*- coding: utf-8 -*- show_zgy_section.py 读取 ZGY(Schlumberger/OpenZGY)格式三维地震数据体并显示 inline / crossline / 时间(深度)切片剖面。 依赖安装 pip install openzgy numpy matplotlib 用法示例 python show_zgy_section.py data.zgy # 打印数据体信息 python show_zgy_section.py data.zgy --inline 500 # 显示某条 inline 剖面 python show_zgy_section.py data.zgy --crossline 300 # 显示某条 crossline 剖面 python show_zgy_section.py data.zgy --timeslice 2000 # 显示某时间切片 python show_zgy_section.py data.zgy --inline 500 --cmap seismic --save inline500.png import argparse import sys import numpy as np import matplotlib.pyplot as plt def open_reader(zgy_file): 打开 ZGY 文件返回 reader 对象失败时给出中文提示。 try: from openzgy.api import ZgyReader except ImportError: sys.exit(缺少依赖库 openzgy请先执行pip install openzgy) try: return ZgyReader(zgy_file) except Exception as e: sys.exit(f无法打开 ZGY 文件 {zgy_file}: {e}) def _geometry(reader): 返回 (il_start, il_step, xl_start, xl_step, z_start, z_step, z_unit)。 il0, xl0 reader.annotstart # 起始真实坐标 (inline, crossline) dil, dxl reader.annotinc # 步长 z0, dz, zunit reader.zstart, reader.zinc, reader.zunitname or ms return il0, dil, xl0, dxl, z0, dz, zunit def print_info(reader): 打印数据体基本 geometry 信息。 il0, dil, xl0, dxl, z0, dz, zunit _geometry(reader) ni, nx, nz reader.size print( * 60) print(f数据体大小 (inline, crossline, samples): {ni} x {nx} x {nz}) print(fInline 范围: {il0} ~ {il0 (ni - 1) * dil} (step {dil})) print(fCrossline 范围: {xl0} ~ {xl0 (nx - 1) * dxl} (step {dxl})) print(fZ 范围: {z0} ~ {z0 (nz - 1) * dz} (step {dz}, 单位 {zunit})) print( * 60) def _real_to_index(real, start, step): 把真实道号/采样值换算成数组下标。 idx int(round((real - start) / step)) return idx def read_section(reader, kind, value): 读取一个剖面返回 (data, x轴标签, y轴标签)。 kind: inline / crossline / timeslice il0, dil, xl0, dxl, z0, dz, zunit _geometry(reader) ni, nx, nz reader.size if kind inline: idx _real_to_index(value, il0, dil) if not (0 idx ni): sys.exit(fInline {value} 超出范围 (数组下标 {idx}共 {ni} 条)) # read(start, data): start 为起点下标data 为预先分配的 float32 缓冲 buf np.zeros((1, nx, nz), dtypenp.float32) reader.read((idx, 0, 0), buf) data buf[0, :, :].T # (nz, nx) x xl0 np.arange(nx) * dxl y z0 np.arange(nz) * dz return data, x, y, fInline {value} if kind crossline: idx _real_to_index(value, xl0, dxl) if not (0 idx nx): sys.exit(fCrossline {value} 超出范围 (数组下标 {idx}共 {nx} 条)) buf np.zeros((ni, 1, nz), dtypenp.float32) reader.read((0, idx, 0), buf) data buf[:, 0, :].T # (nz, ni) x il0 np.arange(ni) * dil y z0 np.arange(nz) * dz return data, x, y, fCrossline {value} if kind timeslice: idx _real_to_index(value, z0, dz) if not (0 idx nz): sys.exit(f时间切片 {value} 超出范围 (数组下标 {idx}共 {nz} 个采样)) buf np.zeros((ni, nx, 1), dtypenp.float32) reader.read((0, 0, idx), buf) data buf[:, :, 0] # (ni, nx) x xl0 np.arange(nx) * dxl y il0 np.arange(ni) * dil return data, x, y, fTime slice {value} {zunit} sys.exit(f未知剖面类型: {kind}) def plot_section(data, x, y, title, xlabel, ylabel, cmap, saveNone, clip99): 绘制二维剖面imshow 色标clip 为振幅截断百分位压制异常值。 vmax np.percentile(np.abs(data[np.isfinite(data)]), clip) if np.isfinite(data).any() else 1 vmax vmax if vmax 0 else 1.0 fig, ax plt.subplots(figsize(12, 6)) extent [x[0], x[-1], y[-1], y[0]] # y 轴默认上小下大深度向下 im ax.imshow(data, cmapcmap, aspectauto, extentextent, vmin-vmax, vmaxvmax, interpolationbilinear) ax.set_title(title) ax.set_xlabel(xlabel) ax.set_ylabel(ylabel) fig.colorbar(im, axax, labelAmplitude) fig.tight_layout() if save: fig.savefig(save, dpi200) print(f剖面已保存到: {save}) else: plt.show() def main(): p argparse.ArgumentParser(description显示 ZGY 地震体剖面) p.add_argument(zgy_file, helpZGY 文件路径) p.add_argument(--inline, typefloat, defaultNone, help要显示的 inline 号) p.add_argument(--crossline, typefloat, defaultNone, help要显示的 crossline 号) p.add_argument(--timeslice, typefloat, defaultNone, help要显示的时间/深度切片值) p.add_argument(--cmap, defaultseismic, helpcolormap默认 seismic) p.add_argument(--save, defaultNone, help保存为图片文件(如 out.png)不填则弹出窗口) args p.parse_args() reader open_reader(args.zgy_file) print_info(reader) _, _, _, _, _, _, zunit _geometry(reader) sections [ (inline, args.inline, Crossline, zunit), (crossline, args.crossline, Inline, zunit), (timeslice, args.timeslice, Crossline, Inline), ] shown False n_sections sum(1 for _, v, _, _ in sections if v is not None) for kind, value, xlab, ylab in sections: if value is None: continue data, x, y, title read_section(reader, kind, value) # 多个剖面同时保存时自动在文件名中加剖面标识避免互相覆盖 save args.save if save and n_sections 1: base, dot, ext save.rpartition(.) base base or save ext dot ext if dot else save f{base}_{kind}{int(value)}{ext} plot_section(data, x, y, title, xlab, ylab, args.cmap, save, clip99) shown True if not shown: print(未指定剖面使用 --inline / --crossline / --timeslice 之一来显示剖面。) print(提示可先只传文件名查看数据体范围再选取合适的剖面号。) if __name__ __main__: main()