getnc_s.m

matlacb程序包

function values = getnc_s(cell_args)
% getnc_s is called by getnc in the script or non-interactive case. The
% arguments passed to getnc have been bundled into the cell named cell_args
% which is passed to getnc_s. For documentation do a help on getnc.
  
%
% AUTHOR:   J. V. Mansbridge, CSIRO
%---------------------------------------------------------------------
%
%     Copyright (C), J.V. Mansbridge, 1992
%     Commonwealth Scientific and Industrial Research Organisation
%     $Id: getnc_s.m Mon, 03 Jul 2006 17:16:40 $
% 
%--------------------------------------------------------------------

% The function getnc_s calls the following external functions:
%      check_nc.m, check_st.m, fill_att.m, get_dods_dds.m, loaddap or
%      loaddods, mexnc.m, pos_cds.m, uncmp_nc.m, y_rescal.m.
% getnc_s also calls the internal functions :
%      remove_quotes, error_handle, parse_args, substitute_new.
%
% This function is called by: getnc.m

% Written by Jim Mansbridge december 10 1991

% In November 1998 some code was added to deal better with byte type data. Note
% that any values greater than 127 will have 256 subtracted from them. This is
% because on some machines (an SGI running irix6.2 is an example) values are
% returned in the range 0 to 255. Note that in the fix the values less than 128
% are unaltered and so we do not have to know whether the particular problem has
% occurred or not; for machines where there is no problem no values will be
% altered. This is applied to byte type attributes (like _FillValue) as well as
% the variable values.

% Check that there are the correct number of arguments. If the 6th, 7th
% or 8th arguments are missing then they are set here.  If the 3rd, 4th
% or 5th arguments are missing then their defaults will be set later
% when we find out the dimensions of the variable.  It produces an error
% if the bl_corner is set but the tr_corner is not.

[file, varid, bl_corner, tr_corner, stride, order, change_miss, new_miss, ...
 squeeze_it, rescale_opts, err_opt] = parse_args(cell_args);

% Set whether values will be automatically rescaled or not.

if isnumeric(rescale_opts)
  if ndims(rescale_opts) ~= 2
    error('rescale_opts must be [], -1 or a 2 element vector');
  end
  len_rescale_opts = length(rescale_opts);
  if len_rescale_opts == 2
    [rescale_var, rescale_att] = y_rescal(rescale_opts);
  elseif len_rescale_opts < 2
    [rescale_var, rescale_att] = y_rescal;
  else
    error('rescale_opts must be [], -1 or a 2 element vector')
  end
else 
  error('rescale_opts must be numeric');
end

% Set some constants.

blank = abs(' ');

% Do some initialisation.

CSIRO_add_jar_file_maybe;
try
  [mex_name, full_name, desc_das, file_status, exe_name] = ...
      choose_mexnc_opendap(file);
catch
  cdfid = -1000;
  mess_str = lasterr;
  rcode = -1001;
  values = error_handle(cdfid, mess_str, rcode, err_opt);
  return
end
  

switch mex_name
 case 'mexnc'
  
  % I make mexnc calls to find the integers that specify the attribute
  % types

  nc_byte = mexnc('parameter', 'nc_byte'); %1
  nc_char = mexnc('parameter', 'nc_char'); %2
  nc_short = mexnc('parameter', 'nc_short'); %3
  nc_long = mexnc('parameter', 'nc_long'); %4
  nc_float = mexnc('parameter', 'nc_float'); %5
  nc_double = mexnc('parameter', 'nc_double'); %6

  % Set the value of imap.  Note that this is used simply as a
  % placeholder in calls to vargetg - its value is never used.

  imap = 0;

  % Open the netcdf file.

  [cdfid, rcode] = mexnc('ncopen', full_name, 'NC_NOWRITE');
  if rcode ~= 0
    mess_str = ['ncopen ' full_name];
    values = error_handle(cdfid, mess_str, rcode, err_opt);
    return
  end

  % Suppress all error messages from netCDF 

  [rcode] = mexnc('setopts', 0);

  % Collect information about the cdf file.

  [ndimens, nvars, ngatts, recdim, rcode] =  mexnc('ncinquire', cdfid);
  if rcode < 0
    mess_str = ['ncinquire: ' full_name];
    values = error_handle(cdfid, mess_str, rcode, err_opt);
    return
  end

  % Determine the netcdf id number for the required variable.  If varid is a
  % string then an appropriate call to mexnc is used to convert it to the
  % relevant integer. If varid is a number then it is decremented.  This is done
  % because inqnc & getnc count the variables from 1 to nvars whereas the calls
  % to the mexnc routines use c-type counting, from 0 to nvars - 1.

  varid_old = varid;
  if ischar(varid)
    [varid, rcode] = mexnc('ncvarid', cdfid, varid_old);
    if rcode < 0
      mess_str = ['ncvarid: ' full_name ': variable named ' varid_old];
      values = error_handle(cdfid, mess_str, rcode, err_opt);
      return
    end
  else
    varid = varid - 1;
  end

  % Check the value of varid

  if varid < 0 | varid >= nvars
    error([ 'getnc_s was passed varid = ' int2str(varid) ])
  end

  % Find out info about the variable, in particular find nvdims, the number of
  % dimensions that the variable has.

  [varnam, vartypv, nvdims, vdims, nvatts, rcode] = ...
      mexnc('ncvarinq', cdfid, varid);
  if rcode < 0
    mess_str = ['ncvarinq: ' full_name ': variable named ' varid_old];
    values = error_handle(cdfid, mess_str, rcode, err_opt);
    return
  end

  % Turn off the rescaling of the byte type data because mexnc does not do this
  % for variables anyway. The rescaling of the VALUES array will be done
  % explicitly.

  if vartypv == nc_byte
    rescale_var = 0;
    rescale_att = 0;
  end
  
  % Do checks on bl_corner, tr_corner and stride. The cases where
  % bl_corner, tr_corner and stride are -1 or -1*ones(nvdims, 1) are checked
  % for and handled here.  Note that stride may also be a scalar 1 when a
  % vector may seem to be required.

  if nvdims == 0
    bl_corner = 1;
    tr_corner = 1;
    stride = 1;
  else
    if length(bl_corner) == 1
      if bl_corner < 0
	bl_corner = ones(nvdims, 1);
	tr_corner = -1*ones(nvdims, 1);
      end
    elseif length(bl_corner) ~= nvdims
      error('The bl_corner vector is the wrong length')
    end
    
    if length(stride) == 1
      if stride < 0 | stride == 1
	stride = ones(nvdims, 1);
      end
    elseif length(stride) == nvdims
      ff = find(stride < 0);
      if ~isempty(ff)
	stride(ff) = 1;
      end
    else
      error('The stride vector is the wrong length')
    end
  end
  
  bl_corner_min = min(size(bl_corner));
  bl_corner_max = max(size(bl_corner));
  tr_corner_min = min(size(tr_corner));
  tr_corner_max = max(size(tr_corner));
  stride_min = min(size(stride));
  stride_max = max(size(stride));
  if bl_corner_min ~= tr_corner_min | bl_corner_min ~= stride_min | ...
	bl_corner_max ~= tr_corner_max | bl_corner_max ~= stride_max
    error('The sizes of bl_corner, tr_corner and stride do not agree')
  end

  % Set take_stride which specifies whether strides need to be taken.

  if max(stride) > 1
    take_stride = 1;
  else
    take_stride = 0;
  end

  % Make bl_corner, tr_corner, stride and order into column vectors.

  bl_corner = bl_corner(:);
  tr_corner = tr_corner(:);
  stride = stride(:);
  order = order(:);

  % Check order

  if length(order) == 1
    if order == 1 % Special case where the netcdf variable is a vector
      order = -1;
    elseif order ~= -1 & order ~= -2
      error('ERROR: if order is a scalar it must be -1 or -2')
    end
  else
    if length(order) ~= nvdims
      error('The order vector is the wrong length')
    elseif sum(abs(sort(order) - (1:nvdims)')) ~= 0
      error(['The order vector must be a rearrangement of the numbers 1 to ' ...
	     num2str(nvdims)])
    end
  end

  % Find out whether to return a scalar, vector or matrix.  It is here
  % that bl_corner is decremented and edge is calculated so that the c-style
  % conventions in mexnc will be followed. We also check whether any of the
  % dimensions have zero length (this is certainly possible for the record
  % dimension).

  if nvdims == 0
    edge = 1;
  else
    edge = ones(nvdims, 1);
    for i = 1:nvdims
      dimid = vdims(i);
      [name, sizem, rcode] = mexnc('ncdiminq', cdfid, dimid);
      if rcode < 0
	mess_str = ['ncdiminq: ' full_name ': dimid = ' dimid];
	values = error_handle(cdfid, mess_str, rcode, err_opt);
	return
      end
      
      % Check for a zero length dimension (this is certainly possible for the
      % record dimension). If so then just return an empty array.
      
      if sizem == 0
	mess_str = [varnam ' has dimension ' name ' which has zero length'];
	values = error_handle(cdfid, mess_str, rcode, err_opt);
	return
      end
      
      % Apply the deaults if required.
      
      if bl_corner(i) < 0
	bl_corner(i) = 1;
	tr_corner(i) = sizem;
	stride(i) = 1;
      end
      if tr_corner(i) < 0
	tr_corner(i) = sizem;
      end
      
      % Check that bl_corner & tr_corner are in the correct range.  If they
      % are then calculate edge.  Note that because I am using the
      % matlab & fortran conventions for counting indices I must
      % subtract 1 from the bl_corner and end point values.
	
      bl_corner(i) = bl_corner(i) - 1;
      tr_corner(i) = tr_corner(i) - 1;
      if bl_corner(i) >= sizem | tr_corner(i) < 0 | tr_corner(i) >= sizem
	s = [ 'getnc_s was passed bl_corner = ' int2str(bl_corner(i)+1) ...
	      ' & tr_corner = ' int2str(tr_corner(i)+1) ...
	      ' for dimension ' name ];
	error(s)
      end
      if stride(i) > 1
	edge(i) = fix( ( tr_corner(i) - bl_corner(i) )/stride(i) ) + 1;
      else
	edge(i) = tr_corner(i) - bl_corner(i) + 1;
      end
    end
  end
  num_edge = length( find(edge ~= 1) );

  if num_edge == 0

    % Get the scalar.

    [values, rcode] = mexnc('ncvarget1', cdfid, varid, bl_corner, rescale_var);
    if rcode < 0
      mess_str = ['ncvarget1: ' full_name ': varid =  ' varid];
      values = error_handle(cdfid, mess_str, rcode, err_opt);
      return
    end
    
    % Do possible byte correction.
    
    if vartypv == nc_byte
      ff = find(values > 127);
      if ~isempty(ff)
	values(ff) = values(ff) - 256;
      end
    end

  else

    % Get the full hyperslab and return it as an array of the appropriate
    % dimensions.  Note that we must allow for the C-type notation where
    % the fastest changing index is the last mentioned.

    if take_stride
      [values, rcode] = mexnc('ncvargetg', cdfid, varid, bl_corner, edge, ...
			      stride, imap, rescale_var);
      if rcode < 0
	mess_str = ['ncvargetg: ' full_name ': varid =  ' varid];
	values = error_handle(cdfid, mess_str, rcode, err_opt);
	return
      end
    else
      [values, rcode] = mexnc('ncvarget', cdfid, varid, bl_corner, edge, ...
			      rescale_var);
      if rcode < 0
	mess_str = ['ncvarget: ' full_name ': varid =  ' varid];
	values = error_handle(cdfid, mess_str, rcode, err_opt);
	return
      end
    end
    
    % Do possible byte correction.
    
    if vartypv == nc_byte
      ff = find(values > 127);
      if ~isempty(ff)
	values(ff) = values(ff) - 256;
      end
    end

    % Permute the array as required.  Note that the default behaviour is to
    % reverse the order of the indices to map between the matlab and C
    % conventions for ordering indices. If order is a vector then this
    % reversing must be done before the permutation using the vector.
    
    if length(order) == 1
      if order == -1
	values = permute(values, (ndims(values):-1:1));
      end
    else
      values = permute(values, (ndims(values):-1:1));
      values = permute(values, order);
    end
    
    % Squeeze the array if required.
    
    if squeeze_it
      values = squeeze(values);
    end

    % After squeezing a vector may be a row or column vector and so
    % turn any row vector into a column vector for consistency.
    
    if ndims(values) == 2
      [m_temp, n_temp] = size(values);
      if m_temp == 1
	values = values(:);
      end
    end
  end

  % If the missing values are to be replaced then do it here.

  scalef = [];
  addoff = [];
  if change_miss ~= 1

    % Find any scale factors or offsets.

    attstring = fill_att(cdfid, varid, nvatts);
    if rescale_att == 1 | vartypv == nc_byte
      pos = check_st('scale_factor', attstring, nvatts);
      if pos > 0
	[scalef, rcode] = mexnc('ncattget', cdfid, varid, 'scale_factor');
	if rcode < 0
	  mess_str = ['ncattget: ' full_name ': varid =  ' varid ...
		      ': scale_factor'];
	  values = error_handle(cdfid, mess_str, rcode, err_opt);
	  return
	end
      end
      pos = check_st('add_offset', attstring, nvatts);
      if pos > 0
	[addoff, rcode] = mexnc('ncattget', cdfid, varid, 'add_offset');
	if rcode < 0
	  mess_str = ['ncattget: ' full_name ': varid =  ' varid ...
		      ': add_offset'];
	  values = error_handle(cdfid, mess_str, rcode, err_opt);
	  return
	end
      end
    end

    % check for missing values.  Note that a
    % missing value is taken to be one less than valid_min, greater than
    % valid_max or 'close to' _FillValue or missing_value.
    % Note 1: valid_min and valid_max may be specified by the attribute
    % valid_range and if valid_range exists than the existence of
    % valid_min and valid_max is not checked.
    % Note 2: a missing value must be OUTSIDE the valid range to be
    % recognised.
    % Note 3: a range does not make sense for character arrays.
    % Note 4: By 'close to' _FillValue I mean that an integer or character
    % must equal _FillValue and a real must be in the range
    % 0.99999*_FillValue tp 1.00001*_FillValue.  This allows real*8 
    % rounding errors in moving the data from the netcdf file to matlab;
    % these errors do occur although I don't know why given that matlab
    % works in double precision.