write_segy_file.m

基于Matlab的地震数据处理显示和测井数据显示于处理的小程序

function write_segy_file(seismic,filename,varargin)
% Function writes seismic data to disk in SEG-Y file format ; if the start 
% time is greater than zero, zeros are prepended. If it is less than zero, a 
% warning message is printed. The only floating-point format supported is 
% IEEE big-endian, which is one of the official SEG-Y standard formats. ProMAX
% recognizes the numeric code associated with this format and reads SeisLab-
% generated SEG-Y files without need for for any special settings.
%
% Written by: E. R.: March 12, 2000
% Last updated: July 23, 2006: Replaced header mnemonics "field_rec_no" by "ffid" 
%                              and  "trc_id" by "trc_type";
%                              write field "sort" with trace-sort code (if it 
%                              exists) to bytes 29-30 of the binary tape header
%
%        write_segy_file(seismic,filename,varargin)
% INPUT
% seismic     structure; the following fields are required
%        traces    seismic traces
%        first
%        last
%        step
%        units
% filename    Full filename; if empty (or invalid) a file-selectin window will
%             pop up to allow interactive file selection
% varargin    Variable number of arguments. Each argument is a cell array whose
%             first element is a keyword and whose other elements can be
%             strings, numeric values, or other cell arrays. Possible keywords  
%             are:  
%        'headers' header values to be written to the binary trace header in 
%             addition to those written by default (if available); 
%             Headers stored by default are:
%               ds_seqno      Trace sequence number within line (1-4)
%               ffid          Original Field record number (9-12)
%               o_trace_no    Trace sequence number within original field record (13-16)
%               source        Energy source point number (17-20)
%               cdp           CDP ensemble number (21-24)
%               seq_cdp       Trace sequence number within CDP ensemble (25-28)
%               trc_type      Trace ID (1=live,2=dead,3=dummy,4=time break,...) (29-30)
%               offset        Distance from source point to receiver group (37-40)
%               depth         Source depth below surface (49-52)
%               sou_h2od      Water depth at source (61-64)
%               rec_h2od      Water depth at receiver group (65-68)
%               sou_x         X coordinate of source (73-76)
%               sou_y         Y coordinate of source (77-80)
%               sou_elev      surface elevation at source (45-48)
%               rec_x         X coordinate of receiver (81-84)
%               rec_y         Y coordinate of receiver (85-88)
%               rec_elev      receiver elevation (41-44);
%               lag           Lag time between shot and recording start in ms (109-110)
%                  (the value of lag is added to the start time of the 
%                  seismic; hence it can be used to simulate non-zero start
%                  time of the data)
%               cdp_x         X coordinate of CDP (181-184)
%               cdp_y         Y coordinate of CDP (185-189)
%               iline_no      In-line number (189-192)
%               xline_no      Cross-line number (193-196)
%               The numbers in parentheses at the end of the line denote the location 
%               of the corresponding bytes in the SEG-Y trace header
%
%               The cell array for the user-specified headers has the form 
%               {'headers',{mnem1,first,bytes},{mnem2,first,bytes},...} where                
%               "mnem1", "mnem2", ... denote header mnemonics (such as ILINE_NO, XLINE_NO), 
%               "first" denotes the first byte in the binary header, 
%               "bytes" denotes the number of bytes occupied by the mnemonic (2 or 4), 
%               Example: {'headers',{'ILINE_NO',181,4},{'XLINE_NO',185,4}}
%                         these headers must, of course, be present in the input data set.
%               Default: no headers other than the standard SEG-Y headers listed above
%                         will be saved.
%        'print'   printout of messages. possible values 0 or 1.
%                Default: {'print',1}   this means that messages are printed.
%        'ascii_header'  ASCII version of the 40-line, 3200 byte EBCDIC header
%                a default EBCDIC header will be inserted if this header is not supplied
% EXAMPLE
%        write_segy_file(seismic,filename,{'headers',{'iline_no',189,4},{'xline_no',193,4}});              

global ABORTED

if ~istype(seismic,'seismic');
   error(' The first input argument must be a seismic data set.')
end

%       Set default parameters
param.format='ieee';
param.headers={};
param.ascii_header='';
param.print=1;

%       Decode and assign input arguments
param=assign_input(param,varargin);

if nargin < 2
   filename='';
end

%       Add zeros to make start time zero
if seismic.first > 0
   seismic=s_select(seismic,{'times',0,seismic.last});
elseif seismic.first < 0
   display([' WARNING! Start time of seismic < 0. Seimic will be shifted by ', ...
      num2str(-seismic.first),' ms']);
end

%       Define file format
if strcmpi(param.format,'ieee')
   param.format='ieee-be';
   o_format='ieee-be';
   datafmt=5;	% Data-format code
   fidx=0;
end

%      Create additional headers for output
[nsamp,ntr]=size(seismic.traces);
dt=seismic.step*1000;
seismic=s_header(seismic,'add','gen_scale_factor',1,'','Scale factor');
seismic=s_header(seismic,'add','number_of_samples',nsamp,'','Number of samples in this trace');
seismic=s_header(seismic,'add','samp_int',dt,'us','Sample interval of this trace in us');
  
 
%   Handle nulls in traces
if isfield(seismic,'null')
   if isnan(seismic.null)
%      idx=find(isnan(seismic.traces));
      seismic.traces(isnan(seismic.traces))=0;
   end
end

%   Handle nulls in headers
if isfield(seismic,'header_null')
   if isnan(seismic.header_null)
%     idx=find(isnan(seismic.headers));
     seismic.headers(isnan(seismic.headers))=0;
   end
end   

%  Open file
if ~isempty(filename)
   fid=fopen(filename,'w',o_format);
   if fid < 0
      disp(['... unable to create requested file "', filename,'"']);
   end
else
   fid=-1;
   filename=[seismic.name,'.sgy'];
end
if fid == -1
   [selected_file,ierr]=get_filename4w('sgy',filename);
   if ierr
      return;  
   end
   fid=fopen(selected_file,'w',o_format);
   if fid < 0
      disp(['... unable to create requested file "', filename,'"']);
      ABORTED=logical(1);
      return
   end
%   [pathname,filename]=fileparts(selected_file);
%   S4M.seismic_path=pathname;
end 

%    Write EBCDI reel header
if isempty(param.ascii_header)
   ascii_header=make_header; 
else
   ascii_header=param.ascii_header;
   [nh1,mh1]=size(ascii_header);
   if nh1*mh1 ~= 3200
      error('ASCII/EBCDIC header nust have 3200 bytes')
   else
      if nh1 > mh1
         ascii_header=ascii_header';
      end
   end
end
fwrite(fid,ascii2ebcdic(ascii_header),'char');
if param.print, disp('EBCDIC reel header written'), end;

%   Write binary reel header
two_bytes=zeros(194,1);

if isfield(seismic,'job_id')
  jobid=seismic.job_id;
else
  jobid=1;
end 

if isfield(seismic,'line_number')
  lineid=seismic.line_number;
else
  lineid=1;
end 

if isfield(seismic,'reel_number')
  reelid=seismic.reel_number;
else
  reelid=1;
end 

if isfield(seismic,'traces_per_record')
  two_bytes(1)=seismic.traces_per_record;
else
  two_bytes(1)=ntr;
end 

if isfield(seismic,'aux_per_record')
   two_bytes(2)=seismic.aux_per_record;
else
   two_bytes(2)=0;
end 

two_bytes(3:7)=[dt,dt,nsamp,nsamp,datafmt]';

if isfield(seismic,'cdp_fold')
   two_bytes(8)=seismic.cdp_fold;
end

if isfield(seismic,'sort')
   two_bytes(9)=seismic.sort;
end

if ~isfield(seismic,'headers')
    two_bytes(22)=0;
else 
  idx=find(ismember(seismic.header_info(:,1), ...
        {'offset','sou_x','rec_x','cdp_x','sou_y','rec_y','cdp_z'}));
  if isempty(idx)
    two_bytes(22)=0;
  else
    units=unique(seismic.header_info(idx,2));
    if length(units) > 1
       disp(' Units of measurement for distance:')
       disp(units)
       warning('Inconsistent distance units in headers; units in SEG-Y file are set to unknown');
    elseif strcmpi(char(units),'m')
      two_bytes(22)=1;
    elseif strcmpi(char(units),'ft')
      two_bytes(22)=2;
    else
      warning('Distance units in headers are neither ft nor m; units in SEG-Y file are set to unknown');