#!/usr/bin/env python
print "WARNING - this code is not yet functional!!!"
"""
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.4 2005/05/25 16:18:43 schwehr Exp schwehr $
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
# 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"
}
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']
}
# 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' }
class Segy:
""" Base class for reading segy seismic files
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
return
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
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]
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 getBinaryFileEntryHeader(self,fieldName):
""" FIX: explain yourself! """
entry = binaryHeaderEntries[fieldName]
#print fieldName, entry
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
#print
#print fieldName," ",start," ",length, structIntCodes[length], " LEN==",len(self.data[start:end+1])
val = struct.unpack(structIntCodes[length],self.data[start:end+1])[0]
#print "val =",val
return val
def printBinaryHeaders(self):
""" Dump all header into to stdout"""
for key in binaryHeaderEntries.keys():
val = self.getBinaryFileEntryHeader(key)
print key, " = ",val
######################################################################
# Trace parsing
######################################################################
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.'],
}
class SegyTrace:
"""
Handle one trace in a segy file.
"""
def __init__(self, mmapData, startOffset):
"""
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)
"""
return None # FIX
def getHeader(self,fieldName):
""" FIX: explain yourself! """
entry = traceHeader[fieldName]
print
print fieldName, entry[0],entry[1]
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
print fieldName," ",start," ",length, structIntCodes[length], " LEN =",len(self.data[start:end+1])
print self.data[start:end+1]
val = struct.unpack(structIntCodes[length],self.data[start:end+1])[0]
#print "val =",val
return val
def printHeader(self):
""" Dump all header info to stdout"""
for key in traceHeader.keys():
val = self.getHeader(key)
print key, " = ",val
######################################################################
# 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 = Segy("bpsio01L00.xstar")
self.ol = Segy("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.getBinaryFileEntryHeader('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 = SegyTrace(self.xstar.data,3600) # Load up the 1st trace... assume no ext text block
t.printHeader()
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
########################################
# 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