Get Started in 10 Minutes

Altay Sansal

Sep 22, 2026

8 min read

In this page we will be showing basic capabilities of MDIO.

For demonstration purposes, we will ingest a public raw shot record from the Soda Lake 2010 3D/3C seismic survey (Geothermal Data Repository / AWS Open Data, CC-BY-4.0).

We use one ~7.5 MiB field record: F7733R1.SGY. No AWS credentials are required. The file is a 3C common-source gather: shot 7733, 958 traces (1 aux + 3 × 319 receivers), 2000 samples at 2 ms.

Warning

For plotting and remote ingestion the notebook requires Matplotlib and aiohttp as a dependency. Please install it before executing using pip install matplotlib aiohttp or conda install matplotlib aiohttp.

Defining the SEG-Y Dataset

Since MDIO 0.8 we can directly ingest remote SEG-Y files. To make the header scan performant we can also set up an environment variable for MDIO. See here for when to use this: Buffered Reads in Ingestion.

The gather is a single live shot, so the live mask is all True. The first trace is an auxiliary channel at (0, 0) and is much louder than the seismic. Plots below drop it and keep one of the three components.

The following environment variables are essential here:

  • MDIO__IMPORT__CLOUD_NATIVE tells MDIO to do buffered reads for headers due to remote file.

  • MDIO__IMPORT__SAVE_SEGY_FILE_HEADER makes MDIO save the SEG-Y specific file headers (text, binary) which is not strictly necessary for consumption and is disabled by default.

import os

os.environ["MDIO__IMPORT__CLOUD_NATIVE"] = "true"
os.environ["MDIO__IMPORT__SAVE_SEGY_FILE_HEADER"] = "true"

input_url = "https://gdr-data-lake.s3.us-west-2.amazonaws.com/soda_lake/raw_seismic/2010/v1.0.0/F7733R1.SGY"

Ingesting to MDIO

To do this, we can use the convenient SEG-Y to MDIO converter.

The shot identifier is the field record number at byte 9, and the channel number is at byte 13.

MDIO uses TGSAI/segy to parse the SEG-Y; the field names conform to its canonical keys defined in SEGY Binary Header Keys and SEGY Trace Header Keys. Since MDIO v1 we also introduced templates for common seismic data types. For this shot gather we use the StreamerShotGathers2D template here, which expects shot_point and channel.

We will also specify the units for the time domain. The spatial units will be automatically parsed from SEG-Y binary header. However, there may be a case where it is corrupt in the file, for that see the Fixing X/Y Units Issues section.

In summary, we will use the byte locations as defined for ingestion.

import matplotlib.pyplot as plt
from segy.schema import HeaderField
from segy.standards import get_segy_standard

from mdio import segy_to_mdio
from mdio.builder.schemas.v1.units import TimeUnitModel
from mdio.builder.template_registry import get_template

shot_trace_headers = [
    HeaderField(name="shot_point", byte=9, format="int32"),
    HeaderField(name="channel", byte=13, format="int32"),
]

rev1_segy_spec = get_segy_standard(1.0)
soda_lake_segy_spec = rev1_segy_spec.customize(trace_header_fields=shot_trace_headers)

mdio_template = get_template("StreamerShotGathers2D")
unit_ms = TimeUnitModel(time="ms")
mdio_template.add_units({"time": unit_ms})

segy_to_mdio(
    input_path=input_url,
    output_path="soda_lake.mdio",
    segy_spec=soda_lake_segy_spec,
    mdio_template=mdio_template,
    overwrite=True,
)

It only took a few seconds to ingest, since this is a very small file.

However, MDIO scales up to TB (that’s ~1,000 GB) sized volumes!

Opening the Ingested MDIO File

Let’s open the MDIO file with the open_mdio function. This will return a pretty xarray representation with our standardized format.

from mdio import open_mdio

dataset = open_mdio("soda_lake.mdio")
dataset
<xarray.Dataset> Size: 8MB
Dimensions:           (shot_point: 1, channel: 958, time: 2000)
Coordinates:
  * shot_point        (shot_point) int32 4B 7733
    source_coord_x    (shot_point) float64 8B ...
    source_coord_y    (shot_point) float64 8B ...
  * channel           (channel) int32 4kB 1 2 3 4 5 6 ... 954 955 956 957 958
    group_coord_y     (shot_point, channel) float64 8kB ...
    group_coord_x     (shot_point, channel) float64 8kB ...
  * time              (time) int32 8kB 0 2 4 6 8 10 ... 3990 3992 3994 3996 3998
Data variables:
    segy_file_header  <U1 4B ...
    trace_mask        (shot_point, channel) bool 958B ...
    headers           (shot_point, channel) [('trace_seq_num_line', '<i4'), ('trace_seq_num_reel', '<i4'), ('energy_source_point_num', '<i4'), ('ensemble_num', '<i4'), ('trace_num_ensemble', '<i4'), ('trace_id_code', '<i2'), ('vertically_summed_traces', '<i2'), ('horizontally_stacked_traces', '<i2'), ('data_use', '<i2'), ('source_to_receiver_distance', '<i4'), ('receiver_group_elevation', '<i4'), ('source_surface_elevation', '<i4'), ('source_depth_below_surface', '<i4'), ('receiver_datum_elevation', '<i4'), ('source_datum_elevation', '<i4'), ('source_water_depth', '<i4'), ('receiver_water_depth', '<i4'), ('elevation_depth_scalar', '<i2'), ('coordinate_scalar', '<i2'), ('source_coord_x', '<i4'), ('source_coord_y', '<i4'), ('group_coord_x', '<i4'), ('group_coord_y', '<i4'), ('coordinate_unit', '<i2'), ('weathering_velocity', '<i2'), ('subweathering_velocity', '<i2'), ('source_uphole_time', '<i2'), ('group_uphole_time', '<i2'), ('source_static_correction', '<i2'), ('receiver_static_correction', '<i2'), ('total_static_applied', '<i2'), ('lag_time_a', '<i2'), ('lag_time_b', '<i2'), ('delay_recording_time', '<i2'), ('mute_time_start', '<i2'), ('mute_time_end', '<i2'), ('samples_per_trace', '<i2'), ('sample_interval', '<i2'), ('instrument_gain_type', '<i2'), ('instrument_gain_const', '<i2'), ('instrument_gain_initial', '<i2'), ('correlated_data', '<i2'), ('sweep_freq_start', '<i2'), ('sweep_freq_end', '<i2'), ('sweep_length', '<i2'), ('sweep_type', '<i2'), ('sweep_taper_start', '<i2'), ('sweep_taper_end', '<i2'), ('taper_type', '<i2'), ('alias_filter_freq', '<i2'), ('alias_filter_slope', '<i2'), ('notch_filter_freq', '<i2'), ('notch_filter_slope', '<i2'), ('low_cut_freq', '<i2'), ('high_cut_freq', '<i2'), ('low_cut_slope', '<i2'), ('high_cut_slope', '<i2'), ('year_recorded', '<i2'), ('day_of_year', '<i2'), ('hour_of_day', '<i2'), ('minute_of_hour', '<i2'), ('second_of_minute', '<i2'), ('time_basis_code', '<i2'), ('trace_weighting_factor', '<i2'), ('group_num_roll_switch', '<i2'), ('group_num_first_trace', '<i2'), ('group_num_last_trace', '<i2'), ('gap_size', '<i2'), ('taper_overtravel', '<i2'), ('cdp_x', '<i4'), ('cdp_y', '<i4'), ('inline', '<i4'), ('crossline', '<i4'), ('shot_point', '<i4'), ('shot_point_scalar', '<i2'), ('trace_value_unit', '<i2'), ('transduction_const_mantissa', '<i4'), ('transduction_const_exponent', '<i2'), ('transduction_unit', '<i2'), ('device_trace_id', '<i2'), ('times_scalar', '<i2'), ('source_type_orientation', '<i2'), ('source_energy_dir_mantissa', '<i4'), ('source_energy_dir_exponent', '<i2'), ('source_measurement_mantissa', '<i4'), ('source_measurement_exponent', '<i2'), ('source_measurement_unit', '<i2'), ('channel', '<i4')] 218kB ...
    amplitude         (shot_point, channel, time) float32 8MB ...
Attributes:
    apiVersion:  1.2.3
    createdOn:   2026-09-22 21:03:51.528919+00:00
    name:        StreamerShotGathers2D
    attributes:  {'surveyType': '2D', 'gatherType': 'common_source', 'default...

Querying Metadata

Now let’s look at the file text header saved in the segy_file_header metadata variable.

print(dataset["segy_file_header"].attrs["textHeader"])
C 1 CLIENT                        COMPANY                       CREW NO         
C 2 LINE            AREA                        MAP ID                          
C 3 REEL NO           DAY-START OF REEL     YEAR      OBSERVER                  
C 4 INSTRUMENT: MFG            MODEL            SERIAL NO                       
C 5 DATA TRACES/RECORD        AUXILIARY TRACES/RECORD         CDP FOLD          
C 6 SAMPLE INTERVAL         SAMPLES/TRACE       BITS/IN      BYTES/SAMPLE       
C 7 RECORDING FORMAT        FORMAT THIS REEL        MEASUREMENT SYSTEM          
C 8 SAMPLE CODE: FLOATING PT     FIXED PT     FIXED PT-GAIN     CORRELATED      
C 9 GAIN  TYPE: FIXED     BINARY     FLOATING POINT     OTHER                   
C10 FILTERS: ALIAS     HZ  NOTCH     HZ  BAND     -     HZ  SLOPE    -    DB/OCT
C11 SOURCE: TYPE            NUMBER/POINT        POINT INTERVAL                  
C12     PATTERN:                           LENGTH        WIDTH                  
C13 SWEEP: START     HZ  END     HZ  LENGTH      MS  CHANNEL NO     TYPE        
C14 TAPER: START LENGTH       MS  END LENGTH       MS  TYPE                     
C15 SPREAD: OFFSET        MAX DISTANCE        GROUP INTERVAL                    
C16 GEOPHONES: PER GROUP     SPACING     FREQUENCY     MFG          MODEL       
C17     PATTERN:                           LENGTH        WIDTH                  
C18 TRACES SORTED BY: RECORD     CDP     OTHER                                  
C19 AMPLITUDE RECOVERY: NONE      SPHERICAL DIV      AGC     OTHER              
C20 MAP PROJECTION                      ZONE ID       COORDINATE UNITS          
C21 PROCESSING:                                                                 
C22 PROCESSING:                                                                 
C23                                                                             
C24                                                                             
C25                                                                             
C26                                                                             
C27                                                                             
C28                                                                             
C29                                                                             
C30                                                                             
C31                                                                             
C32                                                                             
C33                                                                             
C34                                                                             
C35                                                                             
C36                                                                             
C37                                                                             
C38                                                                             
C39                                                                             
C40 END EBCDIC                                                                  

Since we saved the binary header, we can look at that as well.

dataset["segy_file_header"].attrs["binaryHeader"]
{'job_id': 18909,
 'line_num': 1,
 'reel_num': 1,
 'data_traces_per_ensemble': 957,
 'aux_traces_per_ensemble': 1,
 'sample_interval': 2000,
 'orig_sample_interval': 2000,
 'samples_per_trace': 2000,
 'orig_samples_per_trace': 2000,
 'data_sample_format': 1,
 'ensemble_fold': 0,
 'trace_sorting_code': 1,
 'vertical_sum_code': 1,
 'sweep_freq_start': 0,
 'sweep_freq_end': 0,
 'sweep_length': 0,
 'sweep_type_code': 4,
 'sweep_trace_num': 0,
 'sweep_taper_start': 0,
 'sweep_taper_end': 0,
 'taper_type_code': 3,
 'correlated_data_code': 1,
 'binary_gain_code': 2,
 'amp_recovery_code': 1,
 'measurement_system_code': 2,
 'impulse_polarity_code': 0,
 'vibratory_polarity_code': 0,
 'fixed_length_trace_flag': 1,
 'num_extended_text_headers': 0,
 'segy_revision_major': 0,
 'segy_revision_minor': 0}

Fetching Data and Plotting

We keep one component (transduction_const_mantissa == 2), clip by the 98th percentile, and plot against source-receiver offset. Channel order follows the 2D receiver patch, not offset.

Since MDIO v1 we are using Xarray under the hood, so we can use its convenient indexing. It also handles the plotting with proper dimension coordinate labels.

import numpy as np

COMPONENT_A = 2

shot = dataset.sel(shot_point=7733)
mantissa = shot.headers.to_numpy()["transduction_const_mantissa"]
comp = shot.isel(channel=mantissa == COMPONENT_A)

offset = np.hypot(
    comp.group_coord_x.to_numpy() - float(comp.source_coord_x),
    comp.group_coord_y.to_numpy() - float(comp.source_coord_y),
)
comp = comp.assign_coords(offset=("channel", offset)).sortby("offset")

amp = comp["amplitude"]
clip = float(np.percentile(np.abs(amp.to_numpy()), 98))

amp.plot(x="offset", y="time", vmin=-clip, vmax=clip, cmap="gray_r", yincrease=False);
../_images/4b9ca4ae15a67eb520e133a1029a910a34177e96f97678e0d0c5a98df17a6c87.png

Receiver coordinates in map view, then the same points colored by amplitude at 1,000 ms.

fig, ax = plt.subplots(figsize=(6, 5))
ax.scatter(comp.group_coord_x, comp.group_coord_y, s=16, c="0.3", label="receivers")
ax.scatter(
    float(comp.source_coord_x),
    float(comp.source_coord_y),
    marker="*",
    s=180,
    c="tab:red",
    label="source",
    zorder=3,
)
ax.set_xlabel("group_coord_x")
ax.set_ylabel("group_coord_y")
ax.set_aspect("equal", adjustable="box")
ax.legend()
ax.set_title("component A receiver patch");
../_images/7d70530e547cfe03ec6d2b2b8c0968cf251b37fee36a03fcfed4479c888ed8d3.png

Color the same map by amplitude at 1,000 ms.

twt = amp.sel(time=1000)

fig, ax = plt.subplots(figsize=(6, 5))
sc = ax.scatter(
    comp.group_coord_x,
    comp.group_coord_y,
    c=twt,
    s=16,
    cmap="gray_r",
    vmin=-clip,
    vmax=clip,
)
ax.scatter(
    float(comp.source_coord_x),
    float(comp.source_coord_y),
    marker="*",
    s=180,
    c="tab:red",
    label="source",
    zorder=3,
)
ax.set_xlabel("group_coord_x")
ax.set_ylabel("group_coord_y")
ax.set_aspect("equal", adjustable="box")
ax.legend()
plt.colorbar(sc, ax=ax, label="amplitude @ 1000 ms");
../_images/7233345b3c25194323a4f56995615ee084727d517699bf78891280ce650c5eb9.png

Query Headers

We can query headers for the whole dataset very quickly because they are separated from the seismic wavefield.

Let’s get all the headers for source and receiver coordinates in this dataset.

Note that the header maps will still share the geometry/grid of the dataset!

The compute property fetches the lazily opened MDIO values.

dataset.headers["group_coord_x"].compute()
<xarray.DataArray 'group_coord_x' (shot_point: 1, channel: 958)> Size: 8kB
array([[      0., 7734489., 7735044., 7735667., 7736158., 7736714.,
        7737268., 7737825., 7738382., 7738938., 7739494., 7741720.,
        7742273., 7742831., 7743326., 7743945., 7744504., 7745053.,
        7745611., 7746169., 7746725., 7735425., 7735981., 7736537.,
        7737093., 7737650., 7738204., 7738762., 7739319., 7739874.,
        7740432., 7740987., 7741544., 7742098., 7742657., 7743212.,
        7743767., 7744323., 7744879., 7745435., 7745947., 7746549.,
        7747660., 7736361., 7736909., 7737472., 7738028., 7738587.,
        7739142., 7739699., 7740255., 7740810., 7741366., 7741925.,
        7742480., 7743037., 7743591., 7744147., 7744704., 7745258.,
        7745818., 7746372., 7746875., 7737298., 7737856., 7738411.,
        7738969., 7739523., 7740086., 7740636., 7741192., 7741747.,
        7742321., 7742858., 7743416., 7743973., 7744529., 7745083.,
        7745638., 7746197., 7746717., 7748979., 7749532., 7738232.,
        7738789., 7739347., 7739904., 7740457., 7741018., 7741565.,
        7742129., 7742683., 7743236., 7743790., 7744350., 7744908.,
        7745463., 7746023., 7746577., 7747132., 7747690., 7748244.,
        7748800., 7749360., 7749914., 7750469., 7739172., 7739720.,
        7740280., 7740842., 7741404., 7741955., 7742505., 7743068.,
        7743616., 7744177., 7744732., 7745294., 7745847., 7746406.,
...
        7745143., 7745699., 7746254., 7746808., 7747367., 7747923.,
        7748480., 7749035., 7749591., 7750145., 7750704., 7751260.,
        7751817., 7752373., 7752928., 7753486., 7754042., 7754596.,
        7755154., 7743853., 7744409., 7744965., 7745521., 7746078.,
        7746633., 7747192., 7747747., 7748303., 7748859., 7749416.,
        7749972., 7750529., 7751083., 7751639., 7752199., 7752752.,
        7753307., 7753864., 7754421., 7754977., 7755534., 7756089.,
        7744789., 7745346., 7745903., 7746459., 7747015., 7747571.,
        7748128., 7748684., 7749237., 7749796., 7750351., 7750908.,
        7751464., 7752021., 7752577., 7753135., 7753688., 7754246.,
        7754802., 7755357., 7755914., 7756470., 7757026., 7745728.,
        7746283., 7746840., 7747394., 7747951., 7748506., 7750175.,
        7750731., 7751289., 7751845., 7752402., 7752958., 7753514.,
        7754070., 7754625., 7755170., 7755736., 7756295., 7756852.,
        7757405., 7757961., 7747221., 7747775., 7748332., 7748888.,
        7752779., 7753338., 7753893., 7754450., 7755005., 7755564.,
        7756119., 7756676., 7757229., 7757787., 7758343., 7758899.,
        7748711., 7749269., 7749826., 7753161., 7753717., 7754272.,
        7754831., 7755385., 7755940., 7756501., 7757053., 7757613.,
        7758169., 7758724., 7759277., 7759836.]])
Coordinates:
  * shot_point      (shot_point) int32 4B 7733
    source_coord_x  (shot_point) float64 8B 7.735e+06
    source_coord_y  (shot_point) float64 8B 4.534e+07
  * channel         (channel) int32 4kB 1 2 3 4 5 6 ... 953 954 955 956 957 958
    group_coord_y   (shot_point, channel) float64 8kB 0.0 ... 4.535e+07
    group_coord_x   (shot_point, channel) float64 8kB 0.0 7.734e+06 ... 7.76e+06
Attributes:
    unitsV1:  {'length': 'ft'}
dataset.headers["group_coord_y"].compute()
<xarray.DataArray 'group_coord_y' (shot_point: 1, channel: 958)> Size: 8kB
array([[       0., 45338731., 45338358., 45337937., 45337606., 45337232.,
        45336858., 45336483., 45336108., 45335735., 45335358., 45333861.,
        45333488., 45333112., 45332779., 45332367., 45332013., 45331615.,
        45331238., 45330866., 45330488., 45340121., 45339748., 45339376.,
        45338997., 45338622., 45338248., 45337872., 45337501., 45337122.,
        45336752., 45336440., 45336001., 45335627., 45335252., 45334878.,
        45334503., 45334128., 45333754., 45333380., 45332941., 45332630.,
        45331881., 45341515., 45341140., 45340762., 45340388., 45340016.,
        45339638., 45339266., 45338892., 45338516., 45338142., 45337767.,
        45337393., 45337019., 45336643., 45336269., 45335893., 45335518.,
        45335144., 45334769., 45334431., 45342905., 45342531., 45342145.,
        45341778., 45341405., 45341034., 45340657., 45340279., 45339906.,
        45339520., 45339156., 45338783., 45338406., 45338034., 45337660.,
        45337284., 45336907., 45336559., 45335038., 45334661., 45344294.,
        45343919., 45343547., 45343172., 45342797., 45342423., 45342048.,
        45341668., 45341298., 45340919., 45340557., 45340165., 45339796.,
        45339422., 45339048., 45338672., 45338298., 45337925., 45337549.,
        45337174., 45336802., 45336426., 45336052., 45345684., 45345307.,
        45344935., 45344564., 45344186., 45343811., 45343437., 45343065.,
        45342690., 45342315., 45341943., 45341561., 45341191., 45340814.,
...
        45349747., 45349373., 45348997., 45348623., 45348247., 45347873.,
        45347499., 45347124., 45346750., 45346375., 45346001., 45345628.,
        45345251., 45344878., 45344501., 45344127., 45343754., 45343379.,
        45343004., 45352638., 45352261., 45351886., 45351512., 45351138.,
        45350761., 45350389., 45350015., 45349640., 45349265., 45348890.,
        45348515., 45348142., 45347766., 45347392., 45347016., 45346642.,
        45346268., 45345892., 45345519., 45345145., 45344769., 45344393.,
        45354027., 45353652., 45353276., 45352902., 45352527., 45352152.,
        45351779., 45351403., 45351030., 45350656., 45350280., 45349907.,
        45349531., 45349158., 45348781., 45348409., 45348032., 45347657.,
        45347284., 45346909., 45346533., 45346158., 45345785., 45355416.,
        45355043., 45354668., 45354294., 45353918., 45353543., 45352417.,
        45352045., 45351670., 45351295., 45350923., 45350547., 45350171.,
        45349797., 45349424., 45349051., 45348675., 45348301., 45347925.,
        45347550., 45347174., 45356434., 45356059., 45355685., 45355310.,
        45352687., 45352312., 45351937., 45351564., 45351189., 45350814.,
        45350440., 45350064., 45349691., 45349315., 45348941., 45348565.,
        45357449., 45357074., 45356700., 45354450., 45354078., 45353704.,
        45353328., 45352954., 45352578., 45352202., 45351828., 45351458.,
        45351080., 45350706., 45350331., 45349955.]])
Coordinates:
  * shot_point      (shot_point) int32 4B 7733
    source_coord_x  (shot_point) float64 8B 7.735e+06
    source_coord_y  (shot_point) float64 8B 4.534e+07
  * channel         (channel) int32 4kB 1 2 3 4 5 6 ... 953 954 955 956 957 958
    group_coord_y   (shot_point, channel) float64 8kB 0.0 ... 4.535e+07
    group_coord_x   (shot_point, channel) float64 8kB 0.0 7.734e+06 ... 7.76e+06
Attributes:
    unitsV1:  {'length': 'ft'}

As we mentioned before, we can also get specific dataset slices of headers while fetching a slice.

Let’s overlay receiver coordinates on the offset-sorted gather.

group_x_header = comp["group_coord_x"]
group_y_header = comp["group_coord_y"]
image = comp["amplitude"]

gs_kw = {"height_ratios": (1, 8)}
fig, (hdr_ax, img_ax) = plt.subplots(2, 1, figsize=(7, 6), gridspec_kw=gs_kw, sharex="all")
hdr_ax2 = hdr_ax.twinx()

group_x_header.plot(ax=hdr_ax, x="offset", c="red")
group_y_header.plot(ax=hdr_ax2, x="offset", c="black")
image.plot(
    ax=img_ax,
    x="offset",
    y="time",
    vmin=-clip,
    vmax=clip,
    cmap="gray_r",
    yincrease=False,
    add_colorbar=False,
)
img_ax.set_title("")
hdr_ax.set_xlabel("")

plt.tight_layout()
../_images/9c30dbe1fe1e00f75dec72de7a17ac269e929ce2b493f6c98c0649ab8fc5f369.png

MDIO to SEG-Y Conversion

Finally, let’s demonstrate going back to SEG-Y.

We will use the convenient mdio_to_segy function and write it out as a round-trip file.

The output spec can be modified if we want to write things to different byte locations, etc, but we will use the same one as before.

from mdio import mdio_to_segy

mdio_to_segy(
    input_path="soda_lake.mdio",
    output_path="soda_lake_roundtrip.sgy",
    segy_spec=soda_lake_segy_spec,
)

Validate Round-Trip SEG-Y File

We can validate if the round-trip SEG-Y file is matching the original using TGSAI/segy.

Step by step:

  • Open original file

  • Open round-trip file

  • Compare text headers

  • Compare binary headers

  • Compare 100 random headers and traces

import numpy as np
from segy import SegyFile

original_segy = SegyFile(input_url)
roundtrip_segy = SegyFile("soda_lake_roundtrip.sgy")

# Compare text header
assert original_segy.text_header == roundtrip_segy.text_header

# Compare bin header
assert original_segy.binary_header == roundtrip_segy.binary_header

# Compare 100 random trace headers and samples.
# IBM float32 cannot preserve denormals through IEEE decode/encode.
rng = np.random.default_rng()
rand_indices = rng.integers(low=0, high=original_segy.num_traces, size=100)
for idx in rand_indices:
    np.testing.assert_equal(original_segy.header[idx], roundtrip_segy.header[idx])
    np.testing.assert_allclose(original_segy.sample[idx], roundtrip_segy.sample[idx], rtol=0, atol=1e-37)

print("Files identical!")
Files identical!