Dataset :
3D reflection seismic data is acquired offshore of southeast Japan as part of the Nankai Trough Seismogenic Zone Experiment (NanTroSEIZE) . This data can be used to study active accretionary prism processes.
Motivation :
The 3D seismic volume revealed complex interactions between active sedimentation and tectonics within multiple slope basins above the accretionary prism. We require xpensive specialized software packages to understand these interactions .
Implementation :
In this notebook we will be learning how to implement stratal slicing of the 3D volume and co-rendering of multiple attributes in python to better visualize our results. Stratal slicing allows volumetric attributes to be displayed in map view along an arbitrary geologic timeline(~30MB animated gif) by interpolating between interpreted geologic surfaces. This enhances the visibility of subtle changes in stratigraphic architecture through time. Co-rendering coherence on top of seismic amplitudes facilitates fault interpretation in both cross section and map view. This technique allowed us to confidently interpret faults near the limit of seismic resolution.The scientific python ecosystem proved to be an effective platform both for making publication-quality cross sections and for rapidly implementing state-of-the-art seismic visualization techniques.
!pip install geoprobeCollecting geoprobe
Downloading https://files.pythonhosted.org/packages/c5/ff/cd775f26011164c75362ec772b062ada6924b7c53b98355e82ab385cb9cc/geoprobe-0.4.0.tar.gz
Requirement already satisfied: numpy in /opt/conda/lib/python3.6/site-packages (from geoprobe) (1.18.1)
Requirement already satisfied: six in /opt/conda/lib/python3.6/site-packages (from geoprobe) (1.14.0)
Building wheels for collected packages: geoprobe
Building wheel for geoprobe (setup.py) ... - \ done
Created wheel for geoprobe: filename=geoprobe-0.4.0-cp36-none-any.whl size=41055 sha256=0de508ec1c9a65489e072a1eae6194ae72569fdae43b157f15898ab51facd251
Stored in directory: /root/.cache/pip/wheels/0a/4b/3f/e31bd6626dca606594f1c17532dc0267ebf69e8ce391c08a1a
Successfully built geoprobe
Installing collected packages: geoprobe
Successfully installed geoprobe-0.4.0
#importing the libraries
import geoprobe
import matplotlib.pyplot as plt
from matplotlib.widgets import Slider
import numpy as np
from scipy.ndimage import gaussian_filter1d#Utilities Functions for seismic visualizations
def stack_bands(r, g, b):
"""Stacks three 3D volumes into a 4D RGB volume along the last axis."""
slice_ = np.s_[..., None]
return np.concatenate([r[slice_], g[slice_], b[slice_]], axis=r.ndim)
class Explorer(object):
"""A simple interactive slicer that "pans" a series of 2D slices along a
given axis of a 3D volume. Additional kwargs are passed to "imshow"."""
def __init__(self, data, axis=2, **kwargs):
self.data = data
self.axis = axis
self.start, self.stop = 0, data.shape[self.axis] - 1
kwargs['cmap'] = kwargs.get('cmap', 'gray_r')
self.fig, self.ax = plt.subplots()
self.im = self.ax.imshow(self[self.start], **kwargs)
self.ax.axis('off')
self.sliderax = self.fig.add_axes([0.2, 0.02, 0.6, 0.03])
self.slider = Slider(self.sliderax, 'Z-slice', self.start, self.stop,
valinit=self.start, valfmt='%i')
self.slider.on_changed(self.update)
def __getitem__(self, i):
slices = self.data.ndim * [slice(None)]
slices[self.axis] = i
return self.data[slices].swapaxes(0, 1)
def update(self, val):
dat = self[int(val)]
self.im.set(data=dat, clim=[dat.min(), dat.max()])
self.fig.canvas.draw()
def show(self):
plt.show()
def grid(x, y, z, dx=1, dy=1, extent=None, nodata=0):
"""Turn unordered but regularly spaced (i.e. alreadly gridded) points into
a regular 2D grid. "z" can be either a 1d array of numbers or a 2d array
of r,g,b colors."""
if extent is None:
extent = [x.min(), y.min(), x.max(), y.max()]
xmin, ymin, xmax, ymax = extent
nrows = (ymax - ymin) // dy + 1
ncols = (xmax - xmin) // dx + 1
i = ((y - ymin) / dy).astype(int)
j = ((x - xmin) / dx).astype(int)
if z.ndim > 1:
output = nodata * np.ones((nrows, ncols, z.shape[1]), z.dtype)
output[i, j, :] = z
else:
output = nodata * np.ones((nrows, ncols), z.dtype)
output[i, j] = z
return output