mirror of
https://github.com/wassname/segpy.git
synced 2026-07-30 12:31:11 +08:00
Initial revision
This commit is contained in:
cultpenguin
@@ -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<segypy_verbose:
|
||||
print 'SegyPY',segypy_version,': ',txt
|
||||
|
||||
|
||||
|
||||
#def read(filename):
|
||||
# read_segy(filename)
|
||||
#
|
||||
#def write(filename):
|
||||
# write_segy(filename)
|
||||
+42
@@ -0,0 +1,42 @@
|
||||
#
|
||||
#
|
||||
# testsegy.py
|
||||
#
|
||||
import struct
|
||||
import segypy
|
||||
|
||||
filename='mini.sgy'
|
||||
|
||||
data = open(filename).read()
|
||||
|
||||
l_long = struct.calcsize('l')
|
||||
l_ulong = struct.calcsize('L')
|
||||
l_short = struct.calcsize('h')
|
||||
l_ushort = struct.calcsize('H')
|
||||
|
||||
segypy_version=segypy.version()
|
||||
|
||||
print 'filename=',filename
|
||||
|
||||
|
||||
#index=3200;
|
||||
#Job,size = segypy.getValue(data,index,'l','>')
|
||||
#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);
|
||||
|
||||
|
||||
Reference in New Issue
Block a user