diff --git a/segypy2.py b/segypy2.py
new file mode 100644
index 0000000..26f1041
--- /dev/null
+++ b/segypy2.py
@@ -0,0 +1,1047 @@
+"""
+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 pickle
+import struct
+import logging
+
+import numpy
+from numpy import transpose
+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
+STH_def["ReceiverGroupElevation"]={"pos":40 ,"type":"int32"} #'int32');             %40
+STH_def["SourceSurfaceElevation"]={"pos":44 ,"type":"int32"} #'int32');             %44
+STH_def["SourceDepth"]={"pos":48 ,"type":"int32"} #'int32');             %48
+STH_def["ReceiverDatumElevation"]={"pos":52 ,"type":"int32"} #'int32');             %52
+STH_def["SourceDatumElevation"]={"pos":56 ,"type":"int32"} #'int32');             %56
+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
+STH_def["TotalStaticApplied"]={"pos":102 ,"type":"int16"} #'int16');  %102
+STH_def["LagTimeA"]={"pos":104 ,"type":"int16"} #'int16');  %104
+STH_def["LagTimeB"]={"pos":106 ,"type":"int16"} #'int16');  %106
+STH_def["DelayRecordingTime"]={"pos":108 ,"type":"int16"} #'int16');  %108
+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]
+
+STH_def["SweepFrequenceStart"]={"pos":126 ,"type":"int16"} #'int16');  %126
+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 imageSegy(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 getDefaultSegyHeader(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 (tmpkey.has_key('def')):
+            val=tmpkey['def']
+        else:
+            val=0
+        SH[key]=val
+
+    SH["ntraces"]=ntraces;
+    SH["ns"]=ns;
+    
+        
+    return SH
+def getDefaultSegyTraceHeaders(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 (tmpkey.has_key('def')):
+            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 getSegyTraceHeader(SH,THN='cdp',data=None):
+    """
+    getSegyTraceHeader(SH,TraceHeaderName)
+    """
+    logging.info("get trace header {0}".format(THN))
+    bps=getBytePerSample(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 = createValueGetter(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 getAllSegyTraceHeaders(SH,data='none'):
+    SegyTraceHeaders = {'filename': SH["filename"]}
+
+    printverbose('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 = getSegyTraceHeader(SH,key,data)      
+        SegyTraceHeaders[key]=sth      
+
+        txt =  "getAllSegyTraceHeaders :  " + key      
+        printverbose(txt,logging.DEBUG)
+
+    return SegyTraceHeaders
+    
+
+
+def readSegy(filename)    :
+    """
+    Data,SegyHeader,SegyTraceHeaders=getSegyHeader(filename)
+    """
+    
+    printverbose("readSegy : Trying to read "+filename, logging.DEBUG)
+
+    with open(filename, 'rb') as raw:
+        data = raw.read()
+
+    filesize = os.path.getsize(filename)
+
+    segy_header=getSegyHeader(filename)
+
+    bytes_per_sample = getBytePerSample(segy_header)
+
+    num_traces = (filesize - REEL_HEADER_NUM_BYTES) / (segy_header['ns'] * bytes_per_sample + TRACE_HEADER_NUM_BYTES)
+
+    printverbose("readSegy : Length of data : " + str(filesize), logging.DEBUG)
+
+    segy_header["ntraces"] = num_traces
+
+    ndummy_samples = TRACE_HEADER_NUM_BYTES / bytes_per_sample
+    printverbose("readSegy : ndummy_samples="+str(ndummy_samples), logging.DEBUG)
+    printverbose("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 = getSegyTraceHeader(segy_header, "Inline3D", data)
+    crossline_3d_trace_header = getSegyTraceHeader(segy_header, "Crossline3D", data)
+    num_samples_trace_header = getSegyTraceHeader(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")
+    
+    
+    printverbose("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]
+
+    printverbose("readSegy : SEG-Y revision = "+str(revision), logging.INFO)
+    printverbose("readSegy : DataSampleFormat="+str(dsf)+"("+DataDescr+")", logging.INFO)
+
+    if (segy_header["DataSampleFormat"]==1):
+        printverbose("readSegy : Assuming DSF=1, IBM FLOATS", logging.INFO)
+        pseudo_samples = getValue(data,index,'ibm',endian,num_pseudo_samples)
+    elif (segy_header["DataSampleFormat"]==2):
+        printverbose("readSegy : Assuming DSF=" + str(segy_header["DataSampleFormat"]) + ", 32bit INT", logging.INFO)
+        pseudo_samples = getValue(data,index,'l',endian,num_pseudo_samples)
+    elif (segy_header["DataSampleFormat"]==3):
+        printverbose("readSegy : Assuming DSF=" + str(segy_header["DataSampleFormat"]) + ", 16bit INT", logging.INFO)
+        pseudo_samples = getValue(data,index,'h',endian,num_pseudo_samples)
+    elif (segy_header["DataSampleFormat"]==5):
+        printverbose("readSegy : Assuming DSF=" + str(segy_header["DataSampleFormat"]) + ", IEEE", logging.INFO)
+        pseudo_samples = getValue(data,index,'float',endian,num_pseudo_samples)
+    elif (segy_header["DataSampleFormat"]==8):
+        printverbose("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:
+        printverbose("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
+
+    printverbose("readSegy :  read data", logging.DEBUG)
+    SegyTraceHeaders = None # Temp
+    return signed_data_3d, segy_header, SegyTraceHeaders
+
+
+def getSegyTrace(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=getBytePerSample(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 = getValue(data,index,'float',endian,SH['ns'])
+    return SegyTraceHeader,SegyTraceData
+
+def getSegyHeader(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 = getValue(data,pos,format,endian);    
+
+        txt =  str(pos) + " " + str(format) + "  Reading " + key +"="+str(SegyHeader[key])
+        printverbose(txt, logging.INFO)
+    
+    # SET NUMBER OF BYTES PER DATA SAMPLE
+    bps=getBytePerSample(SegyHeader)
+
+    filesize=len(data)
+    ntraces = (filesize-3600)/(SegyHeader['ns']*bps+240)
+    SegyHeader["ntraces"]=ntraces;
+
+    printverbose('getSegyHeader : succesfully read '+filename, logging.INFO)
+
+    return SegyHeader
+
+def writeSegy(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
+    
+    """
+
+    printverbose("writeSegy : Trying to write "+filename, logging.DEBUG)
+
+    N=Data.shape
+    ns=N[0]
+    ntraces=N[1]
+    print ntraces,ns
+
+    SH = getDefaultSegyHeader(ntraces,ns);
+    STH = getDefaultSegyTraceHeaders(ntraces,ns,dt)
+
+    # ADD STHin, if exists...
+    for key in STHin.keys():
+        print key
+        for a in range(ntraces):
+            STH[key]=STHin[key][a]
+    
+
+    writeSegyStructure(filename,Data,SH,STH)
+    
+
+def writeSegyStructure(filename,Data,SH,STH):
+    """
+    writeSegyHeader(filename,Data,SegyHeader,SegyTraceHeaders)
+
+    Write SEGY file using SegyPy data structures
+
+    See also readSegy
+
+    (c) 2005, Thomas Mejer Hansen
+    
+    """
+
+    printverbose("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]
+
+    printverbose("writeSegyStructure : SEG-Y revision = "+str(revision), logging.INFO)
+    printverbose("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);    
+        putValue(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
+        printverbose('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)
+            printverbose(txt, logging.INFO)
+            putValue(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 putValue(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
+    else:
+        printverbose('Bad Ctype : ' +ctype, logging.ERROR)
+
+    cformat=endian + ctype*number
+
+    printverbose('putValue : cformat :  ' + cformat + ' ctype=' + ctype, logging.DEBUG)
+    
+    strVal=struct.pack(cformat, value)
+    fileid.seek(index)
+    fileid.write(strVal);
+
+
+    return 1
+
+def createValueGetter(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:
+        printverbose('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'):
+            printverbose('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 getValue(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:
+        printverbose('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'):
+        printverbose('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 printverbose(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 getBytePerSample(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))
+
+    printverbose("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 InlineX(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 = readSegy(r'C:\Users\rjs\opendtectroot\Blake_Ridge_Hydrates_3D\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)