getnc_s.m

来自「matlacb程序包」· M 代码 · 共 2,160 行 · 第 1/5 页

    % 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] = mexnc('ncattinq', cdfid, varid, 'valid_range');
      if rcode < 0
	mess_str = ['ncattinq: ' full_name ': varid =  ' varid ': valid_range'];
	values = error_handle(cdfid, mess_str, rcode, err_opt);
	return
      end
      [miss, rcode] = mexnc('ncattget', cdfid, varid, 'valid_range');
      if rcode < 0
	mess_str = ['ncattget: ' full_name ': varid =  ' varid ': valid_range'];
	values = error_handle(cdfid, mess_str, rcode, err_opt);
	return
      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;
	  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] = mexnc('ncattinq', cdfid, varid, 'valid_min');
	if rcode < 0
	  mess_str = ['ncattinq: ' full_name ': varid =  ' varid ': vvalid_min'];
	  values = error_handle(cdfid, mess_str, rcode, err_opt);
	  return
	end
	[miss_low, rcode] = mexnc('ncattget', cdfid, varid, 'valid_min');
	if rcode < 0
	  mess_str = ['ncattget: ' full_name ': varid =  ' varid ': valid_min'];
	  values = error_handle(cdfid, mess_str, rcode, err_opt);
	  return
	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] = mexnc('ncattinq', cdfid, varid, 'valid_max');
	if rcode < 0
	  mess_str = ['ncattinq: ' full_name ': varid =  ' varid ': vvalid_max'];
	  values = error_handle(cdfid, mess_str, rcode, err_opt);
	  return
	end
	[miss_up, rcode] = mexnc('ncattget', cdfid, varid, 'valid_max');
	if rcode < 0
	  mess_str = ['ncattget: ' full_name ': varid =  ' varid ': valid_max'];
	  values = error_handle(cdfid, mess_str, rcode, err_opt);
	  return
	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] = mexnc('ncattinq', cdfid, varid, '_FillValue');
      if rcode < 0
	mess_str = ['ncattinq: ' full_name ': varid =  ' varid ': _FillValue'];
	values = error_handle(cdfid, mess_str, rcode, err_opt);
	return
      end
      [miss_val, rcode] = mexnc('ncattget', cdfid, varid, '_FillValue');
      if rcode < 0
	mess_str = ['ncattget: ' full_name ': varid =  ' varid ': _FillValue'];
	values = error_handle(cdfid, mess_str, rcode, err_opt);
	return
      end

      % Check that _FillValue 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;
      
      % Check whether _FillValue is outside the valid range to decide
      % whether to keep going.
      
      keep_going = 1;
      if ~isempty(miss_low_orig)
	if (miss_val < miss_low_orig )
	  keep_going = 0;
	end
      end
      if ~isempty(miss_up_orig)
	if (miss_val > miss_up_orig )
	  keep_going = 0;
	end
      end
      
      if keep_going == 1
	
	% 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
    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.  This is only done if the missing_value exists and is
    % different to the _FillValue
    
    pos_missv = check_st('missing_value', attstring, nvatts);
    if pos_missv > 0
      [attype, attlen, rcode] = mexnc('ncattinq', cdfid, varid, 'missing_value');
      if rcode < 0
	mess_str = ['ncattinq: ' full_name ': varid =  ' varid ': missing_value'];
	values = error_handle(cdfid, mess_str, rcode, err_opt);
	return
      end
      [miss_val, rcode] = mexnc('ncattget', cdfid, varid, 'missing_value');
      if rcode < 0
	mess_str = ['ncattget: ' full_name ': varid =  ' varid ': missing_value'];
	values = error_handle(cdfid, mess_str, rcode, err_opt);
	return
      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
      
      % Check whether missing_value is outside the valid range to decide
      % whether to keep going.  Also check whether it equals the original
      % _FillValue.
      
      keep_going = 1;
      if ~isempty(miss_low_orig)
	if (miss_val < miss_low_orig)
	  keep_going = 0;
	end
      end
      if ~isempty(miss_up_orig)
	if (miss_val > miss_up_orig)
	  keep_going = 0;
	end
      end
      if ~isempty(fill_value_orig)
	if (miss_val == fill_value_orig)
	  keep_going = 0;
	end
      end
      
      if keep_going == 1
	
	% 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
    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 change them to a
    % more convenient value.
    
    if len_index_miss > 0
      if change_miss == 2
	if vartypv == nc_char
	  values(index_miss) = char(0);
	else
	  values(index_miss) = NaN;
	end
      elseif change_miss == 3
	if vartypv == nc_char
	  if isnumeric(new_miss)
	    values(index_miss) = char(new_miss);
	  else
	    values(index_miss) = new_miss;
	  end
	else
	  values(index_miss) = new_miss;
	end
      else
	s = [ 'getnc_s was passed change_miss = ' int2str(change_miss) ];
	error(s)
      end
    end
  end

  % Rescale the byte type data which was not done automatically. If the option
  % to not rescale has been selected then scalef and addoff will be empty and
  % there 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] = mexnc('ncclose', cdfid);
  if rcode < 0