read_las_file.m

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

function [wlog,las_header]=read_las_file(filename,iprint)
% Function reads well logs from a disk file in LAS (Log ASCII Standard) Version 2.0
% format as specified by the Canadian Well Logging Society
%
% Written by: E. R.: February 6, 2000;  
% Last updated: April 13, 2006: Create fielfd "from" with full filename
%
%          [log,las_header]=read_las_file(filename,iprint)
% INPUT
% filename  string with the file name (and full path if desired);
%         if the file is not found a file selection box will pop up to
%         allow interactive file selection
% iprint  Control diagnostic output (optional)
%         iprint = 0  ==> no output; DEFAULT
%         iprint ~= 0 ==> print out progress of file reading
% OUTPUT
% wlog    structure containing the log curves and ancillary information
%    wlog.curves  Matrix of curve values
%		 Each log is in a column and the first column is usually the depth.
%                Thus, if w1log.curves has n columns, then there are n-1 log curves.
%                If necessary, the rows of curve data are resorted so that
%                depth increases from one row to the next
%         
%    wlog.curve_info Cell array (3 x number of curves) with curve mnemonics, 
%                curve units of measurement, and curve descriptions
%    wlog.first	 Start of log (first depth in file)
%    wlog.last   	 End of log (last depth in file)
%    wlog.step	 Depth increment (0 if unequal)
%    wlog.units  Units of measurement for depth
%    wlog.null	 Null value (set to NaN) if there are no-data values
%                otherwise, this field is not set
%    wlog.wellname  Name of well
%    wlog.location  Location of well
%    wlog.company  Oil/gas company which drilled the well
%    wlog.field   Field in which well is located
% 
%                The following components, while desirable, may or may not be present 
%                (depending on LAS file), and others may be present. 
%                They include well identification such as
%    wlog.country Country in which well is located            
%    wlog.wellid	 Well Identifier
%    wlog.api     API number of the well
%    wlog.service Service company which logged the well
%
%                Also included are all parameters in the Parameter Section of the
%                LAS header which have numerical values other than the null value. 
%                Each parameter is stored in a field of the log structure; its 
%                units of measurement and a description are stored in a row of
%                the cell array "wlog.parameter_info" 
%                Examples illustrating such parameters are 
%    wlog.ekb=84	Kelly bushing elevation     
%    wlog.egl=0		Ground level elevation
%    wlog.parameter_info={'ekb','ft','Kelly bushing elevation';
%                         'egl','ft','Ground level elevation'}
% las_header  LAS header of log (provides additional information about the well)

global S4M

run_presets_if_needed

if nargin  <= 1;
   iprint=0;
   if nargin == 0
      filename='';
   end
end  

%  Open the file
if ~isempty(filename)
   fid=fopen(filename,'rt');
   if fid > 0
      filename2S4M(filename)
   end
else
   fid=-1;
end

if fid == -1 
   [filename,ierr]=get_filename4r('las');
   if ierr
      error('No file selected')
   end

   fid=fopen(filename,'rt');
   if fid < 0
      error(['File "',filename,'" could not be opened.'])
   end
end 

lashead=repmat(32,3,100);
% bmax=0;
nhead=1;

buffer=get_buffer(fid);   
lashead(1,1:length(buffer))=buffer;
% if length(buffer) > bmax, 
%   bmax=length(buffer); 
% end

while strcmpi(buffer(1),'#')
   buffer=get_buffer(fid);
   nhead=nhead+1;
   lashead(nhead,1:length(buffer))=buffer;
end

if ~strcmpi(buffer(1:2),'~v')
     error('File is not a LAS file');
end

% 	Determine the version of LAS
buffer=get_buffer(fid);
nhead=nhead+1;
lashead(nhead,1:length(buffer))=buffer;
% if length(buffer) > bmax, bmax=length(buffer); end
ind=find(buffer == ':')-1;
ind1=find(buffer == '.')+1;
%vernum=sscanf(buffer(ind1(1):ind),'%*s%f');
vernum=sscanf(buffer(ind1(1):ind),'%f');

%	Check for wrapping
buffer=get_buffer(fid);
nhead=nhead+1;
lashead(nhead,1:length(buffer))=buffer;
% if length(buffer) > bmax, bmax=length(buffer); end
ind=find(buffer == ':')-1;
wrapflag=sscanf(buffer(1:ind),'%*s%s');
if strcmpi(wrapflag,'YES')
   disp(' LAS file is wrapped')
%  disp([mfilename,' cannot read wrapped LAS files'])
%  error('  use REFORMAT.EXE from http://www.cwls.org to unwrap file');
end

if vernum ~= 2.0
   disp([mfilename,' can only read Version 2.0 LAS files'])
   error('  use REFORMAT.EXE from http://www.cwls.org to convert from version 1.2 to 2.0');
else
   [wlog,las_header] = read_las_v20(fid,iprint);
end

las_header=char(strim(char(lashead)),las_header);

%       Check log start and end time and step size
if isempty(wlog.curves)
   disp(' No log curve values read; there may be non-numeric characters in the data block')
   error(' Check data section of LAS file')
end

wlog.units=wlog.curve_info{1,2};

if wlog.first ~= wlog.curves(1,1)
   fprintf(['Log start depth (',num2str(wlog.curves(1,1)), ...
      ') differs from header information (',num2str(wlog.first), ...
      '); now corrected\n'])
   wlog.first=wlog.curves(1,1);
end
if wlog.last ~= wlog.curves(end,1)
  fprintf(['Log end depth (',num2str(wlog.curves(end,1)), ...
    ') differs from header information (',num2str(wlog.last), ...
    '); now corrected\n']);
  wlog.last=wlog.curves(end,1);
end
if wlog.step ~= 0
   if ~isconstant(diff(wlog.curves(:,1)),S4M.log_step_error)
      fprintf('Log step size not uniform enough; step size changed to 0\n')
   end
end

%   Make sure that log depths are ascending
if wlog.first > wlog.last
   wlog.curves=flipud(wlog.curves);
   temp=wlog.first;
   wlog.first=wlog.last;
   wlog.last=temp;
   wlog.step=-wlog.step;
end

%  Replace null values by NaNs
idx=find(wlog.curves == wlog.null);
if ~isempty(idx)
   wlog.curves(idx)=NaN;
   wlog.null=NaN;
else
   wlog=rmfield(wlog,'null');
end

%	Replade 'UNDEFINED' in curve description when it is obvious from the mnemonic
wlog.curve_info=description_substitution(wlog.curve_info);

%       Replace unadmissible curve mnemonics
wlog.curve_info(:,1)=fix_mnemonics(wlog.curve_info(:,1));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [wlog,las_header]=read_las_v20(fid,iprint)
% Function reads LAS file version 2.0
% Date Feb. 12, 2000;  written by E. R.
%      Feb. 19, 2000;  updated to handle better missing key words 
%
% INPUT
% fid    File id of LAS file
% iprint  Control diagnostic output (optional)
%        iprint = 0  ==> No output; DEFAULT
%        iprint ~= 0 ==> diagnostic output
%     
% OUTPUT
% wlog      structure contining the log curves and ancillary information
%  wlog.curves   	Matrix of curve values
%		 Each log is in a column and the first column is usually the depth.
%                Thus, if wlog.curves has n columns, then there are n-1 log curves.
%                If necessary, the rows of curve data are resorted so that
%                depth increases from one line to the next
%         
%  wlog.curve_info Cell array (3 x number of curves) with curve mnemonics, 
%                curve units of measurement, and curve descriptions
%  wlog.first	 Start of log (first depth in file)
%  wlog.last	 End of log (last depth in file)
%  wlog.step	 Depth increment (0 if unequal)
%  wlog.null	 Null value (set to NaN)
%  wlog.las_header  LAS header of log (provides additional information about the well)
%  wlog.wellname	 Name of well
%  wlog.location	 Location of well
%  wlog.company   Oil/gas company which drilled the well
%  wlog.field     Field in which well is located
% 
%                The following components, while desirable, may or may not be present 
%                (depending on LAS file), and others may be present. 
%                They include well identification such as
%  wlog.country   Country in which well is located
%  wlog.county    County in which well is located
%  wlog.province  Province in which well is located
%  wlog.state     State in whichwell is located
%  wlog.uwi       Universal well identifier (Canada)            
%  wlog.wellid    Well Identifier
%  wlog.api       API number of the well (USA)
%  wlog.service   Service company which logged the well
%
%                Also included are all parameters in the Parameter Section of the
%                LAS header which have numerical values other than the null value. 
%                Each parameter is stored in a field of the log structure; its 
%                units of measurement and a description are stored in a row of
%                the cell array "wlog.parameter_info" 
%                Examples illustrating such parameters are 
%    wlog.ekb=84	Kelly bushing elevation     
%    wlog.egl=0		Ground level elevation
%    wlog.parameter_info={'ekb','ft','Kelly bushing elevation';
%                         'egl','ft','Ground level elevation'}

global S4M

wlog.type='well_log';
wlog.tag='unspecified';
[dummy,wlog.name]=fileparts(S4M.filename);  %#ok The first output argument 
                                            %    is not required
wlog.from=fullfile(S4M.pathname,S4M.filename);

% 	Allocate space for a 200x100-byte header
nchar=100;
% las_header=32*ones(200,nchar);
las_header=repmat(32,200,nchar);

mnem_length=25;	 	% Maximum length of curve mnemonic

%       Accomodate up to 100 curves and up to 100 parameters
%n_curves=100;
%n_param=100;

if nargin == 1; iprint=0; end
nhead=0;
bmax=0;

buffer=get_buffer(fid);


		while isempty(buffer) | ~strcmpi(buffer(1:2),'~a')

unknown_block=1;  % Used to check if an unknown block is in the file

% 	Well information block
if strcmpi(buffer(1:2),'~w')
  if iprint ~= 0, disp('... reading well info'); end
  unknown_block=0;
  nhead=nhead+1;
  lb=min([length(buffer),nchar]);
  bmax=max([lb,bmax]);
  las_header(nhead,1:lb)=buffer(1:lb);
  
  if iprint ~= 0, disp('... reading well info'); end
  buffer=get_buffer(fid);
  while buffer(1:1) ~= '~' 
    nhead=nhead+1;
    lb=min([length(buffer),nchar]);
    bmax=max([lb,bmax]);
    las_header(nhead,1:lb)=buffer(1:lb);
    ia=find(buffer == '.')+1;
    if length(ia) > 1
       ia=ia(1);
    end
    ind=find(buffer==':')-1;
    kw=deblank(buffer(1:ia-2));   
		      
%  	Check for start, stop and step 
    if strcmpi(kw,'STRT')
%      if buffer(ia) ~= ' '
      if ~strcmp(buffer(ia),' ')
        temp=sscanf(buffer(ia:ind),'%s',1);
        nu=length(temp);
      else
        nu=0;