/[MITgcm]/MITgcm/utils/matlab/rdmnc.m
ViewVC logotype

Diff of /MITgcm/utils/matlab/rdmnc.m

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.19 by mlosch, Tue Nov 24 13:09:16 2009 UTC revision 1.20 by mlosch, Fri May 14 11:29:16 2010 UTC
# Line 41  files={}; Line 41  files={};
41  varlist={};  varlist={};
42  iters=[];  iters=[];
43    
44    % find out if matlab's generic netcdf API is available
45    % if not assume that Chuck Denham's netcdf toolbox is installed
46    usemexcdf = isempty(which('netcdf.open'));
47    
48  % Process function arguments  % Process function arguments
49  for iarg=1:nargin;  for iarg=1:nargin;
50      arg=varargin{iarg};      arg=varargin{iarg};
# Line 95  if isempty(iters) Line 99  if isempty(iters)
99      iters = [];      iters = [];
100      for ieachfile=1:length(files)      for ieachfile=1:length(files)
101          eachfile = [filepaths{ieachfile},files{ieachfile}];          eachfile = [filepaths{ieachfile},files{ieachfile}];
102          nc=netcdf(char(eachfile),'read');          if usemexcdf
103          nciters = nc{'iter'}(:);            nc=netcdf(char(eachfile),'read');
104          if isempty(nciters), nciters = nc{'T'}(:); end            nciters = nc{'iter'}(:);
105              if isempty(nciters), nciters = nc{'T'}(:); end
106            else
107              nc=netcdf.open(char(eachfile),'NC_NOWRITE');
108              nciters = ncgetvar(nc,'iter');
109              if isempty(nciters), nciters = ncgetvar(nc,'T'); end
110            end
111          iters = [iters,nciters'];          iters = [iters,nciters'];
112          close(nc);          if usemexcdf
113              close(nc);
114            else
115              netcdf.close(nc);
116            end
117      end      end
118      iters = unique(iters');      iters = unique(iters');
119  end  end
# Line 109  S.attributes=[]; Line 123  S.attributes=[];
123  for ieachfile=1:length(files)  for ieachfile=1:length(files)
124      eachfile = [filepaths{ieachfile},files{ieachfile}];      eachfile = [filepaths{ieachfile},files{ieachfile}];
125      if dBug > 0, fprintf([' open: ',eachfile]); end      if dBug > 0, fprintf([' open: ',eachfile]); end
126      nc=netcdf(char(eachfile),'read');      if usemexcdf
127      S=rdmnc_local(nc,varlist,iters,S,dBug);        nc=netcdf(char(eachfile),'read');
128      close(nc);        S=rdmnc_local(nc,varlist,iters,S,dBug);
129          close(nc);
130        else
131          nc=netcdf.open(char(eachfile),'NC_NOWRITE');
132          S=rdmnc_local_matlabAPI(char(eachfile),nc,varlist,iters,S,dBug);
133          netcdf.close(nc);
134        end
135  end  end
136    
137    return
138    
139  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
140  %                             Local functions                             %  %                             Local functions                             %
# Line 262  function [S] = rdmnc_local(nc,varlist,it Line 283  function [S] = rdmnc_local(nc,varlist,it
283    end    end
284        
285    return    return
286    
287    function [S] = rdmnc_local_matlabAPI(fname,nc,varlist,iters,S,dBug)
288    
289      fiter = ncgetvar(nc,'iter');                           % File iterations present
290      if isempty(fiter), fiter = ncgetvar(nc,'T'); end
291      if isinf(iters); iters = fiter(end); end
292      if isnan(iters); iters = fiter; end
293      [fii,dii] = ismember(fiter,iters);  fii = find(fii); % File iteration index
294      dii = dii(find(dii ~= 0));                           % Data interation index
295      if dBug > 0,
296        fprintf(' ; fii='); fprintf(' %i',fii);
297        fprintf(' ; dii='); fprintf(' %i',dii); fprintf(' \n');
298      end
299    
300      % get basic information about netcdf file
301      [ndims nvars natts dimm] = netcdf.inq(nc);
302    
303      % No variables specified? Default to all
304      if isempty(varlist),
305        for k=0:nvars-1
306          varlist{k+1} = netcdf.inqVar(nc,k);
307        end
308      end
309      
310      % Attributes for structure
311      if iters>0; S.iters_from_file=iters; end
312      S.attributes.global=ncgetatt(nc,'global');
313      [pstr,netcdf_fname,ext] = fileparts(fname);
314      if strcmp(netcdf_fname(end-3:end),'glob')
315        % assume it is a global file produced by gluemnc and change some
316        % attributes
317        S.attributes.global.sNx = S.attributes.global.Nx;
318        S.attributes.global.sNy = S.attributes.global.Ny;
319        S.attributes.global.nPx = 1;
320        S.attributes.global.nSx = 1;
321        S.attributes.global.nPy = 1;
322        S.attributes.global.nSy = 1;
323        S.attributes.global.tile_number = 1;
324        S.attributes.global.nco_openmp_thread_number = 1;
325      end
326        
327      % Read variable data
328      for ivar=1:size(varlist,2)
329        
330        cvar=char(varlist{ivar});
331        varid=ncfindvarid(nc,cvar);
332        if isempty(varid)
333          disp(['No such variable ''',cvar,''' in MNC file ',fname]);
334          continue
335        end
336        
337        [varname,xtype,dimids,natts] = netcdf.inqVar(nc,varid);
338        
339        % Dimensions
340        clear sizVar dims
341        for k=1:length(dimids)
342          [dims{k}, sizVar(k)] = netcdf.inqDim(nc,dimids(k));
343        end
344        if length(sizVar)==1; sizVar(2) = 1; end;
345        nDims=length(sizVar);
346        if dims{1} == 'T'
347          if isempty(find(fii)), error('Iters not found'); end
348          it = length(dims);
349          tmpdata = netcdf.getVar(nc,varid,'double');
350          tmpdata = tmpdata(fii,:);
351          % leading unity dimensions get lost; add them back:
352          tmpdata=reshape(tmpdata,[length(fii) sizVar(2:end)]);
353        else
354          it = 0;
355          tmpdata = netcdf.getVar(nc,varid,'double');
356          % leading unity dimensions get lost; add them back:
357          tmpdata=reshape(tmpdata,sizVar);
358        end
359        
360        if dBug > 1,
361          fprintf(['  var:',cvar,': nDims=%i ('],nDims);fprintf(' %i',sizVar);
362          fprintf('):%i,nD=%i,it=%i ;',length(size(tmpdata)),length(dims),it);
363        end
364        [ni nj nk nm nn no np]=size(tmpdata);
365          
366        [i0,j0,fn]=findTileOffset(S);
367        cdim=dims{1}; if cdim(1)~='X'; i0=0; end
368        cdim=dims{1}; if cdim(1)=='Y'; i0=j0; j0=0; end
369        if length(dims)>1;
370          cdim=dims{2}; if cdim(1)~='Y'; j0=0; end
371        else
372          j0=0;
373        end
374        if dBug > 1,
375          fprintf(' i0,ni= %i %i; j,nj= %i %i; nk=%i :',i0,ni,j0,nj,nk);
376        end
377        
378        Sstr = '';
379        for istr = 1:max(nDims,length(dims)),
380          if     istr == it,  Sstr = [Sstr,'dii,'];
381          elseif istr == 1,   Sstr = [Sstr,'i0+(1:ni),'];
382          elseif istr == 2,   Sstr = [Sstr,'j0+(1:nj),'];
383          elseif istr == 3,   Sstr = [Sstr,'(1:nk),'];
384          elseif istr == 4,   Sstr = [Sstr,'(1:nm),'];
385          elseif istr == 5,   Sstr = [Sstr,'(1:nn),'];
386          elseif istr == 6,   Sstr = [Sstr,'(1:no),'];
387          elseif istr == 7,   Sstr = [Sstr,'(1:np),'];
388          else, error('Can''t handle this many dimensions!');
389          end
390        end
391        eval(['S.(cvar)(',Sstr(1:end-1),')=tmpdata;'])
392        %S.(cvar)(i0+(1:ni),j0+(1:nj),(1:nk),(1:nm),(1:nn),(1:no),(1:np))=tmpdata;
393        if dBug > 1, fprintf(' %i',size(S.(cvar))); fprintf('\n'); end
394        
395        S.attributes.(cvar)=ncgetatt(nc,cvar);
396        % replace missing or FillValues with NaN
397        if isfield(S.attributes.(cvar),'missing_value');
398          misval = S.attributes.(cvar).missing_value;
399          S.(cvar)(S.(cvar) == misval) = NaN;
400        end
401        if isfield(S.attributes.(cvar),'FillValue_');
402          misval = S.attributes.(cvar).FillValue_;
403          S.(cvar)(S.(cvar) == misval) = NaN;
404        end
405        
406      end % for ivar
407        
408      if isempty(S)
409        error('Something didn''t work!!!');
410      end
411      
412      return
413    
414    function vf = ncgetvar(nc,varname)
415    % read a netcdf variable
416      
417    % find out basics about the files
418      [ndims nvars natts dimm] = netcdf.inq(nc);
419      vf = [];
420      varid = [];
421      for k=0:nvars-1
422        if strcmp(netcdf.inqVar(nc,k),varname)
423          varid = netcdf.inqVarId(nc,varname);
424        end
425      end
426      if ~isempty(varid);
427        [varn,xtype,dimids,natts] = netcdf.inqVar(nc,varid);
428        % get data
429        vf = double(netcdf.getVar(nc,varid));
430      else
431        % do nothing
432      end
433    
434      return
435    
436    function misval = ncgetmisval(nc,varid)
437    
438      [varname,xtype,dimids,natts] = netcdf.inqVar(nc,varid);
439      misval = [];
440      for k=0:natts-1
441        attn = netcdf.inqAttName(nc,varid,k);
442        if strcmp(attn,'missing_value') | strcmp(attn,'_FillValue')
443          misval = double(netcdf.getAtt(nc,varid,attname));
444        end
445      end
446      
447    function A = ncgetatt(nc,varname)
448    % get all attributes and put them into a struct
449      
450    % 1. get global properties of file
451      [ndims nvars natts dimm] = netcdf.inq(nc);
452    
453      % get variable ID and properties
454      if strcmp('global',varname)
455        % "variable" is global
456        varid = netcdf.getConstant('NC_GLOBAL');
457      else
458        % find variable ID and properties
459        varid = ncfindvarid(nc,varname);
460        if ~isempty(varid)
461          [varn,xtype,dimids,natts] = netcdf.inqVar(nc,varid);
462        else
463          warning(sprintf('variable %s not found',varname))
464        end
465      end
466    
467      if natts > 1
468        for k=0:natts-1
469          attn = netcdf.inqAttName(nc,varid,k);
470          [xtype attlen]=netcdf.inqAtt(nc,varid,attn);
471          attval = netcdf.getAtt(nc,varid,attn);
472          if ~ischar(attval)
473            attval = double(attval);
474          end
475          A.(char(attn))=attval;
476        end
477      else
478          A = 'none';
479      end
480      
481      return
482    
483      
484    function varid = ncfindvarid(nc,varname)
485    
486      [ndims nvars natts dimm] = netcdf.inq(nc);
487      varid=[];
488      for k=0:nvars-1
489        if strcmp(netcdf.inqVar(nc,k),varname);
490          varid = netcdf.inqVarId(nc,varname);
491          break
492        end
493      end
494    
495      return
Legend:
Removed from v.1.19  
changed lines
  Added in v.1.20
  ViewVC Help
Powered by ViewVC 1.1.22