getnc.m

读取Network Common Data Form (netCDF)数据

	 end

	 s = [ s '  ' stri ') ' varnam ];
         addit = 52 - length(s);
         for j =1:addit
            s = [ s ' '];
         end
      end 

      if i < nvars - 2
         stri = int2str(i+3);
         if length(stri) == 1
            stri = [ ' ' stri];
         end
         [varnam, vartyp, nvdims, vdims, nvatts, rcode] = ...
	     ncmex('ncvarinq', cdfid, i+2);
	 if rcode == -1
	   error(['** ERROR ** ncvarinq: rcode = ' num2str(rcode)])
	 end
         s = [ s '  ' stri ') ' varnam ];
      end 
      disp(s)
   end
   disp(' ')
   s = [ 'Select a menu number: '];
   k = return_v(s, -1);
end

% try to get information about the variable

varid = k - 1;
[varnam, vartypv, nvdims, vdims, nvatts, rcode] = ...
    ncmex('ncvarinq', cdfid, varid);
if rcode == -1
  error(['** ERROR ** ncvarinq: rcode = ' num2str(rcode)])
end
attstring = fill_att(cdfid, varid, nvatts);

% 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

if nvdims > 0
  message = cell(nvdims, 1);
  name_dim = cell(nvdims, 1);
end

for i = 1:nvdims
    dimid = vdims(i);
    [name, sizem, rcode] = ncmex('ncdiminq', cdfid, dimid);
    if rcode == -1
      error(['** ERROR ** ncdiminq: rcode = ' num2str(rcode)])
    end

   name_dim{i, 1} = name;
   ledim(i) = sizem - 1;

   % Test that the dimension name is also a variable name.  If it is then
   % store information about its initial and final values in the string s.

    rhid = check_st(name, varstring, nvars) - 1;

   if rhid >= 0
      [namejunk, dvartyp, dnvdims, vdimsjunk, nvattsjunk, rcode] = ...
	  ncmex('ncvarinq', cdfid, rhid);
      if rcode == -1
	error(['** ERROR ** ncvarinq: rcode = ' num2str(rcode)])
      end
      if sizem <= 6
         [temp, rcode] = ncmex('ncvarget', cdfid, rhid, [0], [sizem], rescale_var);
	 if rcode == -1
	   error(['** ERROR ** ncvarget: rcode = ' num2str(rcode)])
	 end
         s = ' : Elements';
         for j = 1:sizem
            s = [ s ' ' num2str(temp(j)) ];
         end
      else
         [temp1, rcode] = ncmex('ncvarget', cdfid, rhid, [0], [3], rescale_var);
	 if rcode == -1
	   error(['** ERROR ** ncvarget: rcode = ' num2str(rcode)])
	 end
         s = ' : Elements';
         for j = 1:3
            s = [ s ' ' num2str(temp1(j)) ];
         end
         [temp2, rcode] = ncmex('ncvarget', cdfid, rhid, [sizem-3], [3], rescale_var);
	 if rcode == -1
	   error(['** ERROR ** ncvarget: rcode = ' num2str(rcode)])
	 end
	 s= [ s ' ...' ];
         for j = 1:3
            s = [ s ' ' num2str(temp2(j)) ];
         end
      end
   else
      s = [ ' '];
   end
   s = [ '   ' int2str(i) ')  ' name ' : Length ' int2str(sizem) s ];
   message{i, 1} = s;
   disp(s)
end

% initialise the corner, edge and stride vectors.

if nvdims > 0
  corner = -10*ones(1, nvdims);
  edge = ones(1, nvdims);
  stride = ones(1, nvdims);
else
  corner = 0;
  edge = 1;
  stride = 1;
end

% ask for the index at a point or the corners, edges and strides in
% order to retrieve a (possibly generalised) hyperslab.

take_stride = 0;
for i = 1:nvdims

  % first get the starting point

  name = name_dim{i, 1};
  corner(i) = -1;
  while corner(i) < 0 | corner(i) > ledim(i)
    s = message{i, 1};
    disp(' ')
    disp(s)
    s = [ '    ' name ' : Starting index (between 1 and '];
    s = [ s int2str(ledim(i)+1) ')  (cr for all indices)  ' ];
    clear xtemp;
    xtemp = input(s);
    if isempty(xtemp)
      corner(i) = 0;
      edge(i) = ledim(i) + 1;
      notdone = 0;
    else
      corner(i) = xtemp - 1;
      notdone = 1;
    end
  end

  % next, get the finishing and stride point if these are required.

  if notdone == 1
    end_point = -1;
    ste = [];
    for ii = 1:length(name)
      ste = [ ste ' ' ];
    end
    while end_point < corner(i) | end_point > ledim(i)
      s = [ ste '      finishing index (between ' int2str(corner(i)+1) ];
      s = [ s ' and ' int2str(ledim(i)+1) ')  '];
      ret_val = return_v(s, end_point+1);
      end_point = ret_val - 1;
    end

    stride(i) = -1;
    s=[ ste '      stride length (cr for 1)  ' ];
    while stride(i) < 0 | stride(i) > ledim(i)
      clear xtemp;
      stride(i) = return_v(s, 1);
    end
    
    % Decide whether any non-unit strides are to be taken.

    if stride(i) > 1
      take_stride = 1;
    end

    % Calculate the edge length

    edge(i) = fix( ( end_point - corner(i) )/stride(i) ) + 1;
  end
end

% Retrieve the array.

lenstr = prod(edge);
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
  if nvdims == 0
    [values, rcode] = ncmex('ncvarget1', cdfid, varid, corner, rescale_var);
  else 
    [values, rcode] = ncmex('ncvarget', cdfid, varid, corner, ...
			    edge, rescale_var);
    if rcode == -1
      error(['** ERROR ** ncvarget: rcode = ' num2str(rcode)])
    end
  end
end
  
% Do possible byte correction.
  
if vartypv == nc_byte
  ff = find(values > 127);
  if ~isempty(ff)
    values(ff) = values(ff) - 256;
  end
end

% Handle singleton dimensions.

si = size(values);
len_si = length(si);
squeeze_it = 0;
if (len_si > 2) & (min(si) == 1)
  sq_tmp = 0;
  s = 'Do you want singleton dimensions removed?';
  while any(sq_tmp == [1 2]) == 0
    sq_tmp = menu_old(s, 'yes', 'no');
  end
  squeeze_it = 2 - sq_tmp;
end

% If required do the squeeze.  As well, a new cell array, name_dim_rev, is
% defined to contain the names of the dimensions (whether there has been a
% squeezing or not).

for ii = 1:nvdims
  name_dim_rev{ii} = name_dim{nvdims - ii + 1};
end

if squeeze_it == 1
  ff = find(si ~= 1); % Only non-singleton dimensions are interesting
  for ii = 1:length(ff)
    name_dim_rev{ii} = name_dim_rev{ff(ii)};
  end
  values = squeeze(values);
  si = size(values);
  len_si = length(si);
end

% Calculate num_mults which describes the type of array.

if (len_si > 2)
  num_mults = len_si; % multi-dimensional array.
else
  if max(si) == 1
    num_mults = 0; % number
  elseif min(si) == 1
    num_mults = 1; % vector
  else
    num_mults = 2; % matrix
  end
end

% Manipulate the array according to whether it is a vector, matrix or
% multi-dimensional array.  This may involve permuting arrays.

if num_mults == 0
  
  % getting back a constant
  
elseif num_mults == 1
  
  % ask whether the user wants a row or column vector.

  order = 0;
  s = 'Do you want a row vector or a column vector returned?';
  s1 = [ 'row vector' ];
  s2 = [ 'column vector' ];
  while any(order == [1 2]) == 0
    order = menu_old(s, s1, s2);
  end
  num_rows = size(values, 1);
  if num_rows == 1
    if order == 2
      values = values';
    end
  else
    if order == 1
      values = values';
    end
  end
    
elseif num_mults == 2
  
  % Ask about transposing the matrix. Note that ncmex has returned the
  % elements in the most efficient way, i.e., it has not done any
  % permutation.
  
  order = 0;
  s = 'In which order do you want the indices?';
  s1 = [ varnam '(' name_dim_rev{2} ',' name_dim_rev{1} ')' ];
  s2 = [ varnam '(' name_dim_rev{1} ',' name_dim_rev{2} ')' ];
  while any(order == [1 2]) == 0
    order = menu_old(s, s1, s2);
  end
  if order == 1
    values = values';
  end

else
  
  % A multi-dimensional array. Permute the indices so that they will be
  % consistent with the ncdump output and print out information about
  % multi-dimensional array.
  
  values = permute(values, (num_mults:-1:1));
  
  str = ['The array is ' varnam '('];
  for ii = num_mults:-1:2
    str = [str name_dim_rev{ii} ', '];
  end
  str = [str name_dim_rev{1} ')'];
  disp(str)
end

% Find any scale factors or offsets.

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
else
   scalef = [];
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
else
   addoff = [];
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,