getnc_s.m

matlacb程序包

    mess_str = ['ncclose: ' full_name];
    values = error_handle(cdfid, mess_str, rcode, err_opt);
    return
  end

 case {'loaddap', 'loaddods'}
  
  % Dealing with a dods file.
  
  % Check the value of varid and get its structure.
  
  if ischar(varid)
    if ~isfield(desc_das, varid)
      mess_str = [varid ' is not a variable in ' file];
      values = error_handle([], mess_str, [], err_opt);
      return
    end
  else
    error('For dods data varid must be a string')
  end
  
  % Get the structure of varid and then use this to extract required
  % information about attributes. Also get the dds of the file. Note that we
  % use the function remove_quotes to eliminate extraneous quotes returned as
  % part of the DAS. 
  
  %desc_varid = desc_das.(varid);
  desc_varid = getfield(desc_das, varid);
  if isfield(desc_varid, 'scale_factor')
    %scalef = desc_varid.('scale_factor');
    scalef = getfield(desc_varid, 'scale_factor');
    scalef = remove_quotes(scalef);
    if length(scalef) ~= 1
      error(['The scalef attribute must be a scalar'])
    end
  else
    scalef = [];
  end
  if isfield(desc_varid, 'add_offset')
    %addoff = desc_varid.('add_offset');
    addoff = getfield(desc_varid, 'add_offset');
    addoff =  remove_quotes(addoff);
    if length(addoff) ~= 1
      error(['The addoff attribute must be a scalar'])
    end
 else
    addoff = [];
  end
  if isfield(desc_varid, 'ml__FillValue')
    %fillvalue = desc_varid.('ml__FillValue');
    fillvalue = getfield(desc_varid, 'ml__FillValue');
    fillvalue =  remove_quotes(fillvalue);
    if rescale_att
      if ~isempty(scalef)
	fillvalue = fillvalue*scalef;
      end
      if ~isempty(addoff)
	fillvalue = fillvalue + addoff;
      end
    end
  else
    fillvalue = [];
  end
  if isfield(desc_varid, 'missing_value')
    %missing_value = desc_varid.('missing_value');
    missing_value = getfield(desc_varid, 'missing_value');
    missing_value =  remove_quotes(missing_value);
    if rescale_att
      if ~isempty(scalef)
	missing_value = missing_value*scalef;
      end
      if ~isempty(addoff)
	missing_value = missing_value + addoff;
      end
    end
  else
    missing_value = [];
  end
  if isfield(desc_varid, 'valid_min')
    %valid_min = desc_varid.('valid_min');
    valid_min = getfield(desc_varid, 'valid_min');
    valid_min =  remove_quotes(valid_min);
    if rescale_att
      if ~isempty(scalef)
	valid_min = valid_min*scalef;
      end
      if ~isempty(addoff)
	valid_min = valid_min + addoff;
      end
    end
  else
    valid_min = [];
  end
  if isfield(desc_varid, 'valid_max')
    %valid_max = desc_varid.('valid_max');
    valid_max = getfield(desc_varid, 'valid_max');
    valid_max =  remove_quotes(valid_max);
    if rescale_att
      if ~isempty(scalef)
	valid_max = valid_max*scalef;
      end
      if ~isempty(addoff)
	valid_max = valid_max + addoff;
      end
    end
  else
    valid_max = [];
  end
  if isfield(desc_varid, 'valid_range')
    %valid_range = desc_varid.('valid_range');
    valid_range = getfield(desc_varid, 'valid_range');
    valid_range =  remove_quotes(valid_range);
    if rescale_att
      if ~isempty(scalef)
	valid_range = valid_range*scalef;
      end
      if ~isempty(addoff)
	valid_range = valid_range + addoff;
      end
    end
  else
    valid_range = [];
  end
    
  % Get the dds of the file and use it to get some required information.
  
  [dds_text, desc_dds] = get_dods_dds(file, exe_name);
  
  for ii = 1:length(desc_dds.variable)
    if strcmp(desc_dds.variable(ii).name, varid)
      break
    end
  end
  if ~strcmp(desc_dds.variable(ii).name, varid)
    error('This error should not ever occur')
  end
  var_type = desc_dds.variable(ii).type;
  var_dim_statement = desc_dds.variable(ii).dim_statement;
  var_dim_idents = desc_dds.variable(ii).dim_idents;
  
  nvdims = length(var_dim_idents);
  dim_length = zeros(nvdims, 1);
  for ii = 1:nvdims
    dim_length(ii) = desc_dds.dimension(var_dim_idents(ii)).length;
    dim_name{ii} = desc_dds.dimension(var_dim_idents(ii)).name;
  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.

  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: the argument order must be -1 or -2')
    end
  else
    error(['ERROR: the argument order must be a scalar - the vector ' ...
	   ' option is not supported for opendap access'])
  end

  % Generate a string that specifies the hyperslab to be returned. In the
  % process check that the hyperslab is within bounds.
  
  str_h = [];
  for ii = 1:nvdims
    if bl_corner(ii) < 0
      bl_corner(ii) = 1;
      tr_corner(ii) = dim_length(ii);
      stride(ii) = 1;
    elseif bl_corner(ii) > dim_length(ii)
      error(['bl_corner so big that it is outside the hyperslab'])
    end
    if tr_corner(ii) < 0
      tr_corner(ii) = dim_length(ii);
    elseif tr_corner(ii) > dim_length(ii)
      error(['tr_corner so big that it is outside the hyperslab'])
    elseif tr_corner(ii) < bl_corner(ii)
      error(['tr_corner is less than bl_corner'])
    end
    
    if take_stride & (stride(ii) > 1)
      str_h = [str_h '[' num2str(bl_corner(ii) - 1) ':' num2str(stride(ii)) ...
	       ':' num2str(tr_corner(ii) - 1) ']'];
    else
      str_h = [str_h '[' num2str(bl_corner(ii) - 1) ':' ...
	       num2str(tr_corner(ii) - 1) ']'];
    end
  end
  
  % Get the hyperslab of data. Note that we allow for the vector of variable
  % values to be at different "depths" in the structure. It seems that a
  % variable that is also a dimension will be at values_struct.varid but
  % otherwise it will be at values_struct.varid.varid. We also use
  % method_of_call == 2 which, according to the documentation, does the same
  % thing as method_of_call == 1 but would be expected to be slower. However,
  % this is necessary for loaddods to work properly on Windows boxes when
  % an array is filled with characters. I have no idea what the bug might be
  % in loaddods. There is also a hint that the same problem can occur with
  % loaddap on some linux boxes.
  
  method_of_call = 2;
  switch method_of_call
   case 1
    switch mex_name
     case 'loaddap'
      values_struct = loaddap('+v', [full_name '?' varid str_h]);
     case 'loaddods'
      values_struct = loaddods('+v', [full_name '?' varid str_h]);
    end
   case 2
    switch mex_name
     case 'loaddap'
      loaddap('+v', [full_name '?' varid str_h]);
     case 'loaddods'
      loaddods('+v', [full_name '?' varid str_h]);
    end
    eval(['values_struct = ' varid ';'])
  end
  if isstruct(values_struct) 
    %xx = values_struct.(varid);
    xx = getfield(values_struct, varid);
    if isstruct(xx)
      %values = xx.(varid);
      values = getfield(xx, varid);
    else
      values = xx;
    end
  else
    values = values_struct;
  end
  
  % If we are getting an array of character data then loaddap returns a row
  % vector regardless of the original format. The following kluge is designed
  % to recreate the original array structure.
  
  kluge_char = 1;
  if kluge_char & ischar(values)
    dim_new = zeros(1, nvdims);
    for kk = 1:nvdims
      dim_new(nvdims - kk + 1) = length(bl_corner(kk):stride(kk):tr_corner(kk));
    end
    if nvdims == 1
      values = reshape(values, [1 dim_new]);
    elseif nvdims == 2
      values = reshape(values, dim_new);
    else
      values = reshape(values, dim_new);
      %error('Can''t handle code with a character array of more than 2 dimensions')
    end
  end
  
  % Squeeze the array if required.
  
  if squeeze_it
    values = squeeze(values);
  end

  % Now deal with permutation problems in the data.
  
  size_values = size(values);
  num_dim = length(size_values);
  if num_dim == 2
    if min(size_values) == 1
      % After squeezing a vector may be a row or column vector and so
      % turn any row vector into a column vector for consistency.
      
      values = values(:);
    else
      % Apply the order option to transpose the matrix
      if kluge_char & ischar(values)
	if order == -1
	  values = values';
	end
      else
	if order == -2
	  values = values';
	end
      end
    end
  elseif num_dim > 2
    
    % Permute the array as required.  There are 3 issues here.
    % 1) loaddap returns multi-dimensional numerical arrays in a strange order -
    % half-fortran and half-C conventions. First explain this conceptually. We
    % see the variable dimensions in a certain order when we look at the DDS -
    % the order is the same using inqnc or ncdump on a netcdf file. After
    % loaddap puts the data into a matlab variable it looks as though there
    % have been 2 transformations. First the order is totally reversed and
    % then the order of the first 2 dimensions are reversed
    % again. Schematically, suppose that a variable has N dimensions labelled
    % 1, 2, 3, ..., N-4, N-3, N-2, N-1, N in the DDS. The matlab order is then
    % N-1, N, N-2, N-3, N-4, ..., 3, 2, 1.
    % 2) loaddap returns multi-dimensional character arrays as vectors. It
    % simply streams out the data in the same order that it is stored on
    % disk.
    % 3) We must permute the array according to the value of the argument
    % "order" also.
    
    if kluge_char & ischar(values)
      if order == -1
	vec_permute = [num_dim:-1:1];
      elseif order == -2
	vec_permute = [1:num_dim];
      end
    else
      if order == -1
	vec_permute = [num_dim:-1:3 1 2];
      elseif order == -2
	vec_permute = [2 1 3:num_dim];
      end
    end
    values = permute(values, vec_permute);
  end

  % Do possible byte correction.
    
  if strcmp(lower(var_type), 'byte')
    ff = find(values > 127);
    if ~isempty(ff)
      values(ff) = values(ff) - 256;
    end
  end

  % Possibly do a rescale of the variable.
  
  if rescale_var
    if ~isempty(scalef)
      values = values*scalef;
    end
    if ~isempty(addoff)
      values = values + addoff;
    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 = [];
  if ~isempty(findstr(lower(var_type), 'float')) | ...
	~isempty(findstr(lower(var_type), 'int'))
    var_is_number = 1;
  else
    var_is_number = 0;
  end