Read velocity model from segy file and get domain extents.

[krober10nd]

Hello,

I'd like to use segyio to read in a velocity model (https://wiki.seg.org/wiki/2004_BP_velocity_estimation_benchmark_model) and besides getting the velocity model, also get the domain extents.

I tried this:

 def read_segy(fname):
     """Read a velocity model from a SEG-y file"""
     with segyio.open(fname, ignore_geometry=True) as f:
         nz, nx = len(f.samples), len(f.trace)
         vp = np.zeros(shape=(nz, nx))
         for index, trace in enumerate(f.trace):
             vp[:, index] = trace
         return np.flipud(vp) 

But with this, I only can get the nz and nx but I can't determine the grid spacing information which would lead me to figuring out the domain extents. I can obviously get this information about the grid spacings in x and y from the README but the system I'm integrating this into would benefit greatly from being able to obtain this information automatically directly from the seg-y file.

Thanks!

[jokva]

Hi,

It depends on how grid spacing is recorded in the file, but I assume you can get this from the headers, or at least compute it from the coordinates, often in CDP X/Y or Source X/Y

[krober10nd]

Hey thanks, so I read in the file

 f = segyio.open("vel_z6.25m_x12.5m_exact.segy",ignore_geometry=True)

Print the header which I guess is a dictionary.

>>> f.header[0]
{TRACE_SEQUENCE_LINE: 0, TRACE_SEQUENCE_FILE: 0, FieldRecord: 1, TraceNumber: 1, EnergySourcePoint: 0, CDP: 1, CDP_TRACE: 1, TraceIdentificationCode: 1, NSummedTraces: 0, NStackedTraces: 0, DataUse: 0, offset: 0, ReceiverGroupElevation: 0, SourceSurfaceElevation: 0, SourceDepth: 0, ReceiverDatumElevation: 0, SourceDatumElevation: 0, SourceWaterDepth: 0, GroupWaterDepth: 0, ElevationScalar: -10, SourceGroupScalar: -10, SourceX: 125, SourceY: 0, GroupX: 125, GroupY: 0, CoordinateUnits: 0, WeatheringVelocity: 0, SubWeatheringVelocity: 0, SourceUpholeTime: 0, GroupUpholeTime: 0, SourceStaticCorrection: 0, GroupStaticCorrection: 0, TotalStaticApplied: 0, LagTimeA: 0, LagTimeB: 0, DelayRecordingTime: 0, MuteTimeStart: 0, MuteTimeEND: 0, TRACE_SAMPLE_COUNT: 1911, TRACE_SAMPLE_INTERVAL: 24080, GainType: 0, InstrumentGainConstant: 0, InstrumentInitialGain: 0, Correlated: 0, SweepFrequencyStart: 0, SweepFrequencyEnd: 0, SweepLength: 0, SweepType: 0, SweepTraceTaperLengthStart: 0, SweepTraceTaperLengthEnd: 0, TaperType: 0, AliasFilterFrequency: 0, AliasFilterSlope: 0, NotchFilterFrequency: 0, NotchFilterSlope: 0, LowCutFrequency: 0, HighCutFrequency: 0, LowCutSlope: 0, HighCutSlope: 0, YearDataRecorded: 1969, DayOfYear: 365, HourOfDay: 18, MinuteOfHour: 0, SecondOfMinute: 0, TimeBaseCode: 0, TraceWeightingFactor: 0, GeophoneGroupNumberRoll1: 0, GeophoneGroupNumberFirstTraceOrigField: 0, GeophoneGroupNumberLastTraceOrigField: 0, GapSize: 0, OverTravel: 0, CDP_X: 0, CDP_Y: 0, INLINE_3D: 0, CROSSLINE_3D: 0, ShotPoint: 0, ShotPointScalar: 0, TraceValueMeasurementUnit: 0, TransductionConstantMantissa: 0, TransductionConstantPower: 0, TransductionUnit: 0, TraceIdentifier: 0, ScalarTraceHeader: 0, SourceType: 0, SourceEnergyDirectionMantissa: 65536, SourceEnergyDirectionExponent: 0, SourceMeasurementMantissa: 0, SourceMeasurementExponent: 0, SourceMeasurementUnit: 0}
>>> 

Everything appears empty. Where can I find the Source X/Y?

[jokva]

It's a dictionary, yes. It has two fields, segyio.TraceField.SourceX and segyio.TraceField.SourceY, or segyio.su.sx and segyio.su.sy if you prefer. Your sources are in those two words, probably.

It does not appear empty, there are lots of values in there.

[krober10nd]

Thanks..so why does it error when I pass a key into the dictionary to get the values?

>>> h["CDP"]
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/Applications/miniconda3/lib/python3.7/site-packages/segyio/field.py", line 361, in __getitem__
    try: return self.getfield(self.buf, int(key))
ValueError: invalid literal for int() with base 10: 'CDP'

where h is the result of f.header[0]

[jokva]

Because it doesn't take strings, it take numerical constants. Try segyio.su.cdpx or segyio.TraceField.CDP_X.

[krober10nd]

Ok, I'm a little confused. When you say segyio you mean your package namespace or the open segy file?

f = segyio.open("vel_z6.25m_x12.5m_exact.segy",ignore_geometry=True)

The open segy file f in this case has no fields f.su.cdpx or f.TraceField.CDF_X and when I just type segyio.su.cdpx or segyio.TraceField.CDP_X. I just get constants that I think have nothing to do with the actual file.

[jokva]

Ok, I'm a little confused. When you say segyio you mean your package namespace or the open segy file?

Package namespace.

https://segyio.readthedocs.io/en/latest/segyio.html#constants

[krober10nd]

Ah I see. You get the number for the dictionary through this.

[jokva]

Yes.