getnc_s.m

来自「读取Network Common Data Form (netCDF)数据」· M 代码 · 共 1,042 行 · 第 1/3 页

  error(['** ERROR ** ncvarinq: rcode = ' num2str(rcode)])
end

% Turn off the rescaling of the byte type data because ncmex 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 corner, end_point, stride and order.
% If there were 2 input arguments then set the values of corner, end_point
% and stride to their default values.  If there were 4 input arguments
% then set the value of stride to its default value.  Otherwise, check
% that the last input arguments are acceptable.  Also set take_stride
% which specifies whether strides need to be taken.  The cases where
% corner, end_point 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.  Finally check that the value of order is
% acceptable.

if nvdims == 0
  corner = 1;
  end_point = 1;
  stride = 1;
else
  if nargin == 2
    corner = ones(nvdims, 1);
    end_point = -1*ones(nvdims, 1);
    stride = ones(nvdims, 1);
  elseif nargin == 4
    if length(corner) == 1
      if corner < 0
	corner = ones(nvdims, 1);
        end_point = -1*ones(nvdims, 1);
      end
    elseif length(corner) ~= nvdims
      error('The corner vector is the wrong length')
    end
    stride = ones(nvdims, 1);
    if (sum(abs(size(corner) - size(end_point)))) ~= 0
      error('The sizes of corner and end_point do not agree')
    end
  else
    if length(corner) == 1
      if corner < 0
	corner = ones(nvdims, 1);
        end_point = -1*ones(nvdims, 1);
      end
    elseif length(corner) ~= nvdims
      error('The 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
      error('The stride vector is the wrong length')
    end
  end
end

corner_min = min(size(corner));
corner_max = max(size(corner));
end_point_min = min(size(end_point));
end_point_max = max(size(end_point));
stride_min = min(size(stride));
stride_max = max(size(stride));
if corner_min ~= end_point_min | corner_min ~= stride_min | ...
		 corner_max ~= end_point_max | corner_max ~= stride_max
  error('The sizes of corner, end_point and stride do not agree')
end

% Set take_stride.

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

% Make corner, end_point, stride and order into column vectors.

corner = corner(:);
end_point = end_point(:);
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

% Check the values of change_miss, new_miss and squeeze_it.

if length(change_miss) ~= 1
  error('ERROR: change_miss must be a scalar')
end

if all ( change_miss ~= [ 1 2 3 ] )
  if change_miss < 0
    change_miss = 2;
  else
   error([ 'getnc_s was passed change_miss = ' int2str(change_miss) ])
 end
end

if length(new_miss) ~= 1
  error('ERROR: new_miss must be a scalar')
end
  
if length(squeeze_it) ~= 1
  error('ERROR: squeeze_it must be a scalar')
end

if any ( ~isreal(squeeze_it) )
   error([ 'getnc_s was passed squeeze_it = ' num2str(squeeze_it) ])
end
       
% Find out whether to return a scalar, vector or matrix.  It is here
% that corner is decremented and edge is calculated so that the c-style
% conventions in ncmex will be followed.

if nvdims == 0
  edge = 1;
else
  edge = ones(nvdims, 1);
  for i = 1:nvdims
    dimid = vdims(i);
    [name, sizem, rcode] = ncmex('ncdiminq', cdfid, dimid);
    if rcode == -1
      error(['** ERROR ** ncdiminq: rcode = ' num2str(rcode)])
    end
    if ( corner(i) < 0 |  end_point(i) < 0 )
      corner(i) = 0;
      edge(i) = sizem;
    else
      % Check that corner & end_point 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 corner and end point values.
    
      corner(i) = corner(i) - 1;
      end_point(i) = end_point(i) - 1;
      if corner(i) >= sizem | end_point(i) < 0 | end_point(i) >= sizem
	s = [ 'getnc_s was passed corner = ' int2str(corner(i)+1) ...
	      ' & end_point = ' int2str(end_point(i)+1) ...
	      ' for dimension ' name ];
	error(s)
      end
      if stride(i) > 1
	edge(i) = fix( ( end_point(i) - corner(i) )/stride(i) ) + 1;
      else
	edge(i) = end_point(i) - corner(i) + 1;
      end
    end
  end
end
num_edge = length( find(edge ~= 1) );

if num_edge == 0

  % Get the scalar.

  [values, rcode] = ncmex('ncvarget1', cdfid, varid, corner, rescale_var);
  if rcode == -1
    error(['** ERROR ** ncvarget1: rcode = ' num2str(rcode)])
  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] = ncmex('ncvargetg', cdfid, varid, corner, edge, stride, imap, rescale_var);
    if rcode == -1
      error(['** ERROR ** ncvargetg: rcode = ' num2str(rcode)])
    end
  else
    [values, rcode] = ncmex('ncvarget', cdfid, varid, corner, edge, rescale_var);
    if rcode == -1
      error(['** ERROR ** ncvarget: rcode = ' num2str(rcode)])
    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 == -1
    values = permute(values, (ndims(values):-1:1));
  elseif order ~= -2
    values = permute(values, (ndims(values):-1:1));
    values = permute(values, order);
  end
   
  % Squeeze the array if required.
   
  if squeeze_it ~= 0
    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] = ncmex('attget', cdfid, varid, 'scale_factor');
      if rcode == -1
	error(['** ERROR ** ncattget: rcode = ' num2str(rcode)])
      end
    end
    pos = check_st('add_offset', attstring, nvatts);
    if pos > 0
      [addoff, rcode] = ncmex('attget', cdfid, varid, 'add_offset');
      if rcode == -1
	error(['** ERROR ** ncattget: rcode = ' num2str(rcode)])
      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.
  % Note 5: An earlier version of this software checked for an attribute
  % named missing_value.  This check was taken out because,
  % although in common use, missing_value was not given in the netCDF
  % manual list of attribute conventions.  Since it has now appeared in
  % the netCDF manual I have put the check back in.
  
  % The indices of the data points containing missing value indicators
  % will be stored separately in index_miss_low, index_miss_up, 
  % index_missing_value and index__FillValue.
  
  index_miss_low = [];
  index_miss_up = [];
  index__FillValue = [];
  index_missing_value = [];
  miss_low_orig = [];
  miss_up_orig = [];
  fill_value_orig = [];
  
  % First find the indices of the data points that are outside the valid
  % range.
  
  pos_vr = check_st('valid_range', attstring, nvatts);
  if pos_vr > 0
    [attype, attlen, rcode] = ncmex('ncattinq', cdfid, varid, 'valid_range');
    if rcode == -1
      error(['** ERROR ** ncattinq: rcode = ' num2str(rcode)])
    end
    [ miss, rcode] = ncmex('ncattget', cdfid, varid, 'valid_range');
    if rcode == -1
      error(['** ERROR ** ncattget: rcode = ' num2str(rcode)])
    end
      
    % Check that valid_range is a 2 element vector.
    
    if length(miss) ~= 2
      error(['The valid_range attribute must be a vector'])
    end
    
    % Correct for possible faulty handling of byte type
    
    if attype == nc_byte
      if miss(1) > 127; miss(1) = miss(1) - 256; end
      if miss(2) > 127; miss(2) = miss(2) - 256; end
    end

    miss_low = miss(1);
    miss_up = miss(2);
    miss_low_orig = miss_low;
    miss_up_orig = miss_up;
    
    % Rescale & add offsets if required.
    
    if rescale_att == 1
      if isempty(scalef) == 0
	miss_low = miss_low*scalef;
	miss_up = miss_up*scalef;
      end
      if isempty(addoff) == 0
	miss_low = miss_low + addoff;