getnc.m

读取Network Common Data Form (netCDF)数据

% index_missing_value and index__FillValue.

index_miss_low = [];
index_miss_up = [];
index__FillValue = [];
index_missing_value = [];

% 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);
   
   % 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;
       miss_up = miss_up + addoff;
     end
   end
   
   index_miss_low = find ( values < miss_low );
   index_miss_up = find ( values > miss_up );
 
else
  pos_min = check_st('valid_min', attstring, nvatts);
  if pos_min > 0
    [attype, attlen, rcode] = ncmex('ncattinq', cdfid, varid, 'valid_min');
    if rcode == -1
      error(['** ERROR ** ncattinq: rcode = ' num2str(rcode)])
    end
    [miss_low, rcode] = ncmex('ncattget', cdfid, varid, 'valid_min');
    if rcode == -1
      error(['** ERROR ** ncattget: rcode = ' num2str(rcode)])
    end

    % Check that valid_min is a scalar
    
    if length(miss_low) ~= 1
      error(['The valid_min attribute must be a scalar'])
    end
    
    % Correct for possible faulty handling of byte type
    
    if attype == nc_byte
      if miss_low > 127; miss_low = miss_low - 256; end
    end
    miss_low_orig = miss_low;
    
    % Rescale & add offsets if required.
   
    if rescale_att == 1
      if isempty(scalef) == 0
	miss_low = miss_low*scalef;
      end
      if isempty(addoff) == 0
	miss_low = miss_low + addoff;
      end
    end
      
    index_miss_low = find ( values < miss_low );
  end

  pos_max = check_st('valid_max', attstring, nvatts);
  if pos_max > 0
    [attype, attlen, rcode] = ncmex('ncattinq', cdfid, varid, 'valid_max');
    if rcode == -1
      error(['** ERROR ** ncattinq: rcode = ' num2str(rcode)])
    end
    [miss_up, rcode] = ncmex('ncattget', cdfid, varid, 'valid_max');
    if rcode == -1
      error(['** ERROR ** ncattget: rcode = ' num2str(rcode)])
    end

    % Check that valid_max is a scalar
    
    if length(miss_up) ~= 1
      error(['The valid_max attribute must be a scalar'])
    end
    % Correct for possible faulty handling of byte type
    
    if attype == nc_byte
      if miss_up > 127; miss_up = miss_up - 256; end
    end
    miss_up_orig = miss_up;

    % Rescale & add offsets if required.
   
    if rescale_att == 1
      if isempty(scalef) == 0
	miss_up = miss_up*scalef;
      end
      if isempty(addoff) == 0
	miss_up = miss_up + addoff;
      end
    end
    
    index_miss_up = find ( values > miss_up );
  end
end

% Now find the indices of the data points that are 'close to'
% _FillValue.  Note that 'close to' is different according to the
% data type.

pos_missv = check_st('_FillValue', attstring, nvatts);
if pos_missv > 0
   [attype, attlen, rcode] = ncmex('ncattinq', cdfid, varid, '_FillValue');
   if rcode == -1
     error(['** ERROR ** ncattinq: rcode = ' num2str(rcode)])
   end
   [miss_val, rcode] = ncmex('ncattget', cdfid, varid, '_FillValue');
   if rcode == -1
     error(['** ERROR ** ncattget: rcode = ' num2str(rcode)])
   end

    % Check that _FillValue_orig is a scalar
    
    if length(miss_val) ~= 1
      error(['The _FillValue attribute must be a scalar'])
    end
    
    % Correct for possible faulty handling of byte type
    
    if attype == nc_byte
      if miss_val > 127; miss_val = miss_val - 256; end
    end
    fill_value_orig = miss_val;
   
   % Rescale & add offsets if required.
   
   if rescale_att == 1
     if isempty(scalef) == 0
       miss_val = miss_val*scalef;
     end
     if isempty(addoff) == 0
       miss_val = miss_val + addoff;
     end
   end
   
   if attype == nc_byte | attype == nc_char
     index__FillValue = find ( values == miss_val );
   elseif attype == nc_short | attype == nc_long
     need_index_m = 1;
     if pos_vr > 0 | pos_min > 0
       if miss_val < miss_low
	 need_index_m = 0;
       end
     end
     if pos_vr > 0 | pos_max > 0
       if miss_val > miss_up
	 need_index_m = 0;
       end
     end
     if need_index_m
       index__FillValue = find ( values == miss_val );
     end
   elseif attype == nc_float | attype == nc_double
     need_index_m = 1;
     if miss_val < 0
       miss_val_low = 1.00001*miss_val;
       miss_val_up = 0.99999*miss_val;
     else
       miss_val_low = 0.99999*miss_val;
       miss_val_up = 1.00001*miss_val;
     end
     
     if pos_vr > 0 | pos_min > 0
       if miss_val_up < miss_low
	 need_index_m = 0;
       end
     end
     if pos_vr > 0 | pos_max > 0
       if miss_val_low > miss_up
	 need_index_m = 0;
       end
     end
     if need_index_m
       index__FillValue = find ( miss_val_low <= values & ...
	   values <= miss_val_up );
     end
   end
 end

% Now find the indices of the data points that are 'close to'
% missing_value.  Note that 'close to' is different according to the
% data type.

pos_missv = check_st('missing_value', attstring, nvatts);
if pos_missv > 0
  [attype, attlen, rcode] = ncmex('ncattinq', cdfid, varid, 'missing_value');
  if rcode == -1
    error(['** ERROR ** ncattinq: rcode = ' num2str(rcode)])
  end
  [miss_val, rcode] = ncmex('ncattget', cdfid, varid, 'missing_value');
  if rcode == -1
    error(['** ERROR ** ncattget: rcode = ' num2str(rcode)])
  end
   
  % Check that missing_value is a scalar
    
  if length(miss_val) ~= 1
    error(['The missing_value attribute must be a scalar'])
  end
    
  % Correct for possible faulty handling of byte type
    
  if attype == nc_byte
    if miss_val > 127; miss_val = miss_val - 256; end
  end
    
  % Rescale & add offsets if required.
   
  if rescale_att == 1
    if isempty(scalef) == 0
      miss_val = miss_val*scalef;
    end
    if isempty(addoff) == 0
      miss_val = miss_val + addoff;
    end
  end
   
  if attype == nc_byte | attype == nc_char
    index_missing_value = find ( values == miss_val );
  elseif attype == nc_short | attype == nc_long
    need_index_m = 1;
    if pos_vr > 0 | pos_min > 0
      if miss_val < miss_low
	need_index_m = 0;
      end
    end
    if pos_vr > 0 | pos_max > 0
      if miss_val > miss_up
	need_index_m = 0;
      end
    end
    if need_index_m
      index_missing_value = find ( values == miss_val );
    end
  elseif attype == nc_float | attype == nc_double
    need_index_m = 1;
    if miss_val < 0
      miss_val_low = 1.00001*miss_val;
      miss_val_up = 0.99999*miss_val;
    else
      miss_val_low = 0.99999*miss_val;
      miss_val_up = 1.00001*miss_val;
    end

    if pos_vr > 0 | pos_min > 0
      if miss_val_up < miss_low
	need_index_m = 0;
      end
    end
    if pos_vr > 0 | pos_max > 0
      if miss_val_low > miss_up
	need_index_m = 0;
      end
    end
    if need_index_m
      index_missing_value = find ( miss_val_low <= values & ...
	  values <= miss_val_up );
    end
  end
end

%Combine the arrays of missing value indices into one unordered array.
%Note that for real numbers the range of the _FillValue and
%missing_value may intersect both the valid and invalid range and so
%some indices may appear twice; this does not cause any inaccuracy,
%although it will result in some inefficiency.  In particular, rescaling
%is done on the set of indices NOT in index_miss and so is not
%affected.

index_miss = [ index_miss_low(:); index__FillValue(:); ...
    index_missing_value(:); index_miss_up(:) ];
%index_miss = sort(index_miss);
len_index_miss = length(index_miss);

% If there are any missing values then offer to change them to a
% more convenient value.

if len_index_miss > 0
   s = [ varnam ' contains missing values:  Choose an action' ];
   s1 = 'Leave the missing value unchanged';
   s2 = 'Replace the missing value with NaN';
   s3 = 'Replace the missing value with a new value';
   k = -1;
   while any(k == [1 2 3]) == 0
      k = menu_old(s, s1, s2, s3);
      if k == 1
      elseif k == 2
	  values(index_miss) = NaN*ones(size(index_miss));
	  if vartypv == nc_char
	    values = setstr(values);
	  end
      elseif k == 3
	 if vartypv == nc_char
	    s = '   Type in your new missing value marker [*]  ';
	    new_miss = return_v(s, '*');
	    values(index_miss) = new_miss*ones(size(index_miss));
	    values = setstr(values);
	  else
	    s = '   Type in your new missing value marker [0]  ';
	    new_miss = return_v(s, 0);
	    values(index_miss) = new_miss*ones(size(index_miss));
	 end
      else
         disp(' ')
         disp('You have asked for a non-existent option - try again')
      end
   end
end

% Rescale the byte type data which was not done automatically. If the otion
% to not rescale has been selected then scalef and addoff will be empty and
% ther will be no rescaling.

if vartypv == nc_byte
  if isempty(scalef) == 0
    values = values*scalef;
  end
  if isempty(addoff) == 0
    values = values + addoff;
  end
end
    
% Close the netcdf file.

[rcode] = ncmex('ncclose', cdfid);
if rcode == -1
  error(['** ERROR ** ncclose: rcode = ' num2str(rcode)])
end