Note
Get this example as a Jupyter Notebook ipynb
Working with 3D Gathers
Gathers or pre-stack data are common in seismic interpretation and processing. In this example we will load gather data by finding and specifying the offset byte location. Learn different ways to plot and select offset data. As well as perform a simple mute and trace stack to reduce the offset dimension.
[1]:
import pathlib
from IPython.display import display
import pandas as pd
import xarray as xr
import numpy as np
from segysak.segy import segy_loader, segy_header_scan
import matplotlib.pyplot as plt
This example uses a subset of the Penobscot 3D with data exported from the OpendTect project.
First we scan the data to determine which byte locations contain the relevant information. We will need to provide a byte location for the offset variable so a 4th dimension can be created when loaded the data.
[2]:
segy_file = pathlib.Path("data/3D_gathers_pstm_nmo.sgy")
with pd.option_context("display.max_rows", 100):
display(segy_header_scan(segy_file))
100%|██████████| 1.00k/1.00k [00:00<00:00, 11.1k traces/s]
| byte_loc | count | mean | std | min | 25% | 50% | 75% | max | |
|---|---|---|---|---|---|---|---|---|---|
| TRACE_SEQUENCE_LINE | 1 | 1000.0 | 2.797410e+02 | 186.100369 | 1.0 | 125.75 | 250.5 | 421.25 | 671.0 |
| TRACE_SEQUENCE_FILE | 5 | 1000.0 | 3.055600e+01 | 17.664623 | 1.0 | 15.00 | 30.0 | 46.00 | 61.0 |
| FieldRecord | 9 | 1000.0 | 1.290329e+03 | 0.470085 | 1290.0 | 1290.00 | 1290.0 | 1291.00 | 1291.0 |
| TraceNumber | 13 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| EnergySourcePoint | 17 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| CDP | 21 | 1000.0 | 1.154085e+03 | 3.039245 | 1150.0 | 1152.00 | 1154.0 | 1156.00 | 1160.0 |
| CDP_TRACE | 25 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| TraceIdentificationCode | 29 | 1000.0 | 1.000000e+00 | 0.000000 | 1.0 | 1.00 | 1.0 | 1.00 | 1.0 |
| NSummedTraces | 31 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| NStackedTraces | 33 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| DataUse | 35 | 1000.0 | 1.000000e+00 | 0.000000 | 1.0 | 1.00 | 1.0 | 1.00 | 1.0 |
| offset | 37 | 1000.0 | 1.652800e+03 | 883.231146 | 175.0 | 875.00 | 1625.0 | 2425.00 | 3175.0 |
| ReceiverGroupElevation | 41 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| SourceSurfaceElevation | 45 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| SourceDepth | 49 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| ReceiverDatumElevation | 53 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| SourceDatumElevation | 57 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| SourceWaterDepth | 61 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| GroupWaterDepth | 65 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| ElevationScalar | 69 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| SourceGroupScalar | 71 | 1000.0 | -1.000000e+01 | 0.000000 | -10.0 | -10.00 | -10.0 | -10.00 | -10.0 |
| SourceX | 73 | 1000.0 | 7.337142e+06 | 685.656940 | 7336198.0 | 7336642.00 | 7337085.0 | 7337586.00 | 7338472.0 |
| SourceY | 77 | 1000.0 | 4.895109e+07 | 333.157506 | 48950580.0 | 48950812.00 | 48951044.0 | 48951276.00 | 48951739.0 |
| GroupX | 81 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| GroupY | 85 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| CoordinateUnits | 89 | 1000.0 | 1.000000e+00 | 0.000000 | 1.0 | 1.00 | 1.0 | 1.00 | 1.0 |
| WeatheringVelocity | 91 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| SubWeatheringVelocity | 93 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| SourceUpholeTime | 95 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| GroupUpholeTime | 97 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| SourceStaticCorrection | 99 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| GroupStaticCorrection | 101 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| TotalStaticApplied | 103 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| LagTimeA | 105 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| LagTimeB | 107 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| DelayRecordingTime | 109 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| MuteTimeStart | 111 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| MuteTimeEND | 113 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| TRACE_SAMPLE_COUNT | 115 | 1000.0 | 7.510000e+02 | 0.000000 | 751.0 | 751.00 | 751.0 | 751.00 | 751.0 |
| TRACE_SAMPLE_INTERVAL | 117 | 1000.0 | 4.000000e+03 | 0.000000 | 4000.0 | 4000.00 | 4000.0 | 4000.00 | 4000.0 |
| GainType | 119 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| InstrumentGainConstant | 121 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| InstrumentInitialGain | 123 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| Correlated | 125 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| SweepFrequencyStart | 127 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| SweepFrequencyEnd | 129 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| SweepLength | 131 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| SweepType | 133 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| SweepTraceTaperLengthStart | 135 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| SweepTraceTaperLengthEnd | 137 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| TaperType | 139 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| AliasFilterFrequency | 141 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| AliasFilterSlope | 143 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| NotchFilterFrequency | 145 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| NotchFilterSlope | 147 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| LowCutFrequency | 149 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| HighCutFrequency | 151 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| LowCutSlope | 153 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| HighCutSlope | 155 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| YearDataRecorded | 157 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| DayOfYear | 159 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| HourOfDay | 161 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| MinuteOfHour | 163 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| SecondOfMinute | 165 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| TimeBaseCode | 167 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| TraceWeightingFactor | 169 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| GeophoneGroupNumberRoll1 | 171 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| GeophoneGroupNumberFirstTraceOrigField | 173 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| GeophoneGroupNumberLastTraceOrigField | 175 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| GapSize | 177 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| OverTravel | 179 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| CDP_X | 181 | 1000.0 | 7.337142e+06 | 685.656940 | 7336198.0 | 7336642.00 | 7337085.0 | 7337586.00 | 7338472.0 |
| CDP_Y | 185 | 1000.0 | 4.895109e+07 | 333.157506 | 48950580.0 | 48950812.00 | 48951044.0 | 48951276.00 | 48951739.0 |
| INLINE_3D | 189 | 1000.0 | 1.290329e+03 | 0.470085 | 1290.0 | 1290.00 | 1290.0 | 1291.00 | 1291.0 |
| CROSSLINE_3D | 193 | 1000.0 | 1.154085e+03 | 3.039245 | 1150.0 | 1152.00 | 1154.0 | 1156.00 | 1160.0 |
| ShotPoint | 197 | 1000.0 | 3.055600e+01 | 17.664623 | 1.0 | 15.00 | 30.0 | 46.00 | 61.0 |
| ShotPointScalar | 201 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| TraceValueMeasurementUnit | 203 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| TransductionConstantMantissa | 205 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| TransductionConstantPower | 209 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| TransductionUnit | 211 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| TraceIdentifier | 213 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| ScalarTraceHeader | 215 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| SourceType | 217 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| SourceEnergyDirectionMantissa | 219 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| SourceEnergyDirectionExponent | 223 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| SourceMeasurementMantissa | 225 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| SourceMeasurementExponent | 229 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
| SourceMeasurementUnit | 231 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
[3]:
penobscot_3d_gath = segy_loader(
segy_file, iline=189, xline=193, cdpx=181, cdpy=185, offset=37
)
100%|██████████| 7.38k/7.38k [00:00<00:00, 11.3k traces/s]
Loading as 3D
Fast direction is INLINE_3D
Converting SEGY: 100%|██████████| 7.38k/7.38k [00:02<00:00, 3.63k traces/s]
Note that the loaded Dataset has four dimensions with the additional dimension labeled offset. There are 61 offsets in this dataset or 61 traces per inline and xline location.
[4]:
display(penobscot_3d_gath)
print(penobscot_3d_gath.offset.values)
<xarray.Dataset>
Dimensions: (iline: 11, xline: 11, twt: 751, offset: 61)
Coordinates:
* iline (iline) uint16 1290 1291 1292 1293 1294 ... 1297 1298 1299 1300
* xline (xline) uint16 1150 1151 1152 1153 1154 ... 1157 1158 1159 1160
* twt (twt) float64 0.0 4.0 8.0 12.0 ... 2.992e+03 2.996e+03 3e+03
* offset (offset) uint16 175 225 275 325 375 ... 2975 3025 3075 3125 3175
cdp_x (iline, xline, offset) float32 7.336e+05 7.336e+05 ... 7.338e+05
cdp_y (iline, xline, offset) float32 4.895e+06 4.895e+06 ... 4.895e+06
Data variables:
data (iline, xline, twt, offset) float32 0.0 0.0 0.0 ... 1.104e+03 670.1
Attributes: (12/13)
ns: None
sample_rate: 4.0
text: C01\nC02\nC03\nC04\nC05\nC06\nC07 Positions accordi...
measurement_system: m
d3_domain: None
epsg: None
... ...
corner_points_xy: None
source_file: 3D_gathers_pstm_nmo.sgy
srd: None
datatype: None
percentiles: [-4084.319948142755, -3961.308399237772, -1082.87965...
coord_scalar: -10.0
[ 175 225 275 325 375 425 475 525 575 625 675 725 775 825
875 925 975 1025 1075 1125 1175 1225 1275 1325 1375 1425 1475 1525
1575 1625 1675 1725 1775 1825 1875 1925 1975 2025 2075 2125 2175 2225
2275 2325 2375 2425 2475 2525 2575 2625 2675 2725 2775 2825 2875 2925
2975 3025 3075 3125 3175]
Lets check that the data looks OK for a couple of offsets. We’ve only got a small dataset of 11x11 traces so the seismic will look at little odd at this scale.
[5]:
fig, axs = plt.subplots(ncols=2, figsize=(20, 10))
penobscot_3d_gath.isel(iline=5, offset=0).data.T.plot(
yincrease=False, ax=axs[0], vmax=5000
)
penobscot_3d_gath.isel(xline=5, offset=0).data.T.plot(
yincrease=False, ax=axs[1], vmax=5000
)
[5]:
<matplotlib.collections.QuadMesh at 0x7fc4fc13bdc0>
Plotting Gathers Sequentially
Plotting of gathers is often done in a stacked way, displaying sequential gathers along a common dimension, usually inline or crossline. Xarray provides the stack method which can be used to stack labelled dimensions together.
[6]:
fig, axs = plt.subplots(figsize=(20, 10))
axs.imshow(
penobscot_3d_gath.isel(iline=0)
.data.stack(stacked_offset=("xline", "offset"))
.values,
vmin=-5000,
vmax=5000,
cmap="seismic",
aspect="auto",
)
[6]:
<matplotlib.image.AxesImage at 0x7fc4fac6da00>
One can easily create a common offset stack by reversing the stacked dimension arguments "offset" and "xline".
[7]:
fig, axs = plt.subplots(figsize=(20, 10))
axs.imshow(
penobscot_3d_gath.isel(iline=0)
.data.stack(stacked_offset=("offset", "xline"))
.values,
vmin=-5000,
vmax=5000,
cmap="seismic",
aspect="auto",
)
[7]:
<matplotlib.image.AxesImage at 0x7fc4fabf6100>
Arbitrary line extraction on Gathers
Arbitrary line slicing of gathers based upon coordinates is also possible. Lets create a line that crosses the 3D.
[8]:
arb_line = np.array([(733600, 733850), (4895180.0, 4895180.0)])
ax = penobscot_3d_gath.seis.plot_bounds()
ax.plot(arb_line[0, :], arb_line[1, :], label="arb_line")
plt.legend()
[8]:
<matplotlib.legend.Legend at 0x7fc4f8b05310>
Here we need to think carefully about the bin_spacing_hint. We also don’t want to interpolate the gathers, so we use xysel_method="nearest".
[9]:
penobscot_3d_gath_arb = penobscot_3d_gath.seis.interp_line(
arb_line[0, :], arb_line[1, :], bin_spacing_hint=30, xysel_method="nearest"
)
[10]:
fig, axs = plt.subplots(figsize=(20, 10))
axs.imshow(
penobscot_3d_gath_arb.data.stack(
stacked_offset=(
"cdp",
"offset",
)
).values,
vmin=-5000,
vmax=5000,
cmap="seismic",
aspect="auto",
)
[10]:
<matplotlib.image.AxesImage at 0x7fc4f89ef130>
Muting and Stacking Gathers
Using one of our gathers let’s define a mute function before we stack the data.
[11]:
fig, axs = plt.subplots(ncols=2, figsize=(20, 10))
penobscot_3d_gath.isel(iline=5, xline=0).data.plot(
yincrease=False, ax=axs[0], vmax=5000
)
# the mute relates the offset to an expected twt, let's just use a linear mute for this example
mute = penobscot_3d_gath.offset * 0.6 + 300
# and then we can plot it up
mute.plot(ax=axs[0], color="k")
# apply the mute to the volume
penobscot_3d_gath_muted = penobscot_3d_gath.where(penobscot_3d_gath.twt > mute)
# muted
penobscot_3d_gath_muted.isel(iline=5, xline=0).data.plot(
yincrease=False, ax=axs[1], vmax=5000
)
[11]:
<matplotlib.collections.QuadMesh at 0x7fc4f88af040>
Stacking is the process of averaging the gathers for constant time to create a single trace per inline and crossline location.
[12]:
fig, axs = plt.subplots(ncols=3, figsize=(20, 10))
plot_kwargs = dict(vmax=5000, interpolation="bicubic", yincrease=False)
# compare with the zero offset trace and use imshow for interpolation
penobscot_3d_gath.isel(iline=5, offset=0).data.T.plot.imshow(ax=axs[0], **plot_kwargs)
# stack the no mute data
penobscot_3d_gath.isel(iline=5).data.mean("offset").T.plot.imshow(
ax=axs[1], **plot_kwargs
)
# stack the muted data
penobscot_3d_gath_muted.isel(iline=5).data.mean("offset").T.plot.imshow(
ax=axs[2], **plot_kwargs
)
[12]:
<matplotlib.image.AxesImage at 0x7fc4f86a9520>
