writesegy.m

matlab源程序

% WriteSegy : writes data to disk using SEGY REV 1 standard.
%
% EX
% WriteSegy('datacube.segy',data,'dt',.004,'Inline3D',Inline,'Crossline3D',Crossline,'cdpX',X,'cdpY',Y);
%
% to use a specific SEG revision use :
% WriteSegy('test.segy',seisdata,'revision',0); % SEG-Y Revision 0
% WriteSegy('test.segy',seisdata,'revision',1); % SEG-Y Revision 1
%
% to use a specific Data Sampling Format use :
% WriteSegy('test.segy',seisdata,'dsf',1); % IBM FLAOTING POINT
%
% Forice Revision 1 and IEEE Floating point :
% WriteSegy('test.segy',seisdata,'dsf',5,'revision',1);
%
% See also : WriteSegyStructure, WriteSu, WriteSuStructure
%

%
% (C) 2001-2004, Thomas Mejer Hansen, [e-mail]/[e-mail]
%
%    This program is free software; you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation; either version 2 of the License, or
%    (at your option) any later version.
%
%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program; if not, write to the Free Software
%    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
%
% Modifications made by Urs Boeniger allowing to add all possible headers
% directly within this function.

function WriteSegy(filename,data,varargin);

    if nargin<2;
        help WriteSegy;
        return;
    end

    [ns,ntraces]=size(data);
    
    if nargin>2
        if length(varargin)>2;
            varargin=(reshape(varargin,2,(length(varargin)/2)))';
        end
        fields=varargin(:,1);
        headers=cell2struct(varargin(:,2),fields,1);
        if isfield(headers,'dt')~=1
            headers.dt=0.004.*1e+6;
            SegymatVerbose([mfilename,' - No dt set. Using dt=',num2str(headers.dt)],0)
        else
            headers.dt=headers.dt.*1e+6;
        end
    else
        headers.dt=.004.*1e+6;        
        SegymatVerbose([mfilename,' - No dt set. Using dt=',num2str(headers.dt)],0)
    end
    %%%%%%%%%%
    % SET UP SegyHeader structure.
    if isfield(headers,'TextualFileHeader')==1, SegyHeader.TextualFileHeader=headers.TextualFileHeader;end
    if isfield(headers,'Job')==1, SegyHeader.Job=headers.Job;end;
    if isfield(headers,'Line')==1, SegyHeader.Line=headers.Line;end
    if isfield(headers,'Reel')==1, SegyHeader.Reel=headers.Reel;end
    if isfield(headers,'DataTracePerEnsemble')==1, SegyHeader.DataTracePerEnsemble=headers.DataTracePerEnsemble;end
    if isfield(headers,'AuxiliaryTracePerEnsemble')==1, SegyHeader.AuxiliaryTracePerEnsemble=headers.AuxiliaryTracePerEnsemble;end
    if isfield(headers,'dt')==1, SegyHeader.dt=headers.dt;end
    if isfield(headers,'dtOrig')==1, SegyHeader.dtOrig=headers.dtOrig;end
    if exist('ns')==1, SegyHeader.ns=ns; end
    if isfield(headers,'ns')==1, SegyHeader.ns=headers.ns;end
    if isfield(headers,'nsOrig')==1, SegyHeader.nsOrig=headers.nsOrig;end
    if isfield(headers,'EnsembleFold')==1, SegyHeader.EnsembleFold=headers.EnsembleFold;end
    if isfield(headers,'TraceSorting')==1, SegyHeader.TraceSorting=headers.TraceSorting;end
    if isfield(headers,'VerticalSumCode')==1, SegyHeader.VerticalSumCode=headers.VerticalSumCode;end
    
    %%%%%%%%%%
    %% CHOOSE WHICH SEGY REVISION AND DATATYPE TO USE
    %% DEFAULT IS REVISION 1, 4-byte IEEE
    
    % IF A SPECFIC REVISION HAS BEEN CHOSEN, USE THAT
    if isfield(headers,'revision')==1,
        if headers.revision==1,
            SegyHeader.SegyFormatRevisionNumber=100;
        else
            SegyHeader.SegyFormatRevisionNumber=0;
        end
        SegymatVerbose([mfilename,' :  Using user specified SEG Y revision : ',num2str(headers.revision)],1)
    end
    
    
    % IF A SPECFIC DATA SAMPLING FORMAT HAS BEEN SELECTED USE THAT
    if isfield(headers,'dsf')==1,
        SegyHeader.DataSampleFormat=dsf;
        SegymatVerbose([mfilename,' :  Using user specified Data Sample Format : ',num2str(revision)],1)
    end
    
    % UNCOMMENT THE FOLLOWING TWO LINES TO USE REVISION 1 (2002)
    %SegyHeader.SegyFormatRevisionNumber=100; % 2002 SEG Y STYLE
    %SegyHeader.DataSampleFormat=2; % '1'->4-byte IBM floating point
    %                                 '2'->4-byte two's complement integer
    %                                 '3'->2-byte two's complement integer
    %                                 '5'->4-byte IEEE floating point (default)
    %                                 '8'->1-byte two's complement integer
    
    % UNCOMMENT THE FOLLOWING TWO LINES TO USE REVISION 0 (1975)
    %SegyHeader.SegyFormatRevisionNumber=0; % 1975 SEoG Y STYLE
    %SegyHeader.DataSampleFormat=1; % '1'->4-byte IBM Floating Point
    
    % OPEN SEGY FILE HANDLE
    segyid = fopen(filename,'w','b');   % ALL DISK FILES ARE IN BIG
                                        % ENDIAN FORMAT, ACCORDING SEG
                                        % Y rev 1

    % Write segy header
    SegyHeader=PutSegyHeader(segyid,SegyHeader);
    
    % Separate between single and multivalued header entries
    fnames=fieldnames(headers);
    fsize=zeros(size(fnames));
    for k=1:length(fnames)
        fsize(k)=eval(['length(',['headers.',cell2mat(fnames(k))],')']);
    end
    mvals=1:length(fnames);
    svals=find(fsize==1);mvals(svals)=[];

    for i=1:ntraces;
        SegyTraceHeader=headers;
        if (i/100)==round(i/100),
            SegymatVerbose(['writing trace ',num2str(i),' of ',num2str(ntraces)],0);
        end
        % Initialize SegyTraceHeader structure
        if isfield(headers,'headers')==0;
            SegyTraceHeader=InitSegyTraceHeader(ns,0.004.*1e+6);
        end
        if isfield(headers,'TraceNumber')==0
            SegyTraceHeader.TraceNumber=i;
            SegyTraceHeader.TraceSequenceFile=i;
        end
        % ADD additional header entries
        % Write single valued trace headers into SegyTraceHeaders struct
        if isempty(svals)~=1;
            for t=1:length(svals);
                evalc(['SegyTraceHeader.',cell2mat(fnames(svals(t))),'= headers.',cell2mat(fnames(svals(t)))]);
            end
        end
        % Write multi valued trace headers into SegyTraceHeaders struct running
        % over the indices i.
        if isempty(mvals)~=1;
            for tt=1:length(mvals);
                evalc(['SegyTraceHeader.',cell2mat(fnames(mvals(tt))),'= headers.',cell2mat(fnames(mvals(tt))),'(',num2str(i),')']);
            end
        end
        PutSegyTrace(segyid,data(:,i),SegyTraceHeader,SegyHeader);
    end
    
    % CLOSE SEGY FILE HANDLE
    fclose(segyid);