#!/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