Simple interactive viewer for seismic data built with python+streamlit+segyio

seismic_viewer_app_demo.py

# seismic_viewer_app_demo
#--------------------------------------------------------------------
# 
# testing streamlit for a seismic viewer app
# using Near stack subcube from Per Avseth's QSI dataset
#
# to run: 
# $ streamlit run seismic_viewer_app_demo.py
#--------------------------------------------------------------------

import numpy as np
import matplotlib.pyplot as plt
import streamlit as st
import os
import requests
from tempfile import NamedTemporaryFile
import segyio

@st.cache_data
def st_load_seismic():
    remote_file = 'https://github.com/aadm/geophysical_notes/raw/master/3d_nearstack.sgy'
    response = requests.get(remote_file)

    # raise an HTTPError if the download fails
    if response.raise_for_status() is False:
        st.print('error while downloading {}'.format(remote_file))
        return None, None, None, None

    # create a temporary file to save the downloaded SEGY file
    with NamedTemporaryFile(delete=False) as tmp:
        tmp.write(response.content)
        local_file = tmp.name

    with segyio.open(local_file, iline=41, xline=21) as f:
        ntraces = f.tracecount
        n_inl = f.ilines.size
        n_crl = f.xlines.size
        sr = segyio.tools.dt(f)/1e3
        nsamples = f.samples.size
        twt = f.samples
        data = segyio.tools.cube(f)

    os.remove(local_file) # remove local files

    seis = data.T
    z, crl, inl = seis.shape
    return seis, nsamples, n_crl, n_inl

@st.cache_data
def st_get_section(data, line_type, line):
    if line_type == 'inline':
        section = data[:,:,line]
    else:
        section = data[:,line,:]
    return section


@st.cache_data
def st_plot_section(data, line_type, line, z,  ztop=None, zbot=None,
               cmap='bwr_r', clip=None, interpolation='spline16',
               add_colorbar=True, square=True):
    ''' modified from squit.plot.model_seis() '''
    n_traces = data.shape[1]

    if ztop is None: ztop = z.min()
    if zbot is None: zbot = z.max()

    # calculate upper/lower clips value for color palette
    if clip is None:
        clip = np.nanpercentile(data, 100)
        c1, c2 = -clip, +clip
    elif isinstance(clip, list) or isinstance(clip, tuple):
        c1, c2 = clip
    else:
        c1, c2 = -clip, +clip

    # make plot
    fig, ax = plt.subplots(ncols=1, constrained_layout=True)
    extents = [0, n_traces, z[-1], z[0]]

    opt = dict(extent=extents, cmap=cmap,
                vmin=c1, vmax=c2,
                aspect='auto', interpolation=interpolation)
    im = ax.imshow(data, **opt)
    if add_colorbar:
        cbar = plt.colorbar(im, ax=ax, orientation='horizontal')
        cbar.ax.tick_params(labelsize='x-small')

    title = '{} {}'.format(line_type, line)
    ax.set_title(title, fontsize='large')
    ax.tick_params(axis='both', labelsize='x-small')
    ax.set_ylim(ztop, zbot)
    ax.invert_yaxis()
    if square:
        ax.set_aspect('equal', 'box')

    opt = dict(ztop=ztop, zbot=zbot, cmap=colormap, clip=clip)   
    st.pyplot(fig=fig, use_container_width=True)

st.set_page_config(page_title='Simple Seismic Viewer', layout="centered")

colormap_list = [
    'RdBu', 'RdGy', 'bwr_r', 'coolwarm_r', 'seismic_r',
    'gray_r' or 'gist_gray_r'
    ]

#--- BEGIN APP
st.header('Simple Seismic Viewer')

seis, z, crl, inl = st_load_seismic()

widg_load = st.columns(2)
with widg_load[0]:
    line_type = st.selectbox('inline or crossline', ('inline', 'crossline'))
with widg_load[1]:
    min_line = 0
    if line_type == 'inline':
        max_line = inl
    else:
        max_line = crl
    opt_line = dict(min_value=min_line, max_value=max_line, format='%i')
    line = st.slider('{} #'.format(line_type), **opt_line)

section = st_get_section(seis, line_type, line)

z_vect = np.arange(z)
clip000, clip100 = np.percentile(np.abs(section), (1, 100))
clip099 = np.percentile(np.abs(section), 99)

widg_a = st.columns(3)
with widg_a[0]:
    colormap = st.selectbox('colormap', colormap_list, label_visibility="collapsed",placeholder='colormap')
with widg_a[1]:
    opt_clip = dict(min_value=clip000, max_value=clip100, value=clip099)
    clip = st.number_input('clip', **opt_clip,  label_visibility="collapsed")
with widg_a[2]:
    square = st.checkbox('1:1 aspect ratio')

widg_b = st.columns(2)
opt_z = dict(min_value=z_vect.min(), max_value=z_vect.max(), step=10, format='%i')
with widg_b[0]:
    ztop = st.slider('top window', value=z_vect.min(), **opt_z)
with widg_b[1]:
    zbot = st.slider('bot window', value=z_vect.max(), **opt_z)
if zbot <= ztop:
    zbot = ztop+100

opt = dict(ztop=ztop, zbot=zbot, cmap=colormap, clip=clip, square=square)
st_plot_section(section, line_type, line, z_vect, **opt)