Note
Get this example as a Jupyter Notebook ipynb
[ ]:
SEGY-SAK Basics¶
segysak offers a number of utilities to create and load seismic data using xarray
and segyio
. In general segysak uses xarray.Dataset
to store the data and provides and interface to additional seismic specific functionality by adding the .seis
and .seisio
names-spaces to an xarray.Dataset
(just dataset
from now on). That sounds complicated but let us walk through some examples together.
Creating empty 3D geometry¶
In segysak we use the term seisnc
to refer to a dataset
which is compatible with segysak’s functionality and which has the additional names spaces registered with xarray
, for all intensive purposes it is an xarray.Dataset
but with defined dimensions and coordinates and some extended functionality. The seisnc
dimensions are defined depending on what type of seismic it is (2D, 3D, gathers, etc.)
To create an empty 3D instance of seisnc
use the create3d_dataset
. The function creates a new seisnc
based upon definitions for the dimensions, iline
numbering, xline
numbering and the vertical sampling.
[1]:
import matplotlib.pyplot as plt
%matplotlib inline
from segysak import create3d_dataset
# dims iline, xline, vert
dims = (10, 5, 1000)
new_seisnc = create3d_dataset(
dims,
first_iline=1,
iline_step=2,
first_xline=10,
xline_step=1,
first_sample=0,
sample_rate=1,
vert_domain="TWT",
)
new_seisnc
[1]:
<xarray.Dataset> Dimensions: (iline: 10, twt: 1000, xline: 5) Coordinates: * iline (iline) int64 1 3 5 7 9 11 13 15 17 19 * xline (xline) int64 10 11 12 13 14 * twt (twt) int64 0 1 2 3 4 5 6 7 8 ... 992 993 994 995 996 997 998 999 Data variables: *empty* Attributes: ns: None sample_rate: 1 text: SEGY-SAK Create 3D Dataset measurement_system: None d3_domain: TWT epsg: None corner_points: [(1, 10), (19, 10), (19, 14), (1, 14), (1, 10)] corner_points_xy: None source_file: None srd: None datatype: None percentiles: None coord_scalar: None
- iline: 10
- twt: 1000
- xline: 5
- iline(iline)int641 3 5 7 9 11 13 15 17 19
array([ 1, 3, 5, 7, 9, 11, 13, 15, 17, 19])
- xline(xline)int6410 11 12 13 14
array([10, 11, 12, 13, 14])
- twt(twt)int640 1 2 3 4 5 ... 995 996 997 998 999
array([ 0, 1, 2, ..., 997, 998, 999])
- ns :
- None
- sample_rate :
- 1
- text :
- SEGY-SAK Create 3D Dataset
- measurement_system :
- None
- d3_domain :
- TWT
- epsg :
- None
- corner_points :
- [(1, 10), (19, 10), (19, 14), (1, 14), (1, 10)]
- corner_points_xy :
- None
- source_file :
- None
- srd :
- None
- datatype :
- None
- percentiles :
- None
- coord_scalar :
- None
Dimension based selection and transformation¶
As you can see form the print out of the previous cell, we have three dimensions in this dataset. They are iline
, xline
and twt
(although the order, number and names might change depending on the make up of our volume). The ordering isn’t import to xarray
because it uses labels, and accessing data is done using these labels rather than indexing directly into the data like numpy
. xarray
also makes it further convient by allowing us to select based on the dimension values
using the .sel
method with tools for selecting nearest or ranges as well. If necessary you can also select by index using the .isel
method.
[2]:
# select iline value 9
xarray_selection = new_seisnc.sel(iline=9)
# select xline value 12
xarray_selection = new_seisnc.sel(xline=12)
# select iline and xline intersection point
xarray_selection = new_seisnc.sel(iline=9, xline=12)
# key error
# xarray_selection = new_seisnc.sel(twt=8.5)
# select nearest twt slice
xarray_selection = new_seisnc.sel(twt=8.5, method="nearest")
# select a range
xarray_selection = new_seisnc.sel(iline=[9, 11, 13])
# select a subcube
# also slices can be used to select ranges as if they were indices!
xarray_selection = new_seisnc.sel(iline=slice(9, 13), xline=[10, 11, 12])
# index selection principles are similar
xarray_selection = new_seisnc.sel(iline=slice(1, 4))
# putting it altogether to extract a sub-cropped horizon slice at odd interval
xarray_selection = new_seisnc.sel(
twt=8.5, iline=new_seisnc.iline[:4], xline=[10, 12], method="nearest"
)
xarray_selection
[2]:
<xarray.Dataset> Dimensions: (iline: 4, xline: 2) Coordinates: * iline (iline) int64 1 3 5 7 * xline (xline) int64 10 12 twt int64 9 Data variables: *empty* Attributes: ns: None sample_rate: 1 text: SEGY-SAK Create 3D Dataset measurement_system: None d3_domain: TWT epsg: None corner_points: [(1, 10), (19, 10), (19, 14), (1, 14), (1, 10)] corner_points_xy: None source_file: None srd: None datatype: None percentiles: None coord_scalar: None
- iline: 4
- xline: 2
- iline(iline)int641 3 5 7
array([1, 3, 5, 7])
- xline(xline)int6410 12
array([10, 12])
- twt()int649
array(9)
- ns :
- None
- sample_rate :
- 1
- text :
- SEGY-SAK Create 3D Dataset
- measurement_system :
- None
- d3_domain :
- TWT
- epsg :
- None
- corner_points :
- [(1, 10), (19, 10), (19, 14), (1, 14), (1, 10)]
- corner_points_xy :
- None
- source_file :
- None
- srd :
- None
- datatype :
- None
- percentiles :
- None
- coord_scalar :
- None
Coordinates Selection¶
Usually for seismic the X and Y coordinates labelled cdp_x
and cdp_y
in seisnc are rotated and scaled relative to the grid geometry and now seisnc dimensions iline
, xline
and twt
. For xarray
this means you cannot use the .sel
and .isel
methods to select data for cdp_x
and cdp_y
. segysak is developing more natural interfaces to access data using X and Y coordinates and this is available thorugh the seisnc.seis
namespace, covered in other examples.
Adding data to an empty seisnc¶
Because xarray
needs to understand the dimensions of any data you assign it must be explicitly communicated either via labels or creating an xarray.DataArray
first.
[3]:
# get the dims
dims = new_seisnc.dims
print(dims)
dkeys = ("xline", "iline", "twt")
dsize = [dims[key] for key in dkeys]
print("keys:", dkeys, "sizes:", dsize)
Frozen(SortedKeysDict({'iline': 10, 'xline': 5, 'twt': 1000}))
keys: ('xline', 'iline', 'twt') sizes: [5, 10, 1000]
[4]:
import numpy as np
# create some to data the dimension shapes
xline_, iline_, twt_ = np.meshgrid(new_seisnc.iline, new_seisnc.xline, new_seisnc.twt)
data = np.sin(twt_ / 100) + 2 * iline_ * np.cos(xline_ / 20 + 10)
# assign the data to dataset by passing in a tuple of the dimension keys and the new data
new_seisnc["data"] = (dkeys, data)
fig, axs = plt.subplots(ncols=3, figsize=(15, 5))
# axes are wrong for seismic
new_seisnc.data.sel(iline=7).plot(ax=axs[0])
# rotate the cube
new_seisnc.data.transpose("twt", "iline", "xline").sel(iline=7).plot(
ax=axs[1], yincrease=False
)
# xline is the same?
new_seisnc.data.transpose("twt", "iline", "xline").isel(xline=2).plot(
ax=axs[2], yincrease=False
)
[4]:
<matplotlib.collections.QuadMesh at 0x7f4b2a1355e0>
Other Useful Methods¶
Accessing the coordinates or data as a numpy
array.
[5]:
# Data can be dropped out of xarray but you then need to manage the coordinates and the dimension
# labels
new_seisnc.iline.values
[5]:
array([ 1, 3, 5, 7, 9, 11, 13, 15, 17, 19])
xarray.Dataset
and DataArray
have lots of great built in methods in addition to plot.
[6]:
print(new_seisnc.data.mean())
<xarray.DataArray 'data' ()>
array(-10.76364169)
[7]:
print(new_seisnc.data.max())
<xarray.DataArray 'data' ()>
array(0.08882943)
numpy
functions also work on xarray
objects but this returns a new DataArray
not a ndarray
.
[8]:
print(np.sum(new_seisnc.data))
<xarray.DataArray 'data' ()>
array(-538182.08430941)
With xarray
you can apply operations along 1 or more dimensions to reduce the dataset. This could be useful for colapsing gathers for example by applying the mean along the offset
dimension. Here we combine a numpy
operation abs
which returns an DataArray
and then sum along the time dimension to create a grid without the time dimension. Along with using masks this is fundamental building block for performing horizonal sculpting.
[9]:
map_data = np.abs(new_seisnc.data).sum(dim="twt")
img = map_data.plot()
Sometimes we need to modify the dimensions because they were read wrong or do scale them. Modify your dimension from the seisnc and then put it back using assign_coords
.
[10]:
new_seisnc.assign_coords(iline=new_seisnc.iline * 10, twt=new_seisnc.twt + 1500)
[10]:
<xarray.Dataset> Dimensions: (iline: 10, twt: 1000, xline: 5) Coordinates: * iline (iline) int64 10 30 50 70 90 110 130 150 170 190 * xline (xline) int64 10 11 12 13 14 * twt (twt) int64 1500 1501 1502 1503 1504 ... 2495 2496 2497 2498 2499 Data variables: data (xline, iline, twt) float64 -16.22 -16.21 -16.2 ... -1.803 -1.811 Attributes: ns: None sample_rate: 1 text: SEGY-SAK Create 3D Dataset measurement_system: None d3_domain: TWT epsg: None corner_points: [(1, 10), (19, 10), (19, 14), (1, 14), (1, 10)] corner_points_xy: None source_file: None srd: None datatype: None percentiles: None coord_scalar: None
- iline: 10
- twt: 1000
- xline: 5
- iline(iline)int6410 30 50 70 90 110 130 150 170 190
array([ 10, 30, 50, 70, 90, 110, 130, 150, 170, 190])
- xline(xline)int6410 11 12 13 14
array([10, 11, 12, 13, 14])
- twt(twt)int641500 1501 1502 ... 2497 2498 2499
array([1500, 1501, 1502, ..., 2497, 2498, 2499])
- data(xline, iline, twt)float64-16.22 -16.21 ... -1.803 -1.811
array([[[-16.2166637 , -16.20666387, -16.19666503, ..., -16.73527165, -16.7437957 , -16.75226703], [-14.96704323, -14.9570434 , -14.94704456, ..., -15.48565118, -15.49417523, -15.50264656], [-13.56787701, -13.55787718, -13.54787834, ..., -14.08648496, -14.09500901, -14.10348034], ..., [ -4.86226845, -4.85226862, -4.84226978, ..., -5.3808764 , -5.38940045, -5.39787179], [ -2.90121334, -2.8912135 , -2.88121467, ..., -3.41982129, -3.42834533, -3.43681667], [ -0.91117026, -0.90117042, -0.89117159, ..., -1.42977821, -1.43830225, -1.44677359]], [[-17.83833007, -17.82833024, -17.8183314 , ..., -18.35693802, -18.36546207, -18.3739334 ], [-16.46374755, -16.45374772, -16.44374888, ..., -16.9823555 , -16.99087955, -16.99935089], [-14.92466471, -14.91466488, -14.90466604, ..., -15.44327266, -15.45179671, -15.46026805], ... [ -6.32094899, -6.31094915, -6.30095032, ..., -6.83955694, -6.84808099, -6.85655232], [ -3.77157734, -3.7615775 , -3.75157867, ..., -4.29018529, -4.29870933, -4.30718067], [ -1.18452133, -1.1745215 , -1.16452267, ..., -1.70312928, -1.71165333, -1.72012467]], [[-22.70332918, -22.69332935, -22.68333051, ..., -23.22193713, -23.23046118, -23.23893251], [-20.95386052, -20.94386069, -20.93386185, ..., -21.47246847, -21.48099252, -21.48946385], [-18.99502781, -18.98502798, -18.97502915, ..., -19.51363576, -19.52215981, -19.53063115], ..., [ -6.80717583, -6.797176 , -6.78717717, ..., -7.32578378, -7.33430783, -7.34277917], [ -4.06169867, -4.05169884, -4.0417 , ..., -4.58030662, -4.58883067, -4.597302 ], [ -1.27563836, -1.26563852, -1.25563969, ..., -1.79424631, -1.80277036, -1.81124169]]])
- ns :
- None
- sample_rate :
- 1
- text :
- SEGY-SAK Create 3D Dataset
- measurement_system :
- None
- d3_domain :
- TWT
- epsg :
- None
- corner_points :
- [(1, 10), (19, 10), (19, 14), (1, 14), (1, 10)]
- corner_points_xy :
- None
- source_file :
- None
- srd :
- None
- datatype :
- None
- percentiles :
- None
- coord_scalar :
- None