readsegy.m

matlab编写的地震数据格式转换

    SegymatVerbose([mfilename,' - Using supplied SegyHeader'])
  end
  
  % APPLY CHANGES TO SEGY HEADER IF NEEDE
  if isempty(revision)==0, 
    SegyHeader.SegyFormatRevisionNumber=revision; 
    SegymatVerbose([mfilename,' - Manually set SEG-Y revision to ',num2str(revision)])
  end
  if isempty(dsf)==0, 
    SegyHeader.DataSampleFormat=dsf; 
  end
  
  % JUST SOME INFORMATION TO WRITE TO SCREEN :
  
  Revision=SegyHeader.SegyFormatRevisionNumber;
  if Revision>0, Revision=1; end
  
  if (SegyHeader.DataSampleFormat>length(SegyHeader.Rev(Revision+1).DataSampleFormat));
    SegymatVerbose([mfilename,' : WARNING : YOU HAVE SELECTED (OR THE FILE IS FORMATTED SUCH THAT) A DATASAMPLE FORMAT THAT IS NOT DEFINED. \nREMEBER IEEE IS NOT SPECIFIED IN THE SEGY REV0 STANDARD !'])
    
    if (Revision==0)
      SegymatVerbose([mfilename,' : TRYING TO USE REVISION 1 AS OPPOSED TO REVISION 0'])
      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=(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=1;
  if DataEnd>1e+6, dwaitbar=10;  end
  if DataEnd>1e+7, dwaitbar=50;  end
  if DataEnd>1e+8, dwaitbar=200;  end
  traceinfile=0;  
  outtrace=0;
  tic;
  toc_old=toc;
  if doWaitBar==1;
    hw=waitbar(0,['Reading Segy - ',txt]);
  end
  while (~(ftell(segyid)>=DataEnd))
    traceinfile=traceinfile+1;
    
    usetrace=1; % DEFAULT USING TRACE WHEN [1].
    
    ishow=500;
    if ((traceinfile/ishow)==round(traceinfile/ishow)), 
      PerTrace=(toc-toc_old)/ishow;
      TimeLeft=(ntraces-traceinfile)*PerTrace;
      txt=sprintf('Reading trace %d/%d, (%5.0fs left)',traceinfile,ntraces,TimeLeft);
      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
        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);
      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
  if doWaitBar==1
    close(hw);
  end
  SegymatVerbose([mfilename,' : Elapsed time ',num2str(toc)]);
  t=outtrace;
  
  % 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
  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