From eb47898ae301b3812202abce3c1cb938760d7330 Mon Sep 17 00:00:00 2001 From: Robert Smallshire Date: Wed, 27 Aug 2014 08:47:18 +0200 Subject: [PATCH] Removes unused legacy file. --- segypy2.py | 1121 ---------------------------------------------------- 1 file changed, 1121 deletions(-) delete mode 100644 segypy2.py diff --git a/segypy2.py b/segypy2.py deleted file mode 100644 index 3c0f454..0000000 --- a/segypy2.py +++ /dev/null @@ -1,1121 +0,0 @@ -""" -A python module for reading/writing/manipuating -SEG-Y formatted filed - -segy.readSegy : Read SEGY file -segy.getSegyHeader : Get SEGY header -segy.getSegyTraceHeader : Get SEGY Trace header -segy.getAllSegyTraceHeaders : Get all SEGY Trace headers -segy.getSegyTrace : Get SEGY Trace heder and trace data for one trace - -segy.getValue : Get a value from a binary string -segy.ibm2ieee : Convert IBM floats to IEEE - -segy.version : The version of SegyPY -segy.verbose : Amount of verbose information to the scree -""" -# -# segypy : A Python module for reading and writing SEG-Y formatted data -# -# (C) Thomas Mejer Hansen, 2005 -# - - -import os -import struct -import logging - -import numpy -from numpy import zeros -from numpy import arange - - -from asq.initiators import query - -REEL_HEADER_NUM_BYTES = 3600 -TRACE_HEADER_NUM_BYTES = 240 - -logging.basicConfig(level=logging.INFO) -log = logging.getLogger("segypy") - -# SOME GLOBAL PARAMETERS -version = 0.2 -verbose = 1 - -endian = '>' # Big Endian -# endian='<' # Little Endian -# endian='=' # Native - -l_int = struct.calcsize('i') -l_uint = struct.calcsize('I') -l_long = struct.calcsize('l') -l_ulong = struct.calcsize('L') -l_short = struct.calcsize('h') -l_ushort = struct.calcsize('H') -l_char = struct.calcsize('c') -l_uchar = struct.calcsize('B') -l_float = struct.calcsize('f') - - -############## -# INIT - -############## -# Initialize SEGY HEADER -SH_def = {"Job": {"pos": 3200, "type": "int32", "def": 0}} -SH_def["Line"] = {"pos": 3204, "type": "int32", "def": 0} -SH_def["Reel"] = {"pos": 3208, "type": "int32", "def": 0} -SH_def["DataTracePerEnsemble"] = {"pos": 3212, "type": "int16", "def": 0} -SH_def["AuxiliaryTracePerEnsemble"] = {"pos": 3214, "type": "int16", "def": 0} -SH_def["dt"] = {"pos": 3216, "type": "uint16", "def": 1000} -SH_def["dtOrig"] = {"pos": 3218, "type": "uint16", "def": 1000} -SH_def["ns"] = {"pos": 3220, "type": "uint16", "def": 0} -SH_def["nsOrig"] = {"pos": 3222, "type": "uint16", "def": 0} -SH_def["DataSampleFormat"] = {"pos": 3224, "type": "int16", "def": 5} -SH_def["DataSampleFormat"]["descr"] = {0x0000: { - 1: "IBM Float", - 2: "32 bit Integer", - 3: "16 bit Integer", - 8: "8 bit Integer"}} - -SH_def["DataSampleFormat"]["descr"][0x0100] = { - 1: "IBM Float", - 2: "32 bit Integer", - 3: "16 bit Integer", - 5: "IEEE", - 8: "8 bit Integer"} - -SH_def["DataSampleFormat"]["bps"] = {0x0000: { - 1: 4, - 2: 4, - 3: 2, - 8: 1}} -SH_def["DataSampleFormat"]["bps"][0x0100] = { - 1: 4, - 2: 4, - 3: 2, - 5: 4, - 8: 1} -SH_def["DataSampleFormat"]["datatype"] = {0x0000: { - 1: 'ibm', - 2: 'l', - 3: 'h', - 8: 'B'}} -SH_def["DataSampleFormat"]["datatype"][0x0100] = { - 1: 'ibm', - 2: 'l', - 3: 'h', - # 5: 'float', - 5: 'f', - 8: 'B'} - -SH_def["EnsembleFold"] = {"pos": 3226, "type": "int16", "def": 0} -SH_def["TraceSorting"] = {"pos": 3228, "type": "int16", "def": 0} -SH_def["VerticalSumCode"] = {"pos": 3230, "type": "int16", "def": 0} -SH_def["SweepFrequencyEnd"] = {"pos": 3234, "type": "int16", "def": 0} -SH_def["SweepLength"] = {"pos": 3236, "type": "int16", "def": 0} -SH_def["SweepType"] = {"pos": 3238, "type": "int16", "def": 0} -SH_def["SweepChannel"] = {"pos": 3240, "type": "int16", "def": 0} -SH_def["SweepTaperlengthStart"] = {"pos": 3242, "type": "int16", "def": 0} -SH_def["SweepTaperLengthEnd"] = {"pos": 3244, "type": "int16", "def": 0} -SH_def["TaperType"] = {"pos": 3246, "type": "int16", "def": 0} -SH_def["CorrelatedDataTraces"] = {"pos": 3248, "type": "int16", "def": 0} -SH_def["BinaryGain"] = {"pos": 3250, "type": "int16", "def": 0} -SH_def["AmplitudeRecoveryMethod"] = {"pos": 3252, "type": "int16", "def": 0} -SH_def["MeasurementSystem"] = {"pos": 3254, "type": "int16", "def": 0} -SH_def["ImpulseSignalPolarity"] = {"pos": 3256, "type": "int16", "def": 0} -SH_def["VibratoryPolarityCode"] = {"pos": 3258, "type": "int16", "def": 0} -SH_def["Unassigned1"] = {"pos": 3260, "type": "int16", "n": 120, "def": 0} -SH_def["SegyFormatRevisionNumber"] = { - "pos": 3500, "type": "uint16", "def": 100} -SH_def["FixedLengthTraceFlag"] = {"pos": 3502, "type": "uint16", "def": 0} -SH_def["NumberOfExtTextualHeaders"] = {"pos": 3504, "type": "uint16", "def": 0} -SH_def["Unassigned2"] = {"pos": 3506, "type": "int16", "n": 47, "def": 0} - -############## -# Initialize SEGY TRACE HEADER SPECIFICATION -STH_def = {"TraceSequenceLine": {"pos": 0, "type": "int32"}} -STH_def["TraceSequenceFile"] = {"pos": 4, "type": "int32"} -STH_def["FieldRecord"] = {"pos": 8, "type": "int32"} -STH_def["TraceNumber"] = {"pos": 12, "type": "int32"} -STH_def["EnergySourcePoint"] = {"pos": 16, "type": "int32"} -STH_def["cdp"] = {"pos": 20, "type": "int32"} -STH_def["cdpTrace"] = {"pos": 24, "type": "int32"} -STH_def["TraceIdenitifactionCode"] = { - "pos": 28, "type": "uint16"} # 'int16'); % 28 -STH_def["TraceIdenitifactionCode"]["descr"] = {0x0000: { - 1: "Seismic data", - 2: "Dead", - 3: "Dummy", - 4: "Time Break", - 5: "Uphole", - 6: "Sweep", - 7: "Timing", - 8: "Water Break"}} -STH_def["TraceIdenitifactionCode"]["descr"][0x0100] = { - -1: "Other", - 0: "Unknown", - 1: "Seismic data", - 2: "Dead", - 3: "Dummy", - 4: "Time break", - 5: "Uphole", - 6: "Sweep", - 7: "Timing", - 8: "Waterbreak", - 9: "Near-field gun signature", - 10: "Far-field gun signature", - 11: "Seismic pressure sensor", - 12: "Multicomponent seismic sensor - Vertical component", - 13: "Multicomponent seismic sensor - Cross-line component", - 14: "Multicomponent seismic sensor - In-line component", - 15: "Rotated multicomponent seismic sensor - Vertical component", - 16: "Rotated multicomponent seismic sensor - Transverse component", - 17: "Rotated multicomponent seismic sensor - Radial component", - 18: "Vibrator reaction mass", - 19: "Vibrator baseplate", - 20: "Vibrator estimated ground force", - 21: "Vibrator reference", - 22: "Time-velocity pairs"} -STH_def["NSummedTraces"] = {"pos": 30, "type": "int16"} # 'int16'); % 30 -STH_def["NStackedTraces"] = {"pos": 32, "type": "int16"} # 'int16'); % 32 -STH_def["DataUse"] = {"pos": 34, "type": "int16"} # 'int16'); % 34 -STH_def["DataUse"]["descr"] = {0x0000: { - 1: "Production", - 2: "Test"}} -STH_def["DataUse"]["descr"][0x0100] = STH_def["DataUse"]["descr"][0] -STH_def["offset"] = {"pos": 36, "type": "int32"} # 'int32'); %36 -# 'int32'); %40 -STH_def["ReceiverGroupElevation"] = {"pos": 40, "type": "int32"} -# 'int32'); %44 -STH_def["SourceSurfaceElevation"] = {"pos": 44, "type": "int32"} -# 'int32'); %48 -STH_def["SourceDepth"] = {"pos": 48, "type": "int32"} -# 'int32'); %52 -STH_def["ReceiverDatumElevation"] = {"pos": 52, "type": "int32"} -# 'int32'); %56 -STH_def["SourceDatumElevation"] = {"pos": 56, "type": "int32"} -STH_def["SourceWaterDepth"] = {"pos": 60, "type": "int32"} # 'int32'); %60 -STH_def["GroupWaterDepth"] = {"pos": 64, "type": "int32"} # 'int32'); %64 -STH_def["ElevationScalar"] = {"pos": 68, "type": "int16"} # 'int16'); %68 -STH_def["SourceGroupScalar"] = {"pos": 70, "type": "int16"} # 'int16'); %70 -STH_def["SourceX"] = {"pos": 72, "type": "int32"} # 'int32'); %72 -STH_def["SourceY"] = {"pos": 76, "type": "int32"} # 'int32'); %76 -STH_def["GroupX"] = {"pos": 80, "type": "int32"} # 'int32'); %80 -STH_def["GroupY"] = {"pos": 84, "type": "int32"} # 'int32'); %84 -STH_def["CoordinateUnits"] = {"pos": 88, "type": "int16"} # 'int16'); %88 -STH_def["CoordinateUnits"]["descr"] = {0x0000: { - 1: "Length (meters or feet)", - 2: "Seconds of arc"}} -STH_def["CoordinateUnits"]["descr"][0x0100] = { - 1: "Length (meters or feet)", - 2: "Seconds of arc", - 3: "Decimal degrees", - 4: "Degrees, minutes, seconds (DMS)"} -STH_def["WeatheringVelocity"] = {"pos": 90, "type": "int16"} # 'int16'); %90 -STH_def["SubWeatheringVelocity"] = { - "pos": 92, "type": "int16"} # 'int16'); %92 -STH_def["SourceUpholeTime"] = {"pos": 94, "type": "int16"} # 'int16'); %94 -STH_def["GroupUpholeTime"] = {"pos": 96, "type": "int16"} # 'int16'); %96 -STH_def["SourceStaticCorrection"] = { - "pos": 98, "type": "int16"} # 'int16'); %98 -STH_def["GroupStaticCorrection"] = { - "pos": 100, "type": "int16"} # 'int16'); %100 -# 'int16'); %102 -STH_def["TotalStaticApplied"] = {"pos": 102, "type": "int16"} -STH_def["LagTimeA"] = {"pos": 104, "type": "int16"} # 'int16'); %104 -STH_def["LagTimeB"] = {"pos": 106, "type": "int16"} # 'int16'); %106 -# 'int16'); %108 -STH_def["DelayRecordingTime"] = {"pos": 108, "type": "int16"} -STH_def["MuteTimeStart"] = {"pos": 110, "type": "int16"} # 'int16'); %110 -STH_def["MuteTimeEND"] = {"pos": 112, "type": "int16"} # 'int16'); %112 -STH_def["ns"] = {"pos": 114, "type": "uint16"} # 'uint16'); %114 -STH_def["dt"] = {"pos": 116, "type": "uint16"} # 'uint16'); %116 -STH_def["GainType"] = {"pos": 119, "type": "int16"} # 'int16'); %118 -STH_def["GainType"]["descr"] = {0: { - 1: "Fixes", - 2: "Binary", - 3: "Floating point"}} -STH_def["GainType"]["descr"][1] = STH_def["GainType"]["descr"][0] -STH_def["InstrumentGainConstant"] = { - "pos": 120, "type": "int16"} # 'int16'); %120 -STH_def["InstrumentInitialGain"] = { - "pos": 122, "type": "int16"} # 'int16'); %%122 -STH_def["Correlated"] = {"pos": 124, "type": "int16"} # 'int16'); %124 -STH_def["Correlated"]["descr"] = {0: { - 1: "No", - 2: "Yes"}} -STH_def["Correlated"]["descr"][1] = STH_def["Correlated"]["descr"][0] - -# 'int16'); %126 -STH_def["SweepFrequenceStart"] = {"pos": 126, "type": "int16"} -STH_def["SweepFrequenceEnd"] = {"pos": 128, "type": "int16"} # 'int16'); %128 -STH_def["SweepLength"] = {"pos": 130, "type": "int16"} # 'int16'); %130 -STH_def["SweepType"] = {"pos": 132, "type": "int16"} # 'int16'); %132 -STH_def["SweepType"]["descr"] = {0x0000: { - 1: "linear", - 2: "parabolic", - 3: "exponential", - 4: "other"}} -STH_def["SweepType"]["descr"][0x0100] = STH_def["SweepType"]["descr"][0] - -STH_def["SweepTraceTaperLengthStart"] = { - "pos": 134, "type": "int16"} # 'int16'); %134 -STH_def["SweepTraceTaperLengthEnd"] = { - "pos": 136, "type": "int16"} # 'int16'); %136 -STH_def["TaperType"] = {"pos": 138, "type": "int16"} # 'int16'); %138 -STH_def["TaperType"]["descr"] = {0x0000: { - 1: "linear", - 2: "cos2c", - 3: "other"}} -STH_def["TaperType"]["descr"][0x0100] = STH_def["TaperType"]["descr"][0] - -STH_def["AliasFilterFrequency"] = { - "pos": 140, "type": "int16"} # 'int16'); %140 -STH_def["AliasFilterSlope"] = {"pos": 142, "type": "int16"} # 'int16'); %142 -STH_def["NotchFilterFrequency"] = { - "pos": 144, "type": "int16"} # 'int16'); %144 -STH_def["NotchFilterSlope"] = {"pos": 146, "type": "int16"} # 'int16'); %146 -STH_def["LowCutFrequency"] = {"pos": 148, "type": "int16"} # 'int16'); %148 -STH_def["HighCutFrequency"] = {"pos": 150, "type": "int16"} # 'int16'); %150 -STH_def["LowCutSlope"] = {"pos": 152, "type": "int16"} # 'int16'); %152 -STH_def["HighCutSlope"] = {"pos": 154, "type": "int16"} # 'int16'); %154 -STH_def["YearDataRecorded"] = {"pos": 156, "type": "int16"} # 'int16'); %156 -STH_def["DayOfYear"] = {"pos": 158, "type": "int16"} # 'int16'); %158 -STH_def["HourOfDay"] = {"pos": 160, "type": "int16"} # 'int16'); %160 -STH_def["MinuteOfHour"] = {"pos": 162, "type": "int16"} # 'int16'); %162 -STH_def["SecondOfMinute"] = {"pos": 164, "type": "int16"} # 'int16'); %164 -STH_def["TimeBaseCode"] = {"pos": 166, "type": "int16"} # 'int16'); %166 -STH_def["TimeBaseCode"]["descr"] = {0x0000: { - 1: "Local", - 2: "GMT", - 3: "Other"}} -STH_def["TimeBaseCode"]["descr"][0x0100] = { - 1: "Local", - 2: "GMT", - 3: "Other", - 4: "UTC"} -STH_def["TraceWeightningFactor"] = { - "pos": 168, "type": "int16"} # 'int16'); %170 -STH_def["GeophoneGroupNumberRoll1"] = { - "pos": 170, "type": "int16"} # 'int16'); %172 -STH_def["GeophoneGroupNumberFirstTraceOrigField"] = { - "pos": 172, "type": "int16"} # 'int16'); %174 -STH_def["GeophoneGroupNumberLastTraceOrigField"] = { - "pos": 174, "type": "int16"} # 'int16'); %176 -STH_def["GapSize"] = {"pos": 176, "type": "int16"} # 'int16'); %178 -STH_def["OverTravel"] = {"pos": 178, "type": "int16"} # 'int16'); %178 -STH_def["OverTravel"]["descr"] = {0x0000: { - 1: "down (or behind)", - 2: "up (or ahead)", - 3: "other"}} -STH_def["OverTravel"]["descr"][0x0100] = STH_def["OverTravel"]["descr"][0] - - -STH_def["cdpX"] = {"pos": 180, "type": "int32"} # 'int32'); %180 -STH_def["cdpY"] = {"pos": 184, "type": "int32"} # 'int32'); %184 -STH_def["Inline3D"] = {"pos": 188, "type": "int32"} # 'int32'); %188 -STH_def["Crossline3D"] = {"pos": 192, "type": "int32"} # 'int32'); %192 -STH_def["ShotPoint"] = {"pos": 192, "type": "int32"} # 'int32'); %196 -STH_def["ShotPointScalar"] = {"pos": 200, "type": "int16"} # 'int16'); %200 -STH_def["TraceValueMeasurementUnit"] = { - "pos": 202, "type": "int16"} # 'int16'); %202 -STH_def["TraceValueMeasurementUnit"]["descr"] = {0x0100: { - -1: "Other", - 0: "Unknown (should be described in Data Sample Measurement " - "Units Stanza) ", - 1: "Pascal (Pa)", - 2: "Volts (V)", - 3: "Millivolts (v)", - 4: "Amperes (A)", - 5: "Meters (m)", - 6: "Meters Per Second (m/s)", - 7: "Meters Per Second squared (m/&s2)Other", - 8: "Newton (N)", - 9: "Watt (W)"}} -STH_def["TransductionConstantMantissa"] = { - "pos": 204, "type": "int32"} # 'int32'); %204 -STH_def["TransductionConstantPower"] = { - "pos": 208, "type": "int16"} # 'int16'); %208 -STH_def["TransductionUnit"] = {"pos": 210, "type": "int16"} # 'int16'); %210 -STH_def["TransductionUnit"]["descr"] = STH_def[ - "TraceValueMeasurementUnit"]["descr"] -STH_def["TraceIdentifier"] = {"pos": 212, "type": "int16"} # 'int16'); %212 -STH_def["ScalarTraceHeader"] = {"pos": 214, "type": "int16"} # 'int16'); %214 -STH_def["SourceType"] = {"pos": 216, "type": "int16"} # 'int16'); %216 -STH_def["SourceType"]["descr"] = {0x0100: { - -1: "Other (should be described in Source Type/Orientation stanza)", - 0: "Unknown", - 1: "Vibratory - Vertical orientation", - 2: "Vibratory - Cross-line orientation", - 3: "Vibratory - In-line orientation", - 4: "Impulsive - Vertical orientation", - 5: "Impulsive - Cross-line orientation", - 6: "Impulsive - In-line orientation", - 7: "Distributed Impulsive - Vertical orientation", - 8: "Distributed Impulsive - Cross-line orientation", - 9: "Distributed Impulsive - In-line orientation"}} - -STH_def["SourceEnergyDirectionMantissa"] = { - "pos": 218, "type": "int32"} # 'int32'); %218 -STH_def["SourceEnergyDirectionExponent"] = { - "pos": 222, "type": "int16"} # 'int16'); %222 -STH_def["SourceMeasurementMantissa"] = { - "pos": 224, "type": "int32"} # 'int32'); %224 -STH_def["SourceMeasurementExponent"] = { - "pos": 228, "type": "int16"} # 'int16'); %228 -STH_def["SourceMeasurementUnit"] = { - "pos": 230, "type": "int16"} # 'int16'); %230 -STH_def["SourceMeasurementUnit"]["descr"] = {0x0100: { - -1: "Other (should be described in Source Measurement Unit stanza)", - 0: "Unknown", - 1: "Joule (J)", - 2: "Kilowatt (kW)", - 3: "Pascal (Pa)", - 4: "Bar (Bar)", - 4: "Bar-meter (Bar-m)", - 5: "Newton (N)", - 6: "Kilograms (kg)"}} -STH_def["UnassignedInt1"] = {"pos": 232, "type": "int32"} # 'int32'); %232 -STH_def["UnassignedInt2"] = {"pos": 236, "type": "int32"} # 'int32'); %236 - - -############## -# FUNCTIONS - -def image_segy(Data): - """ - imageSegy(Data) - Image segy Data - """ - import pylab - pylab.imshow(Data) - pylab.title('pymat test') - pylab.grid(True) - pylab.show() - - -def wiggle(Data, SH, skipt=1, maxval=8, lwidth=.1): - """ - wiggle(Data,SH) - """ - import pylab - - t = range(SH['ns']) - - for i in range(0, SH['ntraces'], skipt): - trace = Data[:, i] - trace[0] = 0 - trace[SH['ns'] - 1] = 0 - pylab.plot(i + trace / maxval, t, color='black', linewidth=lwidth) - for a in range(len(trace)): - if (trace[a] < 0): - trace[a] = 0 - # pylab.fill(i+Data[:,i]/maxval,t,color='k',facecolor='g') - pylab.fill(i + Data[:, i] / maxval, t, 'k', linewidth=0) - pylab.title(SH['filename']) - pylab.grid(True) - pylab.show() - - -def get_default_segy_header(ntraces=100, ns=100): - """ - SH=getDefaultSegyHeader() - """ - # INITIALIZE DICTIONARY - SH = {"Job": {"pos": 3200, "type": "int32", "def": 0}} - - for key in SH_def.keys(): - - tmpkey = SH_def[key] - if ('def' in tmpkey): - val = tmpkey['def'] - else: - val = 0 - SH[key] = val - - SH["ntraces"] = ntraces - SH["ns"] = ns - - return SH - - -def get_default_segy_trace_headers(ntraces=100, ns=100, dt=1000): - """ - SH=getDefaultSegyTraceHeader() - """ - # INITIALIZE DICTIONARY - STH = {"TraceSequenceLine": {"pos": 0, "type": "int32"}} - - for key in STH_def.keys(): - - # tmpkey = STH_def[key] - # if ('def' in tmpkey): - # val = tmpkey['def'] - # else: - # val = 0 - STH[key] = zeros(ntraces) - - for a in range(ntraces): - STH["TraceSequenceLine"][a] = a + 1 - STH["TraceSequenceFile"][a] = a + 1 - STH["FieldRecord"][a] = 1000 - STH["TraceNumber"][a] = a + 1 - STH["ns"][a] = ns - STH["dt"][a] = dt - return STH - - -def get_segy_trace_header(SH, THN='cdp', data=None): - """ - getSegyTraceHeader(SH,TraceHeaderName) - """ - logging.info("get trace header {0}".format(THN)) - bps = get_byte_per_sample(SH) - - if data is None: - with open(SH["filename"]) as f: - data = f.read() - - # MAKE SOME LOOKUP TABLE THAT HOLDS THE LOCATION OF HEADERS -# THpos=TraceHeaderPos[THN] - THpos = STH_def[THN]["pos"] - THformat = STH_def[THN]["type"] - ntraces = SH["ntraces"] - thv = zeros(ntraces) - getter = create_value_getter(data, THformat, endian, 1) - - for itrace in range(1, ntraces + 1, 1): - # i=itrace - - pos = THpos + 3600 + (SH["ns"] * bps + 240) * (itrace - 1) - - # txt="getSegyTraceHeader : Reading trace header " + THN + - # " " + str(itrace) + " of " + str(ntraces) + " " +str(pos) - - # printverbose(txt, logging.INFO); - thv[itrace - 1], index = getter(pos) - # txt="getSegyTraceHeader : " + THN + "=" + str(thv[itrace-1]) - # printverbose(txt,logging.DEBUG); - - return thv - - -def get_all_segy_trace_headers(SH, data='none'): - SegyTraceHeaders = {'filename': SH["filename"]} - - print_verbose( - 'getAllSegyTraceHeaders : trying to get all segy trace headers', - logging.INFO) - - if (data == 'none'): - data = open(SH["filename"]).read() - - for key in STH_def.keys(): - sth = get_segy_trace_header(SH, key, data) - SegyTraceHeaders[key] = sth - - txt = "getAllSegyTraceHeaders : " + key - print_verbose(txt, logging.DEBUG) - - return SegyTraceHeaders - - -def read_segy(filename): - """ - Data,SegyHeader,SegyTraceHeaders=getSegyHeader(filename) - """ - - print_verbose("readSegy : Trying to read " + filename, logging.DEBUG) - - with open(filename, 'rb') as raw: - data = raw.read() - - filesize = os.path.getsize(filename) - - segy_header = get_segy_header(filename) - - bytes_per_sample = get_byte_per_sample(segy_header) - - num_traces = (filesize - REEL_HEADER_NUM_BYTES) / \ - (segy_header['ns'] * bytes_per_sample + TRACE_HEADER_NUM_BYTES) - - print_verbose("readSegy : Length of data : " + str(filesize), logging.DEBUG) - - segy_header["ntraces"] = num_traces - - ndummy_samples = TRACE_HEADER_NUM_BYTES / bytes_per_sample - print_verbose( - "readSegy : ndummy_samples=" + str(ndummy_samples), logging.DEBUG) - print_verbose("readSegy : ntraces=" + str(num_traces) + - " nsamples=" + str(segy_header['ns']), logging.DEBUG) - - # GET TRACE - index = REEL_HEADER_NUM_BYTES - num_pseudo_samples = (filesize - REEL_HEADER_NUM_BYTES) / bytes_per_sample - - # READ ALL SEGY TRACE HEADRES - # Temporarily removed - # SegyTraceHeaders = getAllSegyTraceHeaders(segy_header,data) - - inline_3d_trace_header = get_segy_trace_header(segy_header, "Inline3D", data) - crossline_3d_trace_header = get_segy_trace_header( - segy_header, "Crossline3D", data) - num_samples_trace_header = get_segy_trace_header(segy_header, "ns", data) - - num_inlines = query(inline_3d_trace_header).distinct().count() - num_crosslines = query(crossline_3d_trace_header).distinct().count() - num_samples_per_trace = query(num_samples_trace_header).max() - # TODO: This assumes the number of samples per trace is constant - num_pseudo_samples_per_trace = num_samples_per_trace + \ - TRACE_HEADER_NUM_BYTES - - if num_inlines * num_crosslines != num_traces: - raise RuntimeError("inlines * crosslines != traces") - - if num_inlines * num_crosslines * num_pseudo_samples_per_trace \ - != num_pseudo_samples: - raise RuntimeError("num_inlines * num_crosslines * num_samples != nd") - - print_verbose("readSegy : reading segy data", logging.DEBUG) - - # READ ALL DATA EXCEPT FOR SEGY HEADER - # Data = zeros((SH['ns'],ntraces)) - - revision = segy_header["SegyFormatRevisionNumber"] - # if (revision==0x100): - # revision=1 - dsf = segy_header["DataSampleFormat"] - - DataDescr = SH_def["DataSampleFormat"]["descr"][revision][dsf] - - print_verbose("readSegy : SEG-Y revision = " + str(revision), logging.INFO) - print_verbose("readSegy : DataSampleFormat=" + str(dsf) + - "(" + DataDescr + ")", logging.INFO) - - if (segy_header["DataSampleFormat"] == 1): - print_verbose("readSegy : Assuming DSF=1, IBM FLOATS", logging.INFO) - pseudo_samples = get_value( - data, index, 'ibm', endian, num_pseudo_samples) - elif (segy_header["DataSampleFormat"] == 2): - print_verbose("readSegy : Assuming DSF=" + - str(segy_header["DataSampleFormat"]) + ", 32bit INT", - logging.INFO) - pseudo_samples = get_value(data, index, 'l', endian, num_pseudo_samples) - elif (segy_header["DataSampleFormat"] == 3): - print_verbose("readSegy : Assuming DSF=" + - str(segy_header["DataSampleFormat"]) + ", 16bit INT", - logging.INFO) - pseudo_samples = get_value(data, index, 'h', endian, num_pseudo_samples) - elif (segy_header["DataSampleFormat"] == 5): - print_verbose("readSegy : Assuming DSF=" + - str(segy_header["DataSampleFormat"]) + ", IEEE", - logging.INFO) - pseudo_samples = get_value( - data, index, 'float', endian, num_pseudo_samples) - elif (segy_header["DataSampleFormat"] == 8): - print_verbose("readSegy : Assuming DSF=" + - str(segy_header["DataSampleFormat"]) + ", 8bit CHAR", - logging.INFO) - with open(filename, 'rb') as f: - f.seek(index) - pseudo_samples = numpy.fromfile( - f, numpy.int8, num_pseudo_samples) - # Data1 = getValue(data,index,'B',endian,nd) - else: - print_verbose( - "readSegy : DSF=" + - str(segy_header["DataSampleFormat"]) + - ", NOT SUPPORTED", - logging.INFO) - - # Reshape the 1D pseudo_samples array into a 3D array (including - # the trace headers) - # 1. pseudo_data_3d = numpy.reshape(pseudo_samples, (num_crosslines, - # num_inlines, num_pseudo_samples_per_trace)) - pseudo_data_3d = numpy.reshape( - pseudo_samples, (num_inlines, - num_crosslines, - num_pseudo_samples_per_trace)) - - # Slice of the dummy trace header data - sliced_data_3d = pseudo_data_3d[:, :, TRACE_HEADER_NUM_BYTES:] - - # Transpose - # data_3d = numpy.transpose(sliced_data_3d, (0, 1, 2)) - signed_data_3d = numpy.transpose(sliced_data_3d, (0, 2, 1)) - # data_3d = numpy.transpose(sliced_data_3d, (1, 0, 2)) - # data_3d = numpy.transpose(sliced_data_3d, (1, 2, 0)) - # data_3d = numpy.transpose(sliced_data_3d, (2, 0, 1)) - # data_3d = numpy.transpose(sliced_data_3d, (2, 1, 0)) - - # data_3d = pseudo_data_3d - - # Data = pseudo_samples[0] - # - # printverbose("readSegy : - reshaping", logging.INFO) - # Data=reshape(Data,(num_traces,segy_header['ns']+ndummy_samples)) - # printverbose("readSegy : - stripping header dummy data", logging.INFO) - # Data=Data[:,ndummy_samples:(segy_header['ns']+ndummy_samples)] - # printverbose("readSegy : - transposing", logging.INFO) - # Data=transpose(Data) - - # SOMEONE NEEDS TO IMPLEMENT A NICER WAY DO DEAL WITH DSF=8 - # if segy_header["DataSampleFormat"] == 8: - # for i in arange(num_traces): - # for j in arange(segy_header['ns']): - # if Data[i][j]>128: - # Data[i][j]=Data[i][j]-256 - - print_verbose("readSegy : read data", logging.DEBUG) - SegyTraceHeaders = None # Temp - return signed_data_3d, segy_header, SegyTraceHeaders - - -def get_segy_trace(SH, itrace): - """ - SegyTraceHeader,SegyTraceData=getSegyTrace(SegyHeader,itrace) - itrace : trace number to read - THIS DEF IS NOT UPDATED. NOT READY TO USE - """ - data = open(SH["filename"]).read() - - bps = get_byte_per_sample(SH) - - # GET TRACE HEADER - index = 3200 + (itrace - 1) * (240 + SH['ns'] * bps) - SegyTraceHeader = [] - # print index - - # GET TRACE - index = 3200 + (itrace - 1) * (240 + SH['ns'] * bps) + 240 - SegyTraceData = get_value(data, index, 'float', endian, SH['ns']) - return SegyTraceHeader, SegyTraceData - - -def get_segy_header(filename): - """ - SegyHeader=getSegyHeader(filename) - """ - with open(filename, 'rb') as f: - data = f.read() - - SegyHeader = {'filename': filename} - for key in SH_def.keys(): - pos = SH_def[key]["pos"] - format = SH_def[key]["type"] - - SegyHeader[key], index = get_value(data, pos, format, endian) - - txt = str(pos) + " " + str(format) + " Reading " + \ - key + "=" + str(SegyHeader[key]) - print_verbose(txt, logging.INFO) - - # SET NUMBER OF BYTES PER DATA SAMPLE - bps = get_byte_per_sample(SegyHeader) - - filesize = len(data) - ntraces = (filesize - 3600) / (SegyHeader['ns'] * bps + 240) - SegyHeader["ntraces"] = ntraces - - print_verbose('getSegyHeader : succesfully read ' + filename, logging.INFO) - - return SegyHeader - - -def write_segy(filename, Data, dt=.001, STHin={}): - """ - writeSegy(filename,Data,dt) - - Write SEGY - - See also readSegy - - (c) 2005, Thomas Mejer Hansen - - MAKE OPTIONAL INPUT FOR ALL SEGYHTRACEHEADER VALUES - - """ - - print_verbose("writeSegy : Trying to write " + filename, logging.DEBUG) - - N = Data.shape - ns = N[0] - ntraces = N[1] - print ntraces, ns - - SH = get_default_segy_header(ntraces, ns) - STH = get_default_segy_trace_headers(ntraces, ns, dt) - - # ADD STHin, if exists... - for key in STHin.keys(): - print key - for a in range(ntraces): - STH[key] = STHin[key][a] - - write_segy_structure(filename, Data, SH, STH) - - -def write_segy_structure(filename, Data, SH, STH): - """ - writeSegyHeader(filename,Data,SegyHeader,SegyTraceHeaders) - - Write SEGY file using SegyPy data structures - - See also readSegy - - (c) 2005, Thomas Mejer Hansen - - """ - - print_verbose( - "writeSegyStructure : Trying to write " + filename, logging.DEBUG) - - f = open(filename, 'w') - - # VERBOSE INF - revision = SH["SegyFormatRevisionNumber"] - dsf = SH["DataSampleFormat"] - # if (revision==100): - # revision=1 - DataDescr = SH_def["DataSampleFormat"]["descr"][revision][dsf] - - print_verbose( - "writeSegyStructure : SEG-Y revision = " + str(revision), logging.INFO) - print_verbose("writeSegyStructure : DataSampleFormat=" + - str(dsf) + "(" + DataDescr + ")", logging.INFO) - - # WRITE SEGY HEADER - - for key in SH_def.keys(): - pos = SH_def[key]["pos"] - format = SH_def[key]["type"] - value = SH[key] - -# SegyHeader[key],index = putValue(value,f,pos,format,endian); - put_value(value, f, pos, format, endian) - - txt = str(pos) + " " + str(format) + \ - " Reading " + key + "=" + str(value) -# +"="+str(SegyHeader[key]) - # printverbose(txt,-1) - - # SEGY TRACES - - ctype = SH_def['DataSampleFormat']['datatype'][revision][dsf] - bps = SH_def['DataSampleFormat']['bps'][revision][dsf] - - sizeT = 240 + SH['ns'] * bps - - for itrace in range(SH['ntraces']): - index = 3600 + itrace * sizeT - print_verbose( - 'Writing Trace #' + - str(itrace + 1) + - '/' + - str(SH['ntraces']), - logging.INFO) - # WRITE SEGY TRACE HEADER - for key in STH_def.keys(): - pos = index + STH_def[key]["pos"] - format = STH_def[key]["type"] - value = STH[key][itrace] - txt = str(pos) + " " + str(format) + \ - " Writing " + key + "=" + str(value) - print_verbose(txt, logging.INFO) - put_value(value, f, pos, format, endian) - - # Write Data - cformat = endian + ctype - for s in range(SH['ns']): - strVal = struct.pack(cformat, Data[s, itrace]) - f.seek(index + 240 + s * struct.calcsize(cformat)) - f.write(strVal) - - f.close - - # return segybuffer - - -def put_value(value, fileid, index, ctype='l', endian='>', number=1): - """ - putValue(data,index,ctype,endian,number) - """ - if (ctype == 'l') | (ctype == 'long') | (ctype == 'int32'): - # size = l_long - ctype = 'l' - elif (ctype == 'L') | (ctype == 'ulong') | (ctype == 'uint32'): - # size = l_ulong - ctype = 'L' - elif (ctype == 'h') | (ctype == 'short') | (ctype == 'int16'): - # size = l_short - ctype = 'h' - elif (ctype == 'H') | (ctype == 'ushort') | (ctype == 'uint16'): - # size = l_ushort - ctype = 'H' - elif (ctype == 'c') | (ctype == 'char'): - # size = l_char - ctype = 'c' - elif (ctype == 'B') | (ctype == 'uchar'): - # size = l_uchar - ctype = 'B' - elif (ctype == 'f') | (ctype == 'float'): - # size = l_float - ctype = 'f' - elif (ctype == 'ibm'): - # size = l_float - pass - else: - print_verbose('Bad Ctype : ' + ctype, logging.ERROR) - - cformat = endian + ctype * number - - print_verbose( - 'putValue : cformat : ' + cformat + ' ctype=' + ctype, logging.DEBUG) - - strVal = struct.pack(cformat, value) - fileid.seek(index) - fileid.write(strVal) - - return 1 - - -def create_value_getter(data, ctype='l', endian='>', number=1): - """ - getValue(data,index,ctype,endian,number) - """ - if (ctype == 'l') | (ctype == 'long') | (ctype == 'int32'): - size = l_long - ctype = 'l' - elif (ctype == 'L') | (ctype == 'ulong') | (ctype == 'uint32'): - size = l_ulong - ctype = 'L' - elif (ctype == 'h') | (ctype == 'short') | (ctype == 'int16'): - size = l_short - ctype = 'h' - elif (ctype == 'H') | (ctype == 'ushort') | (ctype == 'uint16'): - size = l_ushort - ctype = 'H' - elif (ctype == 'c') | (ctype == 'char'): - size = l_char - ctype = 'c' - elif (ctype == 'B') | (ctype == 'uchar'): - size = l_uchar - ctype = 'B' - elif (ctype == 'f') | (ctype == 'float'): - size = l_float - ctype = 'f' - elif (ctype == 'ibm'): - size = l_float - else: - print_verbose('Bad Ctype : ' + ctype, logging.ERROR) - - cformat = endian + ctype * number - - # printverbose('getValue : cformat : ' + cformat, logging.DEBUG) - - if (ctype == 'ibm'): - # ASSUME IBM FLOAT DATA - - def ibm(index): - index_end = index + size * number - Value = range(number) - for i in arange(number): - index_ibm = i * 4 + index - Value[i] = ibm2ieee2(data[index_ibm:index_ibm + 4]) - # this return an array as opposed to a tuple - if number == 1: - return Value[0], index_end - - return Value, index_end - - return ibm - - else: - if (ctype == 'B'): - print_verbose( - 'getValue : Inefficient use of 1 byte Integer...', - logging.WARNING) - - # ALL OTHER TYPES OF DATA - def other(index): - index_end = index + size * number - Value = struct.unpack(cformat, data[index:index_end]) - - if number == 1: - return Value[0], index_end - - return Value, index_end - - return other - - -def get_value(data, index, ctype='l', endian='>', number=1): - """ - getValue(data,index,ctype,endian,number) - """ - if (ctype == 'l') | (ctype == 'long') | (ctype == 'int32'): - size = l_long - ctype = 'l' - elif (ctype == 'L') | (ctype == 'ulong') | (ctype == 'uint32'): - size = l_ulong - ctype = 'L' - elif (ctype == 'h') | (ctype == 'short') | (ctype == 'int16'): - size = l_short - ctype = 'h' - elif (ctype == 'H') | (ctype == 'ushort') | (ctype == 'uint16'): - size = l_ushort - ctype = 'H' - elif (ctype == 'c') | (ctype == 'char'): - size = l_char - ctype = 'c' - elif (ctype == 'B') | (ctype == 'uchar'): - size = l_uchar - ctype = 'B' - elif (ctype == 'f') | (ctype == 'float'): - size = l_float - ctype = 'f' - elif (ctype == 'ibm'): - size = l_float - else: - print_verbose('Bad Ctype : ' + ctype, logging.ERROR) - - cformat = endian + ctype * number - - # printverbose('getValue : cformat : ' + cformat, logging.DEBUG) - - index_end = index + size * number - - if (ctype == 'ibm'): - # ASSUME IBM FLOAT DATA - Value = range(number) - for i in arange(number): - index_ibm = i * 4 + index - Value[i] = ibm2ieee2(data[index_ibm:index_ibm + 4]) - # this resturn an array as opposed to a tuple - else: - # ALL OTHER TYPES OF DATA - Value = struct.unpack(cformat, data[index:index_end]) - - if (ctype == 'B'): - print_verbose( - 'getValue : Inefficient use of 1 byte Integer...', logging.WARNING) - - # vtxt = 'getValue : '+'start='+str(index)+' size='+str(size)+ ' - # number='+str(number)+' Value='+str(Value)+' - # cformat='+str(cformat) - # printverbose(vtxt, logging.DEBUG) - - if number == 1: - return Value[0], index_end - - return Value, index_end - - -def print_version(): - print 'SegyPY version is ', version - - -def print_verbose(txt, level=logging.DEBUG): - pass - # log.log(level, txt) - - -############## -# MISC FUNCTIONS -def ibm2Ieee(ibm_float): - """ - ibm2Ieee(ibm_float) - Used by permission - (C) Secchi Angelo - with thanks to Howard Lightstone and Anton Vredegoor. - """ - i = struct.unpack('>I', ibm_float)[0] - sign = [1, -1][bool(i & 0x100000000L)] - characteristic = ((i >> 24) & 0x7f) - 64 - fraction = (i & 0xffffff) / float(0x1000000L) - return sign * 16 ** characteristic * fraction - - -def ibm2ieee2(ibm_float): - """ - ibm2ieee2(ibm_float) - Used by permission - (C) Secchi Angelo - with thanks to Howard Lightstone and Anton Vredegoor. - """ - dividend = float(16 ** 6) - - if ibm_float == 0: - return 0.0 - istic, a, b, c = struct.unpack('>BBBB', ibm_float) - if istic >= 128: - sign = -1.0 - istic = istic - 128 - else: - sign = 1.0 - mant = float(a << 16) + float(b << 8) + float(c) - return sign * 16 ** (istic - 64) * (mant / dividend) - - -def get_byte_per_sample(SH): - revision = SH["SegyFormatRevisionNumber"] - - # if revision == 100: - # revision=1 - - log.debug("revision = {0}".format(revision)) - - dsf = SH["DataSampleFormat"] - - log.debug("data sample format = {0}".format(dsf)) - - bps = SH_def["DataSampleFormat"]["bps"][revision][dsf] - - log.debug("bytes per sample = {0}".format(bps)) - - print_verbose("getBytePerSample : bps=" + str(bps), logging.DEBUG) - - return bps - - -############## -# segy class -class SegyTraceheaderClass: - - def __init__(self): - self.cdp = 0 - - -class SegyHeaderClass: - - def __str__(self): - return "SegyHeaderClass " - - def __init__(self): - self.filename = 0 - self.Trace = version - - def cdp(self): - return "Getting CDP trace header" - - def inline_x(self): - return "Getting CDP trace header" - - -class SegyClass: - STH_def = STH_def - SH_def = SH_def - STH = SegyTraceheaderClass() - SH = SegyHeaderClass() - - def __init__(self): - self.THOMAS = 'Thomas' - -if __name__ == '__main__': - signed_data_3d, sh, sth = read_segy( - r'C:\Users\rjs\opendtectroot\Blake_Ridge_Hydrates_3D' - r'\stack_final_scaled50_int8.sgy') - unsigned_data_3d = signed_data_3d.view(numpy.uint8) - unsigned_data_3d += 128 - numpy.save('blake_ridge.npy', unsigned_data_3d)