commit 333130f6ade684aebfd1f33fc31a1aac00945884 Author: cultpenguin Date: Thu Apr 14 05:25:07 2005 +0000 Initial revision diff --git a/segypy.py b/segypy.py new file mode 100644 index 0000000..35fbc36 --- /dev/null +++ b/segypy.py @@ -0,0 +1,207 @@ +""" +A python module for reading/writing/manipuating +SEG-Y formatted filed + +segy.segypy_version : The version of SegyPY +segy.segypy_verbose : A,punt of verbose information to the screen. +segy.getValue : Get a value from a binary string + +""" +# +# segypy : A Python module for reading and writing SEG-Y formatted data +# +# (C) Thomas Mejer Hansen, 2005 + +import struct + +# SOME GLOBAL PARAMETERS +segypy_version=0.1 +segypy_verbose=10; + +#endian='>' # Big Endian +#endian='<' # Little Endian +#endian='=' # Native + + +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') + +def readSegy(filename): + """ + getSegyHeader + """ + from pylab import * + data = open(filename).read() + + filesize=len(data) + vtxt = 'Length of data ; ',filesize + printverbose(vtxt,2) + + SH=getSegyHeader(filename) + + ntraces = (filesize-3600)/(SH['ns']*4+240) + printverbose(vtxt,2) + + vtxt = "readSegy : ntraces=",ntraces,"nsamples=",SH['ns'] + printverbose(vtxt,2) + + index=3600 + + + # ntraces=500; + Data = zeros((SH['ns'],ntraces)) + + printverbose("readSegy : reading data",2) + + for itrace in range(1,ntraces,1): + i=itrace +# print "Reading trace ",itrace," of ",ntraces + SegyTraceHeader,SegyTraceData=getSegyTrace(SH,itrace) + + for iss in range(1,SH['ns'],1): + Data[iss-1][itrace-1]=SegyTraceData[0][iss-1] + + printverbose("readSegy : read data",2) + + if (segypy_verbose>2): + imshow(Data) + title('pymat test') + grid(True) + show() + + + return Data,SH + +def getSegyTrace(SH,itrace): + data = open(SH["filename"]).read() + # GET TRACE HEADER + index=3200+(itrace-1)*(240+SH['ns']*4) + SegyTraceHeader=[]; + #print index + + # GET TRACE + index=3200+(itrace-1)*(240+SH['ns']*4)+240 + SegyTraceData = getValue(data,index,'float','>',SH['ns']) + return SegyTraceHeader,SegyTraceData + + +def getSegyHeader(filename,offset=3200): + + data = open(filename).read() + + printverbose('getSegyHeader : trying to read from '+filename,1) + + # START INDEX IN FILE + index=offset; + index=0; + + SegyHeader = {'filename': filename} + + TextualFileHeader,index = getValue(data,index,'c','>',3200); + + SegyHeader['Job'],index = getValue(data,index,'l','>'); + SegyHeader['Line'],index = getValue(data,index,'l','>') + SegyHeader['Reel'],index = getValue(data,index,'l','>') + SegyHeader['DataTracePerEnsemble'],index = getValue(data,index,'h','>') + SegyHeader['AuxiliaryTracePerEnsemble'],index = getValue(data,index,'short','>') + SegyHeader['dt'],index = getValue(data,index,'uint16'); + SegyHeader['dtOrig'],index = getValue(data,index,'uint16'); + SegyHeader['ns'],index = getValue(data,index,'uint16'); + SegyHeader['nsOrig'],index = getValue(data,index,'uint16'); + SegyHeader['DataSampleFormat'],index = getValue(data,index,'int16'); + SegyHeader['EnsembleFold'],index = getValue(data,index,'int16'); + SegyHeader['TraceSorting'],index = getValue(data,index,'int16'); + SegyHeader['VerticalSumCode'],index = getValue(data,index,'int16'); + + SegyHeader['SweepFrequencyStart'],index = getValue(data,index,'int16'); + SegyHeader['SweepFrequencyEnd'],index = getValue(data,index,'int16'); + SegyHeader['SweepLength'],index = getValue(data,index,'int16'); + SegyHeader['SweepType'],index = getValue(data,index,'int16'); + SegyHeader['SweepChannel'],index = getValue(data,index,'int16'); + SegyHeader['SweepTaperlengthStart'],index = getValue(data,index,'int16'); + SegyHeader['SweepTaperLengthEnd'],index = getValue(data,index,'int16'); + SegyHeader['TaperType'],index = getValue(data,index,'int16'); + SegyHeader['CorrelatedDataTraces'],index = getValue(data,index,'int16'); + SegyHeader['BinaryGain'],index = getValue(data,index,'int16'); + SegyHeader['AmplitudeRecoveryMethod'],index = getValue(data,index,'int16'); + SegyHeader['MeasurementSystem'],index = getValue(data,index,'int16'); + SegyHeader['ImpulseSignalPolarity'],index = getValue(data,index,'int16'); + SegyHeader['VibratoryPolarityCode'],index = getValue(data,index,'int16'); + + index=3500; + SegyHeader['SegyFormatRevisionNumber'],index= getValue(data,index,'uint16'); + SegyHeader['FixedLengthTraceFlag'],index=getValue(data,index,'uint16'); + SegyHeader['NumberOfExtTextualHeaders'],index= getValue(data,index,'uint16'); + + printverbose('getSegyHeader : succesfully read '+filename,1) + + return SegyHeader + + + +def getValue(data,index,ctype='l',endian='>',number=1): + """ + GET + """ + if (ctype=='l')|(ctype=='long'): + size=l_long + ctype='l' + elif (ctype=='L')|(ctype=='ulong'): + 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' + else: + printverbose('Bad Ctype : ' +ctype,-1) + + cformat=endian + ctype*number + + printverbose('cformat : ' + cformat,20) + + index_end=index+size*number + HeaderValue=struct.unpack(cformat, data[index:index_end]) + + str = 'getSegyHeaderValue','start=',index,' size=',size, 'number=',number,'Value=',HeaderValue,'cformat=',cformat + printverbose(str,20) + + if number==1: + return HeaderValue[0], index_end + else: + return HeaderValue,index_end + + +def version(): + return segypy_version + +def print_version(): + print 'SegyPY version is ', segypy_version + +def printverbose(txt,level): + if level') +#Job,size = segypy.getValue(data,index,'long','>') +#index=index+size +#Line,size = segypy.getValue(data,index,'l','>') +#index =index+size +#Reel,size = segypy.getValue(data,index,'l','>') +#index=index+size +#DataTracePerEnsemble, size = segypy.getValue(data,index,'h','>') +#DataTracePerEnsemble, size = segypy.getValue(data,index,'short','>') + +print "GET SEGY GEADER" +SH=segypy.getSegyHeader(filename) + +print SH +#print SH['Reel'] + +print "TESTING READING OF SEGY FILE" +segypy.readSegy(filename); + +