#!/usr/bin/env python
"""
Module to read segy files. No writing at this time.
Copyright (C) 2005 Kurt Schwehr
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
$Id: segy.py,v 1.13 2005/05/26 14:12:28 schwehr Exp $
Eventually, each segy version (e.g. xstar/edgetech, and other
instrument tweaks) will be loaded on the fly. So if a field is
somewhere else, it will be easy to tweak a new config in just a few
minutes
"""
import mmap
import os
import struct # Convert two and from binary formats
import unittest
traceHdrSize = 240 # bytes
# These are the format characters for integers with the struct module for binary conversions
# Force big endian (default on osx, right?)
structIntCodes = { 1: '>b', 2: '>h', 4: '>i', 8: '>q' }
# Name -> entry number and starting byte position
# each entry is 80 bytes long
# name, line, line offset (note: does not start from 0!!!!)
textFileHeaderEntries = {
'CLIENT': [1,5],
'COMPANY': [1,35],
'CREW NO': [1,65],
'LINE': [2,5],
'AREA': [2,11],
'MAP ID': [2,49],
'REAL NO': [3,5],
}
dataFormats = {
1: "4-byte IBM floating point",
2: "4-byte, two's complement integer",
3: "2-byte, two's complement integer",
4: "4-byte, fixed-point with gain(obsolete)",
5: "4-byte IEEE floating-point",
6: "Not currently used",
7: "Not currently used",
8: "1-byte, two's complement integer"
}
# module-struct code. See module-struct.html for the call to get size
dataFormatStruct = {
1: None, # What to do for IBM floats?
2: structIntCodes[4],
3: structIntCodes[2],
4: None, # structintCode[4], * something...
5: "f", # float
6: None,
7: None,
8: structIntCodes[1],
}
traceSortingCodes = {
# Page 10
-1: "Other",
0: "Unknown",
1: "As recorded", # No sorting
2: "CDP ensemble",
3: "Single fold continuous profile",
4: "Horizontally stacked",
5: "Common source point",
6: "Common receiver point",
7: "Common offset point",
8: "Common mid-point",
9: "Common conversion point"
}
sweetTypeCodes = {
1: "linear",
2: "parabolic",
3: "exponential",
4: "other"
}
# Use the index from the segy r1 specification
binaryHeaderEntries = {
'JobId': [3201,3204,'Job identification number'],
'LineNo': [3205,3208,'Line number'],
'ReelNo': [3209,3212,'Reel number'],
'TracesPerEnsemble': [3213,3214,'Number of data traces per ensemble'],
'AuxTracesPerEnsemble': [3215,3216,'Number of auxiliary traces per ensemble'],
'SampleInterval': [3217,3218, 'Sample interval in microseconds'],
'OrigSampleInterval': [3219,3220, 'Sample interval in microsends of original field recording'],
'SamplesPerTrace': [3221,3222, 'Number of samples per data trace'],
'OrigSamplesPerTrace': [3223,3224, 'Number of samples per data trace for original field record'],
'SampleFormat': [3225,3226,'Data sample format code'],
'EnsembleFold': [3227,3228,'The expected number of data traces per trace ensemble - the CMP fold'],
'TraceSorting': [3229,3230,'The sorting code or type of ensemble'],
'VerticalSum': [3231,3232,'Vertical sum code N=M-1 sum'],
'SweepStart': [3233,3234,'Sweep frequence at start (Hz)'],
'SweepEnd': [3235,3236,'Sweep frequence at end (Hz)'],
'SweepLen': [3237,3238,'Sweep length (ms)'],
'SweepType': [3239,3240,'Sweep type code'],
'TraceSweepChannel': [3241,3242, "Trace number of sweep channel"],
'SweepTaperLenStart': [3243,3244,'Sweep trace taper length in ms at start if tapered'],
'SweepTaperLenStart': [3245,3246,'Sweep trace taper length in ms at end if tapered'],
'TaperType': [3247,3248,'Taper type'],
'CorrelatedTraces': [3249,3250,'Correlated data traces. 1=no, 2=yes'],
# NOTICE: yes and no are different between the two!
'BinGainRecovered': [3251,3252,'Binary gain recovered. 1=yes, 2=no'],
'AmpRecovMeth': [3253,3254,'Amplitude recovery method'],
'MeasurementSystem': [3255,3256,'1 = Meters, 2 = Feet'],
'ImpulsePolarity': [3257,3258,'Impulse signal polarity. increasing pres or move give negative value'],
'VibPolarity': [3259,3260,'Vibratory polarity code'],
#'Unassigned3261': [3261,3500, 'Unassigned'],
'SegYRev': [3501,3502,'SEG Y Format Revision Number. Rev 1.0 is 0100'],
'FixedLenTraceFlag': [3503,3504,'1 indeicates that all trances in this file are all the same length'],
'NumExtendedTextHeaders': [3505,3506,'Number of 3200B ext text headers after the bin header'],
#'Unassigned3507': [3507,3600,'Unassigned']
}
class Segy:
"""
Base class that contains the SEGY header and all the traces
http://seg.org/publications/tech-stand/
Then look at segy rev 0 and rev 1
"""
def __init__(self,filename=""):
"""
Open up a segy file for reading using mmap for speed.
"""
file = open(filename,"r+")
self.size = os.path.getsize(filename)
self.data = mmap.mmap(file.fileno(),self.size,access=mmap.ACCESS_READ)
file.close() # No longer need the file handle
self.header=SegyHdr(self.data,self.size)
return
def getHeader(self):
return self.header
def getTraceHdr(self,traceNumber=1):
if 1==traceNumber:
return SegyTraceHdr(segyHdr=self.header,mmapData=self.header.data,startOffset=3600)
extTextSize = 0 # FIX!!!
traceStart = 3600 + extTextSize + (traceNumber-1)*self.getTraceSize()
return SegyTraceHdr(segyHdr=self.header,mmapData=self.header.data,startOffset=traceStart)
def getTraceDataSize(self,traceNumber=1):
""" FIX: handle variable length traces"""
format = self.header.getBinEntry('SampleFormat')
f = dataFormatStruct[format]
size = struct.calcsize(f)
samples = self.header.getBinEntry('SamplesPerTrace')
datasize = size * samples
return datasize
def getTraceSize(self, traceNumber=1):
""" returns the number of bytes for the particular trace
fix: cope with variable length traces
"""
return 240 + self.getTraceDataSize()
def getNumberOfTraces(self):
""" FIX: deal with variable length traces """
extTextSize = 0 # FIX
tracesSize = self.size - 3600 - extTextSize
numTraces = tracesSize / self.getTraceSize()
if 0 != (tracesSize % self.getTraceSize()):
print "WARNING: seems you have variable length traces that will not be handled!!!"
assert(False)
return numTraces
def getTraceData(self,traceNumber=1):
""" return a list of samples
FIX: cope with variable length traces
"""
extTextSize = 0 # Calc this!
traceStart = 3600 + extTextSize + (traceNumber-1)*self.getTraceSize()
print "traceStart = ", traceStart, traceNumber, self.getTraceSize(), (traceNumber-1)*self.getTraceSize()
startData = 240 + traceStart
endData = startData+self.getTraceDataSize()
format = self.header.getBinEntry('SampleFormat')
formatStr = dataFormatStruct[format]
parseStr = formatStr[0] + formatStr[1] * self.header.getBinEntry('SamplesPerTrace')
rawData = self.data[startData:endData]
traceData = struct.unpack(parseStr,rawData)
return traceData
def writeTraceToFile(self,filename='default.dat',traceNumber=1):
""" Write out pairs of samplenumber and value to a file
"""
data = self.getTraceData(traceNumber)
o = open(filename,"w")
for i in range(0,len(data)):
o.write(str(i)+" "+str(data[i])+"\n")
o.close()
######################################################################
class SegyHdr:
""" class for reading segy seismic file header
http://seg.org/publications/tech-stand/
Then look at segy rev 0 and rev 1
"""
def __init__(self,mmapData=None,size=None, filename=None):
"""
Open up a segy file for reading using mmap for speed.
"""
if filename and (not (mmapData and size)):
file = open(filename,"r+")
self.size = os.path.getsize(filename)
self.data = mmap.mmap(file.fileno(),self.size,access=mmap.ACCESS_READ)
file.close() # No longer need the file handle
return
assert(mmapData and size)
self.data=mmapData
self.size=size
return
# FIX: move this to Segy class
def hasR1TextFileHeader(self):
"""Turn True if it has a R1 header with C1..C40 sections as ASCII"""
if self.data[0:2] == 'C1' or self.data[0:3] == 'C 1':
return True
return False
# FIX: move this to Segy class
def getTextHeaderLine(self,lineNumber):
""" Return the string for the line number """
line = self.data[(lineNumber-1)*80:(lineNumber)*80]
nul = line.find("\0")
if -1 == nul: return line
return line[0:nul]
# FIX: move this to Segy class or make a text header class
def getTextHeaderField(self, fieldName):
"""
BROKEN
e.g.:
C2 LINE:5t6 AREA:Goleta
"""
print "FIX: broken"
entry = textFileHeaderEntries[fieldName]
line = self.getTextHeaderLine(entry[0])
if len(line) < entry[1]:
print "Field "+fieldName+" not present"
return None
print "FIX: figure out the end of the entry"
return line[entry[1]:]
def getBinEntry(self,fieldName):
""" FIX: explain yourself! """
entry = binaryHeaderEntries[fieldName]
start = entry[0] - 1
end = entry[1] - 1
length = entry[1]-entry[0] + 1
if length > 8:
print "get bin header... too large. skipping ",fieldName
return None
val = struct.unpack(structIntCodes[length],self.data[start:end+1])[0]
return val
def printBinaryHeaders(self):
""" Dump all header into to stdout"""
for key in binaryHeaderEntries.keys():
val = self.getBinEntry(key)
print key, " = ",val
key = 'SampleFormat'; print key,' => ',dataFormats[self.getBinEntry(key)]
######################################################################
# Trace parsing
######################################################################
traceIdCode = {
-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'
#23 - N = optional use, (maximum N = 32,767)
}
dataUse = {
1: 'Production',
2: 'Test',
}
coordUnits = {
1: 'Length (meters or feet)',
# An arc second is therefore an angular measurement that equals
# 1/60th of an arc minute, or 1/3600 of a single degree of arc.
2: 'Seconds of arc',
3: 'Decimal degrees',
4: 'Degrees, minutes, seconds (DMS)Note: To encode +-DDDMMSS bytes 89-90 equal = +-DDD*104 + MM*102 + SS with bytes 71-72 set to 1; To encode +-DDDMMSS.ss bytes 89-90 equal = +-DDD*106 + MM*104 + SS*102 with bytes 71-72 set to -100.',
}
gainType = {
1: 'fixed',
2: 'binary',
3: 'floating point'
#4 ... N = optional use
}
sweepType = {
1: 'linear ',
2: 'parabolic',
3: 'exponential',
4: 'other',
}
taperType = {
1: 'linear ',
2: 'cos2',
3: 'other'
}
timeBasis = {
1: 'Local',
2: 'GMT (Greenwich Mean Time)',
3: 'Other, should be explained in a user defined stanza in the Extended Textual File Header',
4: 'UTC (Coordinated Universal Time)'
}
overTravel = {
1: 'down (or behind) ',
2: 'up (or ahead)'
}
traceUnits = {
-1: 'Other (should be described in Data Sample Measurement Units Stanza)',
0: 'Unknown',
1: 'Pascal (Pa)',
2: 'Volts (v)',
3: 'Millivolts (mV)',
4: 'Amperes (A)',
5: 'Meters (m)',
6: 'Meters per second (m/s)',
7: 'Meters per second squared (m/s2)',
8: 'Newton (N)',
9: 'Watt (W)'
}
transdUnits = {
-1: 'Other (should be described in Data Sample Measurement Unit stanza, page 36)',
0: 'Unknown',
1: 'Pascal (Pa)',
2: 'Volts (v)',
3: 'Millivolts (mV)',
4: 'Amperes (A)',
5: 'Meters (m)',
6: 'Meters per second (m/s)',
7: 'Meters per second squared (m/s2)',
8: 'Newton (N)',
9: 'Watt (W)',
}
srcTypeOrient = {
#-1 to -n = Other (should be described in Source Type/Orientation stanza, page 38)'
-1:'Other (should be described in Source Type/Orientation stanza, page 38)',
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',
}
srcMeasUnits = {
-1: 'Other (should be described in Source Measurement Unit stanza, page 39)',
0: 'Unknown',
1: 'Joule (J)',
2: 'Kilowatt (kW)',
3: 'Pascal (Pa)',
4: 'Bar (Bar)',
5: 'Bar-meter (Bar-m)', # FIX: What the??? This is not right!!!!! WTK this was listed as 4 too
6: 'Newton (N)',
7: 'Kilograms (kg)'
}
traceHeader = {
'LineSeqNo': [1,4,'Trace sequence number within line - Numbers continue to increase if the same line continues across multiple SEG Y files. Highly recommended for all types of data.'],
'FileSeqNo': [5,8,'Trace sequence number within SEG Y file - Each file starts with trace sequence one.'],
'FieldRecNo': [9,12,'Original field record number. Highly recommended for all types of data.'],
'TraceNo': [13,16,'Trace number within the original field record. Highly recommended for all types of data.'],
'ESrcPtNo': [17,20,'Energy source point number - Used when more than one record occurs at the same effective surface location. It is recommended that the new entry defined in Trace Header bytes 197-202 be used for shotpoint number.'],
'EnsembleNo': [21,24,'Ensemble number (i.e. CDP, CMP, CRP, etc)'],
'TraceNoEnsemble': [25,28,'Trace number within the ensemble - Each ensemble starts with trace number one.'],
'TraceId': [29,30,'Trace identification code'],
'VertSumNo': [31,32,'Number of vertically summed traces yielding this trace. (1 is one trace, 2 is two summed traces, etc.)'],
'HorzStackedNo': [33,34,'Number of horizontally stacked traces yielding this trace. (1 is one trace, 2 is two stacked traces, etc.)'],
'DataUse': [35,36,'Data use: '],
'DistanceCenter': [37,40,'Distance from center of the source point to the center of the receiver group (negative if opposite to direction in which line is shot).'],
'RecvGrpElev': [41,44,'Receiver group elevation (all elevations above the Vertical datum are positive and below are negative). The scalar in Trace Header bytes 69-70 applies to these values. The units are feet or meters as specified in Binary File Header bytes 3255-3256). The Vertical Datum should be defined through a Location Data stanza (see section D-1).'],
'SurfElevSrc': [45,48,'Surface elevation at source.'],
'SrcDepth': [49,52,'Source depth below surface (a positive number).'],
'ElevRcvGrp': [53,56,'Datum elevation at receiver group.'],
'ElevSrc': [57,60,'Datum elevation at source.'],
'WaterDepthSrc': [61,64,'Water depth at source.'],
'WaterDepthGrp': [65,68,'Water depth at group.'],
'ScalerElevDepth': [69,70,'Scalar to be applied to all elevations and depths specified in Trace Header bytes 41-68 to give the real value. Scalar = 1, +10, +100, +1000, or +10,000. If positive, scalar is used as a multiplier; if negative, scalar is used as a divisor.'],
'ScalerCoord': [71,72,'Scalar to be applied to all coordinates specified in Trace Header bytes 73-88 and to bytes Trace Header 181-188 to give the real value. Scalar = 1, +10, +100, +1000, or +10,000. If positive, scalar is used as a multiplier; if negative, scalar is used as divisor.'],
'X': [73,76,'Source coordinate - X. The coordinate reference system should be identified through an extended header Location Data stanza (see section D-1).If the coordinate units are in seconds of arc, decimal degrees or DMS, the X values represent longitude and the Y values latitude. A positive value designates east of Greenwich Meridian or north of the equator and a negative value designates south or west. '],
'Y': [77,80,'Source coordinate - Y.'],
'GrpX': [81,84,'Group coordinate - X.'],
'GrpY': [85,88,'Group coordinate - Y.'],
'CoordUnits': [89,90,'Coordinate units: '],
'WxVel': [91,92,'Weathering velocity. (ft/s or m/s as specified in Binary File Header bytes 3255-3256)'],
'SubWxVel': [93,94,'Subweathering velocity. (ft/s or m/s as specified in Binary File Header bytes 3255-3256)'],
'UpholeTSrc': [95,96,'Uphole time at source in milliseconds. Time in milliseconds as scaled by the scalar specified in Trace Header bytes 215-216.'],
'UpholeTGrp': [97,98,'Uphole time at group in milliseconds.'],
'SrcCorr': [99,100,'Source static correction in milliseconds.'],
'GrpCorr': [101,102,'Group static correction in milliseconds.'],
'TotStatic': [103,104,'Total static applied in milliseconds. (Zero if no static has been applied,)'],
'LagTA': [105,106,'Lag time A - Time in milliseconds between end of 240-byte trace identification header and time break. The value is positive if time break occurs after the end of header; negative if time break occurs before the end of header. Time break is defined as the initiation pulse that may be recorded on an auxiliary trace or as otherwise specified by the recording system.'],
'LatTB': [107,108,'Lag Time B - Time in milliseconds between time break and the initiation time of the energy source. May be positive or negative.'],
'Delay': [109,110,'Delay recording time - Time in milliseconds between initiation time of energy source and the time when recording of data samples begins. In SEG Y rev 0 this entry was intended for deep,water work if data recording does not start at zero time. The entry can be negative to accommodate negative start times (i.e. data recorded before time zero, presumably as a result of static application to the data trace). If a non,zero value (negative or positive) is recorded in this entry, a comment to that effect should appear in the Textual File Header.'],
'MuteStart': [111,112,'Mute time - Start time in milliseconds.'],
'MuteEnd': [113,114,'Mute time - End time in milliseconds.'],
'TraceSamples': [115,116,'Number of samples in this trace. Highly recommended for all types of data.'],
'SampleInterval': [117,118,'Sample interval in microseconds for this trace. The number of bytes in a trace record must be consistent with the number of samples written in the trace header. This is important for all recording media; but it is particularly crucial for the correct processing of SEG Y data in disk files (see Appendix C). If the fixed length trace flag in bytes 3503-3504 of the Binary File Header is set, the sample interval and number of samples in every trace in the SEG Y file must be the same as the values recorded in the Binary File Header. If the fixed length trace flag is not set, the sample interval and number of samples may vary from trace to trace.'],
'GainType': [119,120,'Gain type of field instruments: '],
'instrGainConst': [121,122,'Instrument gain constant (dB).'],
'instrInitGain': [123,124,'Instrument early or initial gain (dB).'],
'correlated': [125,126,'Correlated: 1 = no, 2 = yes'],
'SweepFreqStart': [127,128,'Sweep frequency at start (Hz).'],
'SweepFreqEnd': [129,130,'Sweep frequency at end (Hz).'],
'SweepLen': [131,132,'Sweep length in milliseconds.'],
'SweepType': [133,134,'Sweep type: '],
'SweepTaperLenStart': [135,136,'Sweep trace taper length at start in milliseconds.'],
'SweepTaperLenEnd': [137,138,'Sweep trace taper length at end in milliseconds.'],
'TaperType': [139,140,'Taper type'],
'AliasFiltFreq': [141,142,'Alias filter frequency (Hz), if used.'],
'AliasFiltSlope': [143,144,'Alias filter slope (dB/octave).'],
'NotchFiltFreq': [145,146,'Notch filter frequency (Hz), if used.'],
'NotchFiltSlope': [147,148,'Notch filter slope (dB/octave).'],
'LowCutFreq': [149,150,'Low-cut frequency (Hz), if used.'],
'HighCutFreq': [151,152,'High-cut frequency (Hz), if used.'],
'LowCutSlope': [153,154,'Low-cut slope (dB/octave)'],
'HighCutSlope': [155,156,'High-cut slope (dB/octave)'],
'Year': [157,158,'Year data recorded - The 1975 standard is unclear as to whether this should be recorded as a 2-digit or a 4-digit year and both have been used. For SEG Y revisions beyond rev 0, the year should be recorded as the complete 4-digit Gregorian calendar year (i.e. the year 2001 should be recorded as 200110 (7D116)).'],
'Day': [159,160,'Day of year (Julian day for GMT and UTC time basis).'],
'Hour': [161,162,'Hour of day (24 hour clock).'],
'Min': [163,164,'Minute of hour.'],
'Sec': [165,166,'Second of minute.'],
'TimeBasis': [167,168,'Time basis code: '],
'TraceWeight': [169,170,'Trace weighting factor - Defined as 2-N volts for the least signifi-cant bit. (N = 0, 1, ..., 32767)'],
'GeopGrpRlSw1': [171,172,'Geophone group number of roll switch position one.'],
'GeopGrpTr1': [173,174,'Geophone group number of trace number one within original field record.'],
'GeopGrpNoLstTr': [175,176,'Geophone group number of last trace within original field record.'],
'GapSize': [177,178,'Gap size (total number of groups dropped).'],
'OverTravel': [179,180,'Over travel associated with taper at beginning or end of line:'],
'Xensemble': [181,184,'X coordinate of ensemble (CDP) position of this trace (scalar in Trace Header bytes 71-72 applies). The coordinate reference system should be identified through an extended header Location Data stanza (see section D-1).'],
'Yensemble': [185,188,'Y coordinate of ensemble (CDP) position of this trace (scalar in bytes Trace Header 71-72 applies). The coordinate reference system should be identified through an extended header Location Data stanza (see section D-1).'],
'PostInLineNo': [189,192,'For 3-D poststack data, this field should be used for the in-line number. If one in-line per SEG Y file is being recorded, this value should be the same for all traces in the file and the same value will be recorded in bytes 3205-3208 of the Binary File Header.'],
'PostCrossLineNo': [193,196,'For 3-D poststack data, this field should be used for the cross-line number. This will typically be the same value as the ensemble (CDP) number in Trace Header bytes 21-24, but this does not have to be the case.'],
'Shotpoint': [197,200,'Shotpoint number - This is probably only applicable to 2-D poststack data. Note that it is assumed that the shotpoint number refers to the source location nearest to the ensemble (CDP) location for a particular trace. If this is not the case, there should be a comment in the Textual File Header explaining what the shotpoint number actually refers to.'],
'ShotpointScaler': [201,202,'Scalar to be applied to the shotpoint number in Trace Header bytes 197-200 to give the real value. If positive, scalar is used as a multiplier; if negative as a divisor; if zero the shotpoint number is not scaled (i.e. it is an integer. A typical value will be -10, allowing shotpoint numbers with one decimal digit to the right of the decimal point).'],
'TraceUnits': [203,204,'Trace value measurement unit:'],
'TransdConstMant': [205,208,'Transduction Constant - The multiplicative constant used to convert the Data Trace samples to the Transduction Units MANTISA'],
'TransdConstExp': [209,210,'Transduction Constant - The multiplicative constant used to convert the Data Trace samples to the Transduction Units EXPONENT'],
# 'TransdConst': [205,210,'Transduction Constant - The multiplicative constant used to convert the Data Trace samples to the Transduction Units (specified in Trace Header bytes 211-212). The constant is encoded as a four-byte, two\'s complement integer (bytes 205-208) which is the mantissa and a two-byte, two\'s complement integer (bytes 209-210) which is the power of ten exponent (i.e. Bytes 205-208 * 10**Bytes 209-210). '],
'TransdUnits': [211,212,'Transduction Units - The unit of measurement of the Data Trace samples after they have been multiplied by the Transduction Constant specified in Trace Header bytes 205-210.'],
'DevTrId': [213,214,'Device/Trace Identifier - The unit number or id number of the device associated with the Data Trace (i.e. 4368 for vibrator serial number 4368 or 20316 for gun 16 on string 3 on vessel 2). This field allows traces to be associated across trace ensembles independently of the trace number (Trace Header bytes 25-28).'],
'TimeScaler': [215,216,'Scalar to be applied to times specified in Trace Header bytes 95-114 to give the true time value in milliseconds. Scalar = 1, +10, +100, +1000, or +10,000. If positive, scalar is used as a multiplier; if negative, scalar is used as divisor. A value of zero is assumed to be a scalar value of 1.'],
'SrcTypeOrient': [217,218,'Source Type/Orientation - Defines the type and the orientation of the energy source. The terms vertical, cross-line and in-line refer to the three axes of an orthogonal coordinate system. The absolute azimuthal orientation of the coordinate system axes can be defined in the Bin Grid Definition Stanza (page 27).'],
# FIX: is this the right split for SrcEDir and SrcEDirTenths?
'SrcEDir': [219,222,'Source Energy Direction with respect to the source orientation - The positive orientation direction is defined in Bytes 217-218 of the Trace Header. The energy direction is encoded in tenths of degrees (i.e. 347.8 degrees is encoded as 3478).'],
'SrcEDirTenths': [223,224,'Source Energy Direction with respect to the source orientation - The positive orientation direction is defined in Bytes 217-218 of the Trace Header. The energy direction is encoded in tenths of degrees (i.e. 347.8 degrees is encoded as 3478).'],
'SrcMeasMantisa': [225,228,'Source Measurement - Describes the source effort used to generate the trace. The measurement can be simple, qualitative measurements such as the total weight of explosive used or the peak air gun pressure or the number of vibrators times the sweep duration. Although these simple measurements are acceptable, it is preferable to use true measurement units of energy or work.'],
'SrcMeasExp': [229,230,'The constant is encoded as a four-byte, two\'s complement integer (bytes 225-228) which is the mantissa and a two-byte, two\'s complement integer (bytes 209-230) which is the power of ten exponent (i.e. Bytes 225-228 * 10**Bytes 229-230).'],
'SrcMeasUnits': [231,232,'Source Measurement Unit - The unit used for the Source Measurement, Trace header bytes 225-230.'],
'Unassigned233': [233,240,'Unassigned - For optional information.']
}
##############################
# END OF TRACE DICTIONARIES
##############################
##############################
class SegyTraceHdr:
"""
Handle one trace in a segy file.
"""
def __init__(self, segyHdr, startOffset,mmapData=None):
"""
Load up a trace from an mmap'ed data file by jumping to the startOffset
"""
self.data = mmapData
self.start = startOffset
return
def size(self):
"""
Return the size of the whole trace datablock in the file (incl. trace header)
"""
assert(False) # FIX: nuke this function?
return
def getHeader(self,fieldName):
""" FIX: explain yourself! """
entry = traceHeader[fieldName]
start = entry[0] - 1 + self.start
end = entry[1] - 1 + self.start
length = entry[1]-entry[0] + 1
if length > 8:
print "get bin header... too large. skipping ",fieldName
return None
try:
val = struct.unpack(structIntCodes[length],self.data[start:end+1])[0]
except KeyError:
print "ERROR with key: ",key, " Length not found in table",length
assert(False)
return val
def printHeader(self):
""" Dump all header info to stdout"""
for key in traceHeader.keys():
val = self.getHeader(key)
print key, " = ",val
key = 'TraceId'; print key, ' := ', traceIdCode[self.getHeader(key)]
key = 'DataUse'; print key, ' := ', dataUse[self.getHeader(key)]
key = 'CoordUnits'; print key, ' := ',coordUnits[self.getHeader(key)]
key = 'GainType'; print key, ' := ',gainType[self.getHeader(key)]
key = 'SweepType'; print key, ' := ',sweepType[self.getHeader(key)]
key = 'TimeBasis'; print key, ' := ',timeBasis[self.getHeader(key)]
try:
key = 'OverTravel'; print key, ' := ',overTravel[self.getHeader(key)]
except:
print "UNKNOWN OVER TRAVEL!!!!!"
try:
key = 'TraceUnits'; print key, ' := ',traceUnits[self.getHeader(key)]
except:
print "UNKNOWN TRACE UNITS!!!!!"
try:
key = 'TransdUnits'; print key, ' := ',transdUnits[self.getHeader(key)]
except:
print "UNKOWN TransdUnits"
key = 'SrcTypeOrient'; print key, ' := ',srcTypeOrient[self.getHeader(key)]
key = 'SrcMeasUnits'; print key, ' := ',srcMeasUnits[self.getHeader(key)]
#key = ''; print key, ' := ',[self.getHeader(key)]
# FIX: Make a dictionary of the above, so it can be done automatically!
#traceFieldValuesEnglish = {
#
# }
######################################################################
# Test code
class TestSegy(unittest.TestCase):
"""
Run the segy class through its paces
Use one of the BPSIO xstar segy files. These are not proper segy
"""
def setUp(self):
#self.xstar = Segy("bpsio04L03.xstar")
self.xstar = SegyHdr(filename="bpsio01L00.xstar")
self.ol = SegyHdr(filename="9315MIG.sgy")
return
def tearDown(self):
#print "tearDown: in the end of the test case"
return
def testIntro(self):
"""Test the first few bytes of the file - magic number C1 or C 1 """
#print self.xstar.getTextHeaderLine(1)
#print self.xstar.getTextHeaderLine(2)
#print 'CLIENT: ',self.xstar.getTextHeaderField('CLIENT')
#print 'AREA: ',self.xstar.getTextHeaderField('AREA')
f = self.xstar.getBinEntry('SampleFormat')
#print f,"==",dataFormats[f]
#self.xstar.printBinaryHeaders()
#print "------------------------------"
#self.ol.printBinaryHeaders()
self.failUnless(self.xstar.hasR1TextFileHeader())
self.failUnless(not self.ol.hasR1TextFileHeader())
return
def testSegyTrace(self):
""" work on the trace class """
print "------------------------------\n TRACE 1 \n"
t = SegyTraceHdr(segyHdr=self.xstar,mmapData=self.xstar.data,startOffset=3600) # Load up the 1st trace... assume no ext text block
#t.printHeader()
key='TraceId'; print key,': ',traceIdCode[t.getHeader(key)]
# look! If we divide by 3600, we get degrees!
x = t.getHeader('X'); print x, x/3600.
y = t.getHeader('Y'); print y, y/3600.
# e.g.
# -432364 --> -120.101111111
# 123862 --> 34.4061111111
#print "------------------------------\n TRACE 1 OL \n"
#t = SegyTrace(self.ol.data,3600) # Load up the 1st trace... assume no ext text block
#t.printHeader()
return
def testSegyTrace2(self):
""" work on the trace class - load in two and make sure they are different"""
print "------------------------------\n TRACE 1 \n"
segy = Segy('bpsio01L00.xstar')
t = segy.getTraceHdr()
print "CHECK: tw getTraceHdr() start:", t.start
key='Shotpoint'; print key,': ',t.getHeader(key)
# look! If we divide by 3600, we get degrees!
x = t.getHeader('X'); print x, x/3600.
y = t.getHeader('Y'); print y, y/3600.
t2 = segy.getTraceHdr(1000)
print "CHECK: tw getTraceHdr(100) start:", t2.start
key='Shotpoint'; print key,': ',t2.getHeader(key)
# look! If we divide by 3600, we get degrees!
x2 = t2.getHeader('X'); print x, x/3600.
y2 = t2.getHeader('Y'); print y, y/3600.
self.failUnless ( x != x2)
self.failUnless ( y != y2)
def testTotalSegy(self):
print " ************************************"
segy = Segy('bpsio01L00.xstar')
#print segy.getTraceDataSize(), segy.getTraceSize()
#print segy.getTraceData()
print "Number of traces = ", segy.getNumberOfTraces()
segy.writeTraceToFile('1.dat',traceNumber=1)
segy.writeTraceToFile('2.dat',traceNumber=2)
segy.writeTraceToFile('3.dat',traceNumber=3)
segy.writeTraceToFile('100.dat',traceNumber=100)
########################################
# Run the test code
if __name__ == '__main__':
print "Testing segy"
#import doctest, segy
#doctest.testmod(segy)
unittest.main()
syntax highlighted by Code2HTML, v. 0.9.1