readsegy.m

这是用matlab对segy数据进行处理

        Revision=1;

        if (SegyHeader.DataSampleFormat>length(SegyHeader.Rev(Revision+1).DataSampleFormat));
            SegymatVerbose([mfilename,' : FATAL ERROR : STILL THE DATASAMPLE FORMAT IS NOT SUPPRTED - EXITING (Report error to [e-mail])'])
        else
            SegymatVerbose([mfilename,' : APPARENT SUCCES CHANING FROM Revision 0 to 1 - Continuing'])
            SegyHeader.SegyFormatRevisionNumber=1; % FORCING REVISION TO BE 1 !!!
        end

    end

end


FormatName=SegyHeader.Rev(Revision+1).DataSampleFormat(SegyHeader.DataSampleFormat).name;
Format=SegyHeader.Rev(Revision+1).DataSampleFormat(SegyHeader.DataSampleFormat).format;
BPS=SegyHeader.Rev(Revision+1).DataSampleFormat(SegyHeader.DataSampleFormat).bps;
txt=['SegyRevision ',sprintf('%0.4g',Revision),', ',FormatName,'(',num2str(SegyHeader.DataSampleFormat),')'];


ns=SegyHeader.ns;
% YOU CAN FORCE FixedLengthTraceFlag=1;
% This will make the code much faster (especially when using
% the 'jump' option) but reading data with varying trace lengt will fail.
% It is here since many old data sets with Constant trace length
% has FixedLengthTraceFlag=0;
%
% As of version 1.01 this has been enable by default.
% Change the variable below to '0' if you do not want this behaviour
%
SegyHeader.FixedLengthTraceFlag=1;


SegymatVerbose([mfilename,' : Reading Data'],90);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% READ DATA
%Segy=fread(segyid,4000,'float32');
fseek(segyid,0,'eof'); DataEnd=ftell(segyid);
fseek(segyid,0,'eof'); DataEnd=ftell(segyid);

DataStart=3600+3200*SegyHeader.NumberOfExtTextualHeaders;
fseek(segyid,DataStart,'bof');       % Go to the beginning of the file
fseek(segyid,DataStart,'bof');       % Go to the beginning of the file


ntraces=round((DataEnd-DataStart)./(240+(SegyHeader.ns)*(BPS/8)));

SegymatVerbose(['Number of Samples Per Trace=',num2str(SegyHeader.ns)])
SegymatVerbose(['Number of Traces=',num2str(ntraces)])
if (ntraces~=round(ntraces))
    SegymatVerbose(['Trace lengths seems to vary. trying to read the file anyway'])
end

existJump=~isempty(jump);
existHeader=~isempty(header);
existTmin=~isempty(tmin);
existTmax=~isempty(tmax);

if existJump==1,
    out_ntraces=ceil(ntraces/jump);
else
    out_ntraces=ntraces;
end



dwaitbar=0;
if DataEnd./jump>1e+6, dwaitbar=10;  end
if DataEnd./jump>1e+7, dwaitbar=50;  end
if DataEnd./jump>1e+8, dwaitbar=200;  end
traceinfile=0;
outtrace=0;
tic;
toc_old=toc;
if doWaitBar==1;
    hw=waitbar(0,['Reading Segy - ',txt]);
end


% LOOP OVER TRACES

t0=now;
tlast=t0;
pos0=ftell(segyid);

while (~(ftell(segyid)>=DataEnd))
   
    traceinfile=traceinfile+1;

    usetrace=1; % DEFAULT USING TRACE WHEN [1].

    ishow=10000;
    itime=1/(24*3600)*2; % Min time between verbose info to screen  
    if (((traceinfile/ishow)==round(traceinfile/ishow))&((now-tlast)>itime)),

        
        tnow=now;
        tlast=tnow;
        posnow=ftell(segyid);
        tend=t0+DataEnd.*(tnow-t0)./(posnow-pos0);
        tleft=(tend-tnow)*3600*24;
        txt=sprintf('Reading trace %d/%d, (%5.0fs left) (est end %s)',traceinfile,ntraces,tleft,datestr(tend));
        toc_old=toc;
        SegymatVerbose(txt)
    end
    TraceStart=ftell(segyid);

    % IF 'JUMP' IS SET THEN CHECK IF THIS TRACE SHOULD BE SKIPPED
    if existJump==1
        if (traceinfile/jump)~=round(traceinfile/jump), usetrace=0; end
    end

    if ((usetrace==0)&(SegyHeader.FixedLengthTraceFlag==1)),
        % SKIP FORWARD IN FILE'
        skip=240+(BPS/8)*SegyHeader.ns;
        fseek(segyid,skip,'cof');
        %SegymatVerbose([num2str(traceinfile),' - SKIPPING TRACE ... ',num2str(outtrace)])
    else
        SingleSegyTraceHeaders=GetSegyTraceHeader(segyid,TraceStart,Format,SegyHeader.ns,[]);
        SingleSegyData.data=GetSegyTraceData(segyid,SingleSegyTraceHeaders.ns,SegyHeader);

        if SingleSegyTraceHeaders.TraceNumber<1
            SingleSegyTraceHeaders.TraceNumber=traceinfile;
            SegymatVerbose(sprintf('TraceNumber malformatetd. Setting TraceNumber=%d',traceinfile),10);
        end

        SegymatVerbose(sprintf('ns=%d, Trace in line : %d, Trace in file : %d, ns=%10.5f dt=%10.5f',SingleSegyTraceHeaders.ns,SingleSegyTraceHeaders.TraceSequenceLine,SingleSegyTraceHeaders.TraceSequenceFile,SingleSegyTraceHeaders.ns,SingleSegyTraceHeaders.dt),10)


    end


    % IF HEADER MIN MAX HAS BEEN CHOSEN, THEN CHECK THAT TRACE IS GOOD ENOUGH
    if ((existHeader==1)&(usetrace==1))
        headervalue=getfield(SingleSegyTraceHeaders,header);
        if ((headervalue<headermin)|(headervalue>headermax))
            usetrace=0;
        end
    end

    % USE THIS TRACE IF usetrace=1 !!
    if usetrace==1,
        %% IF TIME RANGE IS SPECIFIED, THEN EXTRACT THIS
        if (existTmin==1)&(existTmax==1)
            % NEXT LINE SHOULD CONSIDER THAT ns in Trace and Segy Header could vary !!!
            origtrange=[1:1:SegyHeader.ns].*SegyHeader.dt.*1e-6-SingleSegyTraceHeaders.DelayRecordingTime.*1e-3;
            gooddata=find(origtrange>tmin & origtrange<tmax);
            SingleSegyData.data=SingleSegyData.data(gooddata);
            % CHECK NEXT LINE TAHT DelatRec... is in micro seconds
            SingleSegyTraceHeaders.DelayRecordingTime=tmin;
            SingleSegyTraceHeaders.ns=length(gooddata);
            ns=length(gooddata); %  for use below
        end

        outtrace=outtrace+1;

        if (outtrace==1),
            % Preallocate RAM
            ta1=now;
            SegymatVerbose(sprintf('Pre allocating RAM ntraces=%d out_traces=%d',ntraces,out_ntraces));
            SegyData=repmat(struct('data',zeros(ns,1)),1,out_ntraces);
            SegyTraceHeaders=repmat(SingleSegyTraceHeaders,1,out_ntraces);           
            %whos SegyData SegyTraceHeaders
            %save T1
            ta2=now;
            t0=t0+ta2-ta1;
        end

        
        SegyData(outtrace).data=SingleSegyData.data;
        SegyTraceHeaders(outtrace)=SingleSegyTraceHeaders;
    
        
        if doWaitBar==1,
            if ((outtrace/dwaitbar)==round(outtrace/dwaitbar))
                waitbar(ftell(segyid)/DataEnd,hw);
            end
        end

    end

end
%save T2
%whos SegyData SegyTraceHeaders

if outtrace==0
    SegymatVerbose(sprintf('%s : No traces read!',mfilename));
    SegyTraceHeaders=[];
    Data=[];
    return
end

if doWaitBar==1
    try
        close(hw);
    end
end
SegymatVerbose([mfilename,' : Elapsed time ',num2str(toc),'  ended at ',datestr(now)]);
t=outtrace;

% Write time to SegyHeader
SegyHeader.ns=ns;
%  SegyHeader.time=[1:1:SegyHeader.ns].*SegyHeader.dt./1e+6+SegyTraceHeaders(1).DelayRecordingTime./1e+3;
SegyHeader.time=[1:1:SegyHeader.ns].*SegyHeader.dt./1e+6+SegyTraceHeaders(1).DelayRecordingTime;

% CHANGE SEGY HEADER IF TIME RANGE WAS SET
%if (existTmin==1)&(existTmax==1)
%  SegyHeader.ns=ns;
%  SegyHeader.time=[1:1:SegyHeader.ns].*SegyHeader.dt./1e+6+SegyTraceHeaders(1).DelayRecordingTime;
%end

% Make sure that only read SegyTraceHEaders are returned
if outtrace~=out_ntraces
    SegyTraceHeaders=SegyTraceHeaders(1:outtrace);
end

% MOVE DATA from SegyData.data to a regular variable
% THIS STEP COULD BE AVOIDED WHEN FixedTraceLength=1, which is almost
% allways the case... Should speed up things and reduce memory reqs.

if SkipData==1,
    Data=[];
else

    try
        Data=[SegyData(1:outtrace).data];
    catch


        Data=zeros(ns,outtrace);
        for i=1:outtrace
            try
                Data(:,i)=SegyData(i).data;
            catch
                errmsg=lasterr;
                if isempty(SegyData(i).data)
                    errmsg='Empty data in trace';
                elseif (strfind(errmsg, 'In an assignment  A(:,matrix) = B, the number of rows in A and B'))
                    nns=length(SegyData(i).data);
                    if nns<ns
                        errmsg='Length of trace varies - padding with zeros';
                        Data(1:nns,i)=SegyData(i).data;
                    else
                        errmsg='Length of trace varies - truncating';
                        Data(:,i)=SegyData(i).data(1:ns);
                    end
                end
                SegymatVerbose(sprintf('Had a problem at trace %d : %s',i,errmsg))
            end
        end
    end
end