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.
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.
segy_file = pathlib.Path("data/3D_gathers_pstm_nmo.sgy")
with pd.option_context("display.max_rows", 100):
display(segy_header_scan(segy_file))
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 |
UnassignedInt1 | 233 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
UnassignedInt2 | 237 | 1000.0 | 0.000000e+00 | 0.000000 | 0.0 | 0.00 | 0.0 | 0.00 | 0.0 |
penobscot_3d_gath = xr.open_dataset(
segy_file, dim_byte_fields={"iline":189, "xline":193, "offset":37}, extra_byte_fields={"cdp_x":181, "cdp_y":185}
)
# coordinates are not scaled automatically, to use the coordinate scalar found in the trace headers, we run scale_coords()
penobscot_3d_gath.segysak.scale_coords()
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.
display(penobscot_3d_gath)
print(penobscot_3d_gath.offset.values)
<xarray.Dataset> Size: 22MB Dimensions: (iline: 11, xline: 11, offset: 61, samples: 751) Coordinates: * iline (iline) int16 22B 1290 1291 1292 1293 1294 ... 1297 1298 1299 1300 * xline (xline) int16 22B 1150 1151 1152 1153 1154 ... 1157 1158 1159 1160 * offset (offset) int16 122B 175 225 275 325 375 ... 3025 3075 3125 3175 * samples (samples) float32 3kB 0.0 4.0 8.0 ... 2.992e+03 2.996e+03 3e+03 Data variables: cdp_x (iline, xline, offset) float64 59kB 7.336e+05 ... 7.338e+05 cdp_y (iline, xline, offset) float64 59kB 4.895e+06 ... 4.895e+06 data (iline, xline, offset, samples) float32 22MB ... Attributes: seisnc: {"coord_scalar": -10.0, "coord_scaled": true}
[ 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.
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
)
<matplotlib.collections.QuadMesh at 0x7fae250baa40>
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.
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",
)
<matplotlib.image.AxesImage at 0x7fae248985e0>
One can easily create a common offset stack by reversing the stacked dimension arguments "offset"
and "xline"
.
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",
)
<matplotlib.image.AxesImage at 0x7fae25c71f00>
Arbitrary line extraction on Gathers¶
Arbitrary line slicing of gathers based upon coordinates is also possible. Lets create a line that crosses the 3D.
penobscot_3d_gath
<xarray.Dataset> Size: 22MB Dimensions: (iline: 11, xline: 11, offset: 61, samples: 751) Coordinates: * iline (iline) int16 22B 1290 1291 1292 1293 1294 ... 1297 1298 1299 1300 * xline (xline) int16 22B 1150 1151 1152 1153 1154 ... 1157 1158 1159 1160 * offset (offset) int16 122B 175 225 275 325 375 ... 3025 3075 3125 3175 * samples (samples) float32 3kB 0.0 4.0 8.0 ... 2.992e+03 2.996e+03 3e+03 Data variables: cdp_x (iline, xline, offset) float64 59kB 7.336e+05 ... 7.338e+05 cdp_y (iline, xline, offset) float64 59kB 4.895e+06 ... 4.895e+06 data (iline, xline, offset, samples) float32 22MB ... Attributes: seisnc: {"coord_scalar": -10.0, "coord_scaled": true}
penobscot_3d_gath.drop_dims('offset')
<xarray.Dataset> Size: 3kB Dimensions: (iline: 11, xline: 11, samples: 751) Coordinates: * iline (iline) int16 22B 1290 1291 1292 1293 1294 ... 1297 1298 1299 1300 * xline (xline) int16 22B 1150 1151 1152 1153 1154 ... 1157 1158 1159 1160 * samples (samples) float32 3kB 0.0 4.0 8.0 ... 2.992e+03 2.996e+03 3e+03 Data variables: *empty* Attributes: seisnc: {"coord_scalar": -10.0, "coord_scaled": true}
penobscot_3d_gath.isel(offset=0).segysak.calc_corner_points()
[(733625.6000000001, 4895058.0), (733847.2000000001, 4895173.9), (733789.3, 4895284.600000001), (733567.7000000001, 4895168.8), (733625.6000000001, 4895058.0)]
arb_line = np.array([(733600, 4895180.0), (733850, 4895180.0)])
# we must drop the offset axis somehow, xy locations were loaded for each header position but
# they must be reduced, we could take the mean over the dimension but here we just select the
# first value using `offset=0`
ax = penobscot_3d_gath.isel(offset=0).segysak.plot_bounds()
ax.plot(arb_line[:, 0], arb_line[:, 1], label="arb_line")
plt.legend()
<matplotlib.legend.Legend at 0x7fae25008eb0>
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"
.
penobscot_3d_gath_arb = penobscot_3d_gath.segysak.xysel(
arb_line, method="nearest"
)
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",
)
<matplotlib.image.AxesImage at 0x7fae2482be80>
Muting and Stacking Gathers¶
Using one of our gathers let's define a mute function before we stack the data.
fig, axs = plt.subplots(ncols=2, figsize=(20, 10))
penobscot_3d_gath.isel(iline=5, xline=0).transpose('samples', ...).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.samples > mute)
# muted
penobscot_3d_gath_muted.isel(iline=5, xline=0).data.transpose('samples', ...).plot(
yincrease=False, ax=axs[1], vmax=5000
)
<matplotlib.collections.QuadMesh at 0x7fae24e34d00>
Stacking is the process of averaging the gathers for constant time to create a single trace per inline and crossline location.
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
)
<matplotlib.image.AxesImage at 0x7fae24d35d50>