timenc.m

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

      if bl_corner(ii) < 0
	bl_corner(ii) = 1;
      elseif bl_corner(ii) > sizem(ii)
	error(['bl_corner so big that it is outside the hyperslab'])
      end
      if tr_corner(ii) < 0
	tr_corner(ii) = sizem(ii);
      elseif tr_corner(ii) > sizem(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      
      str_h = [str_h '[' num2str(bl_corner(ii) - 1) ':' ...
	       num2str(tr_corner(ii) - 1) ']'];
    end
    
    switch method_of_call
     case 1
      switch mex_name
       case 'loaddap'
	values_struct = loaddap('+v', [full_name '?' time_var str_h]);
       case 'loaddods'
	values_struct = loaddods('+v', [full_name '?' time_var str_h]);
      end
     case 2
      switch mex_name
       case 'loaddap'
	loaddap('+v', [full_name '?' time_var str_h]);
       case 'loaddods'
	loaddods('+v', [full_name '?' time_var str_h]);
      end
      eval(['values_struct = ' time_var ';'])
    end
  end
  
  % Allow for the various weird ways that the data may be returned.
  
  if isstruct(values_struct)
    %xx = values_struct.(time_var);
    xx = getfield(values_struct, time_var);
    if isstruct(xx)
      %serial_rel = xx.(time_var);
      serial_rel = getfield(xx, time_var);
    else
      serial_rel = xx;
    end
  else
    serial_rel = values_struct;
  end

  % Allow for serial_rel to be multi-dimensional in which case we have to
  % work around the weird way that loaddap returns multi-dimensional
  % data. Then make it into a column vector for backwards compatibility with
  % the code that Rose put in for a multi-dimensional time variable in netcdf
  % files.
  
  if num_dim > 2
    vec_permute = [num_dim:-1:3 1 2];
    serial_rel = permute(serial_rel, vec_permute);
  end
  serial_rel = serial_rel(:)';
 case 'java'
  
  % Get the file object
  
  try
    ncdJ = ucar.nc2.dataset.NetcdfDataset.openDataset(file);
  catch
    ss = lasterror;
    mess_str = ['Failed trying to open file: ' file ss.message];      
    rcode = -1000000;
    [gregorian_time, serial_time, gregorian_base, serial_base, ...
     sizem, serial_time_jd, serial_base_jd] = error_handle('java', ...
						  mess_str, rcode, err_opt);
    return
  end
  
  % Get the variable object
  
  try
    varJ = ncdJ.findVariable(time_var);
  catch
    ss = lasterror;
    mess_str = ['Failed getting ' time_var ' in file: ' file ss.message];      
    rcode = -1000000;
    [gregorian_time, serial_time, gregorian_base, serial_base, ...
     sizem, serial_time_jd, serial_base_jd] = error_handle('java', ...
						  mess_str, rcode, err_opt);
    return
  end
  
  % Get the actual data as serial_rel.
  
  sizem = varJ.getShape();
  if any(sizem == 0)
    serial_rel = []; % Presumably time is an unlimited dimension of length 0.
    disp(['Warning: There are apparently no ''time'' records'])
  else
    if get_all == 1
      serial_rel = squeeze(copyToNDJavaArray(varJ.read()));
    else
      % Check that the bl_corner points are acceptable.
      
      varRank = varJ.getRank;
      for ii = 1:varRank
	if (bl_corner(ii) < 1) || (bl_corner(ii) > tr_corner(ii)) || ...
	      (tr_corner(ii) > sizem(ii))
	  mess_str = 'hyperslab is badly specified';
	  rcode = -1000000;
	  [gregorian_time, serial_time, gregorian_base, serial_base, ...
	   sizem, serial_time_jd, serial_base_jd] = error_handle('java', ...
						  mess_str, rcode, err_opt);
	  return
	end
      end
      
      % Construct the string that specifies the hyperslab and then get the data.
      readSpec = '';
      for ii = 1:varRank
	readSpec = [readSpec num2str(bl_corner(ii) - 1) ':' ...
		    num2str(tr_corner(ii) - 1) ':1,'];
      end
      readSpec = readSpec(1:(end - 1));
      serial_rel = squeeze(copyToNDJavaArray(varJ.read(readSpec)));
    end
  end
    
  % Get the string describing the base date.
  
  try
    att = varJ.findAttribute('units');
    base_str = char(att.getStringValue());
  catch
    ss = lasterror;
    mess_str = ['Failed getting units attribute in file: ' file ss.message];
    rcode = -1000000;
    [gregorian_time, serial_time, gregorian_base, serial_base, ...
     sizem, serial_time_jd, serial_base_jd] = error_handle('java', ...
						  mess_str, rcode, err_opt);
    return
  end  

  % Close the file object
  
  ncdJ.close();
 case 'none'
  error(['Couldn''t find a suitable mex-file for reading ' file])
end

% Find out what calendar we are using. If there is no calendar attribute then
% set it to 'gregorian';

if isempty(calendar)
  [att_val, att_name_list] = attnc(full_name, time_var);
  calendar = 'gregorian';
  for ii = 1:length(att_name_list)
    if strcmp(lower(att_name_list{ii}), 'calendar')
      calendar = att_val{ii};
      break;
    end
  end
end

% Parse the string containing the base date to get its constituents and
% then find its serial and gregorian dates. Also rescale the relative serial
% time vector to turn it into days since the base time.

[gregorian_base, rescale_serial_rel, serial_base_jd, serial_base] = ...
    parsetnc(base_str);
if rescale_serial_rel ~= 1
  serial_rel = rescale_serial_rel*serial_rel;
end

% Find the absolute serial date and resultant gregorian date of the time
% vector.

serial_time_jd = serial_rel + serial_base_jd;
if isempty(serial_time_jd)
  gregorian_time = [];
  serial_time = [];
else
  switch lower(calendar)
   case {'standard', 'gregorian'}
    gregorian_time = get_calendar_date(serial_time_jd);
    serial_time = datenum(gregorian_time(:, 1), gregorian_time(:, 2), ...
			  gregorian_time(:, 3), gregorian_time(:, 4), ...
			  gregorian_time(:, 5), gregorian_time(:, 6));

   case 'proleptic_gregorian'
    serial_time = serial_rel + serial_base;
    serial_time = serial_time(:);
    gregorian_time = datevec(serial_time);
   case {'noleap', '365_day'}
    % We use serial_base to give us a proper starting time and work from
    % there in steps of 365 days per year.
    days_per_month = [31 28 31 30 31 30 31 31 30 31 30 31];
    days_ref = [0 cumsum(days_per_month)];
    [year_b, month_b, day_b, hour_b, minute_b, sec_b] = datevec(serial_base);
    days_from_year_base = days_ref(month_b) + day_b - 1 + hour_b/24 + ...
       minute_b/1440 + sec_b/86400;
    day_full = serial_rel + days_from_year_base;
    year_rel = floor(day_full/365);
    year_abs = year_b + year_rel;
    rem_1 = day_full - year_rel*365;
    month_abs = zeros(1, sizem);
    for ii = 1:sizem
       ff = find(days_ref <= rem_1(ii));
       month_abs(ii) = ff(end);
       rem_2(ii) = rem_1(ii) - days_ref(month_abs(ii));
    end
    day_rel = floor(rem_2);
    day_abs = day_rel + 1;
    rem_3 = (rem_2 - day_rel)*24;
    hour_abs = floor(rem_3);
    rem_4 = (rem_3 - hour_abs)*60;
    minute_abs = floor(rem_4);
    second_abs = (rem_4 - minute_abs)*60;
    gregorian_time = [year_abs(:) month_abs(:) day_abs(:) hour_abs(:) ...
       minute_abs(:) second_abs(:)];
    serial_time = datenum(gregorian_time);
   case {'all_leap', '366_day'}
    % We use serial_base to give us a proper starting time and work from
    % there in steps of 366 days per year.
    days_per_month = [31 29 31 30 31 30 31 31 30 31 30 31];
    days_ref = [0 cumsum(days_per_month)];
    [year_b, month_b, day_b, hour_b, minute_b, sec_b] = datevec(serial_base);
    days_from_year_base = days_ref(month_b) + day_b - 1 + hour_b/24 + ...
       minute_b/1440 + sec_b/86400;
    day_full = serial_rel + days_from_year_base;
    year_rel = floor(day_full/366);
    year_abs = year_b + year_rel;
    rem_1 = day_full - year_rel*366;
    month_abs = zeros(1, sizem);
    for ii = 1:sizem
       ff = find(days_ref <= rem_1(ii));
       month_abs(ii) = ff(end);
       rem_2(ii) = rem_1(ii) - days_ref(month_abs(ii));
    end
    day_rel = floor(rem_2);
    day_abs = day_rel + 1;
    rem_3 = (rem_2 - day_rel)*24;
    hour_abs = floor(rem_3);
    rem_4 = (rem_3 - hour_abs)*60;
    minute_abs = floor(rem_4);
    second_abs = (rem_4 - minute_abs)*60;
    gregorian_time = [year_abs(:) month_abs(:) day_abs(:) hour_abs(:) ...
       minute_abs(:) second_abs(:)];
    serial_time = datenum(gregorian_time);
   case '360_day'
    % We use serial_base to give us a proper starting time and work from
    % there in steps of 366 days per year.
    [year_b, month_b, day_b, hour_b, minute_b, sec_b] = datevec(serial_base);
    days_from_year_base = 30*(month_b - 1) + day_b - 1 + hour_b/24 + ...
       minute_b/1440 + sec_b/86400;
    day_full = serial_rel + days_from_year_base;
    year_rel = floor(day_full/360);
    year_abs = year_b + year_rel;
    rem_1 = day_full - year_rel*360;
    month_abs = floor(rem_1/30) + 1;
    rem_2 = rem_1 - (month_abs - 1)*30;
    day_rel = floor(rem_2);
    day_abs = day_rel + 1;
    rem_3 = (rem_2 - day_rel)*24;
    hour_abs = floor(rem_3);
    rem_4 = (rem_3 - hour_abs)*60;
    minute_abs = floor(rem_4);
    second_abs = (rem_4 - minute_abs)*60;
    gregorian_time = [year_abs(:) month_abs(:) day_abs(:) hour_abs(:) ...
       minute_abs(:) second_abs(:)];
    serial_time = datenum(gregorian_time);
   otherwise
    disp(['!! timenc cannot handle the calendar attribute **' calendar '**'])
    disp('!! which may have been found in the original netCDF file. The help')
    disp('!! message for timenc tells you how to specify a different calendar')
    error('strange calendar')
  end
end

function [gregorian_time, serial_time, gregorian_base, serial_base, ...
	  sizem, serial_time_jd, serial_base_jd] = error_handle(fid, ...
						  mess_str, rcode, err_opt)

% error_handle is called after a mexnc or java call has failed. It ensures
% that an open netcdf file is closed. The value of err_opt determines what
% else is done. For a mexnc call fid is cdfid, the handle to the open
% file. For a java call fid is the opened file object.
%    err_opt == 1 prints an error message and then aborts
%            == 2 prints a warning message and then returns an empty
%                 array. This is the default.
%            == 3 returns an empty array. This is a very dangerous option and
%                 should only be used with caution. It might be used when
%                 getnc_s is called in a loop and you want to do your own
%                 error handling without being bothered by warning messages.

% Decide what part of the code made the call.
  
  if isempty(fid)
    called_by = 'loadd';
  else
    if isnumeric(fid)
      called_by = 'mexnc';
    else
      called_by = 'java';
    end
  end

% Close an open netcdf of java file.
  
  switch called_by
   case 'mexnc'
    if fid >= 0
      [rcode_sub] = mexnc('ncclose', fid);
    end
   case 'java'
    if isjava(fid)
      fid.close();
    end
  end
  
  % Handle the errors according to the value of err_opt. If rcode is empty
  % then this is probably because loaddap or loaddods was called.
  
  if ~exist('gregorian_time')
    gregorian_time = [];
  end
  if ~exist('serial_time')
    serial_time = [];
  end
  if ~exist('gregorian_base')
    gregorian_base = [];
  end
  if ~exist('serial_base')
    serial_base = [];
  end
  if ~exist('sizem')
    sizem = [];
  end
  if ~exist('serial_time_jd')
    serial_time_jd = [];
  end
  if ~exist('serial_base_jd')
    serial_base_jd = [];
  end
  
  switch err_opt
   case 1
    if isempty(rcode)
      str = ['ERROR: ' mess_str];
    else
      str = ['ERROR: ' mess_str ' : rcode = ' num2str(rcode)];
    end
    error(str)
   case 2
    if isempty(rcode)
      str = ['WARNING: ' mess_str];
    else
      str = ['WARNING: ' mess_str ' : rcode = ' num2str(rcode)];
    end
    disp(str)
   case 3
    return
   otherwise
    error(['error_handle was called with err_opt = ' num2str(err_opt)])
  end