/[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.7 by enderton, Fri May 12 16:52:09 2006 UTC revision 1.25 by mlosch, Thu Jun 23 08:40:09 2011 UTC
# Line 30  function [S] = rdmnc(varargin) Line 30  function [S] = rdmnc(varargin)
30  %  Author:  Alistair Adcroft  %  Author:  Alistair Adcroft
31  %  Modifications:  Daniel Enderton  %  Modifications:  Daniel Enderton
32    
33    % $Header$
34    % $Name$
35    
36  % Initializations  % Initializations
37    dBug=0;
38  file={};  file={};
39  filepaths={};  filepaths={};
40  files={};  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 51  for iarg=1:nargin; Line 59  for iarg=1:nargin;
59              iters=arg;              iters=arg;
60          else          else
61              error(['The only allowed numeric argument is iterations',...              error(['The only allowed numeric argument is iterations',...
62                     ' to read in as a vector for the last arguement']);                     ' to read in as a vector for the last argument']);
63          end          end
64      end      end
65  end  end
66    if isempty(file)
67        if isempty(varlist),
68           fprintf( 'No file name in argument list\n');
69        else
70           fprintf(['No file in argument list:\n ==> ',char(varlist(1))]);
71           for i=2:size(varlist,2), fprintf([' , ',char(varlist(i))]); end
72           fprintf(' <==\n');
73        end
74        error(' check argument list !!!');
75    end
76    
77  % Create list of filenames  % Create list of filenames
78  for eachfile=file  for eachfile=file
79          filepathtemp=eachfile{:};          filepathtemp=eachfile{:};
80          indecies = find(filepathtemp=='/');          indecies = find(filepathtemp==filesep);
81          if ~isempty(indecies)          if ~isempty(indecies)
82          filepathtemp = filepathtemp(1:indecies(end));          filepathtemp = filepathtemp(1:indecies(end));
83          else          else
# Line 81  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          iters = [iters,nciters];            if isempty(nciters), nciters = nc{'T'}(:); end
106          close(nc);            close(nc);
107            else
108              % the parser complains about "netcdf.open" when the matlab netcdf
109              % API is not available, even when it is not used so we have to
110              % avoid the use of "netcdf.open", etc in this function
111              nciters = ncgetvar(char(eachfile),'iter');
112              if isempty(nciters), nciters = ncgetvar(char(eachfile),'T'); end
113            end
114            iters = [iters,nciters'];
115      end      end
116      iters = unique(iters);      iters = unique(iters');
117  end  end
118    
119  % Cycle through files for data extraction.  % Cycle through files for data extraction.
120  S.attributes=[];  S.attributes=[];
121  for ieachfile=1:length(files)  for ieachfile=1:length(files)
122      eachfile = [filepaths{ieachfile},files{ieachfile}];      eachfile = [filepaths{ieachfile},files{ieachfile}];
123      nc=netcdf(char(eachfile),'read');      if dBug > 0, fprintf([' open: ',eachfile]); end
124      S=rdmnc_local(nc,varlist,iters,S);      if usemexcdf
125      close(nc);        nc=netcdf(char(eachfile),'read');
126          S=rdmnc_local(nc,varlist,iters,S,dBug);
127          close(nc);
128        else
129          % the parser complains about "netcdf.open" when the matlab netcdf
130          % API is not available, even when it is not used so we have to
131          % avoid the use of "netcdf.open", etc in this function
132          S=rdmnc_local_matlabAPI(char(eachfile),varlist,iters,S,dBug);
133        end
134  end  end
135    
136    return
137    
138  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
139  %                             Local functions                             %  %                             Local functions                             %
# Line 131  function [i0,j0,fn] = findTileOffset(S); Line 166  function [i0,j0,fn] = findTileOffset(S);
166      end      end
167      %[snx,sny,nsx,nsy,npx,npy,ntx,nty,i0,j0,fn];      %[snx,sny,nsx,nsy,npx,npy,ntx,nty,i0,j0,fn];
168    
169  function [S] = rdmnc_local(nc,varlist,iters,S)  function [S] = rdmnc_local(nc,varlist,iters,S,dBug)
170    
171      fiter = nc{'iter'}(:);                               % File iterations present    fiter = nc{'iter'}(:);                               % File iterations present
172      if isempty(fiter), fiter = nc{'T'}(:); end    if isempty(fiter), fiter = nc{'T'}(:); end
173      [fii,dii] = ismember(fiter,iters);  fii = find(fii); % File iteration index    if isinf(iters); iters = fiter(end); end
174      dii = dii(find(dii ~= 0));                           % Data interation index    if isnan(iters); iters = fiter; end
175          [fii,dii] = ismember(fiter,iters);  fii = find(fii); % File iteration index
176      % No variables specified? Default to all    dii = dii(find(dii ~= 0));                           % Data interation index
177      if isempty(varlist), varlist=ncnames(var(nc)); end    if dBug > 0,
178            fprintf(' ; fii='); fprintf(' %i',fii);
179      % Attributes for structure      fprintf(' ; dii='); fprintf(' %i',dii); fprintf(' \n');
180      if iters>0; S.iters_read_from_file=iters; end    end
181      S.attributes.global=read_att(nc);      
182          % No variables specified? Default to all
183          % Read variable data    if isempty(varlist), varlist=ncnames(var(nc)); end
184          for ivar=1:size(varlist,2)    
185              % Attributes for structure
186          cvar=char(varlist{ivar});    if iters>0; S.iters_from_file=iters; end
187          if isempty(nc{cvar})    S.attributes.global=read_att(nc);
188              disp(['No such variable ''',cvar,''' in MNC file ',name(nc)]);    [pstr,netcdf_fname,ext] = fileparts(name(nc));
189              continue    if strcmp(netcdf_fname(end-3:end),'glob')
190          end      % assume it is a global file produced by gluemnc and change some
191                % attributes
192          dims = ncnames(dim(nc{cvar}));        % Dimensions      S.attributes.global.sNx = S.attributes.global.Nx;
193          adj = 0;      S.attributes.global.sNy = S.attributes.global.Ny;
194          if dims{1} == 'T'      S.attributes.global.nPx = 1;
195              if isempty(find(fii)), return, end      S.attributes.global.nSx = 1;
196                    S.attributes.global.nPy = 1;
197              tmpdata = nc{cvar}(fii,:);      S.attributes.global.nSy = 1;
198              if ismember('Zd000001' ,dims), adj = adj - 1; end      S.attributes.global.tile_number = 1;
199              if ismember('Zmd000001',dims), adj = adj - 1; end      S.attributes.global.nco_openmp_thread_number = 1;
200              it = length(dims) + adj;    end
201          else      
202              tmpdata = nc{cvar}(:);    % Read variable data
203              it = 0;    for ivar=1:size(varlist,2)
204        
205        cvar=char(varlist{ivar});
206        if isempty(nc{cvar})
207          disp(['No such variable ''',cvar,''' in MNC file ',name(nc)]);
208          continue
209        end
210        % code by Bruno Deremble: if you do not want to read all tiles these
211        % modifications make the output field smaller, let us see, if it is
212        % robust
213        if (isfield(S,cvar) == 0); firstiter = 1; else firstiter = 0; end
214        % end code by Bruno Deremble
215    
216        dims = ncnames(dim(nc{cvar}));        % Dimensions
217        sizVar = size(nc{cvar}); nDims=length(sizVar);
218        if dims{1} == 'T'
219          if isempty(find(fii)), error('Iters not found'); end
220          it = length(dims);
221          tmpdata = nc{cvar}(fii,:);
222          % leading unity dimensions get lost; add them back:
223          tmpdata=reshape(tmpdata,[length(fii) sizVar(2:end)]);
224        else
225          it = 0;
226          tmpdata = nc{cvar}(:);
227          % leading unity dimensions get lost; add them back:
228          tmpdata=reshape(tmpdata,sizVar);
229        end
230        
231        if dBug > 1,
232          fprintf(['  var:',cvar,': nDims=%i ('],nDims);fprintf(' %i',size(nc{cvar}));
233          fprintf('):%i,nD=%i,it=%i ;',length(size(tmpdata)),length(dims),it);
234        end
235        if length(dims) > 1,
236          tmpdata=permute(tmpdata,[nDims:-1:1]);
237        end
238        if dBug > 1,
239    %         fprintf('(tmpdata:');fprintf(' %i',size(tmpdata)); fprintf(')');
240        end
241        [ni nj nk nm nn no np]=size(tmpdata);
242          
243        [i0,j0,fn]=findTileOffset(S);
244        cdim=dims{end}; if cdim(1)~='X'; i0=0; end
245        cdim=dims{end}; if cdim(1)=='Y'; i0=j0; j0=0; end
246        if length(dims)>1;
247          cdim=dims{end-1}; if cdim(1)~='Y'; j0=0; end
248        else
249          j0=0;
250        end
251        if dBug > 1,
252          fprintf(' i0,ni= %i %i; j,nj= %i %i; nk=%i :',i0,ni,j0,nj,nk);
253        end
254        % code by Bruno Deremble: if you do not want to read all tiles these
255        % modifications make the output field smaller, let us see, if it is
256        % robust
257        if (firstiter)
258          S.i_first.(cvar) = i0;
259          S.j_first.(cvar) = j0;
260        end
261        i0 = i0 - S.i_first.(cvar);
262        j0 = j0 - S.j_first.(cvar);
263        % end code by Bruno Deremble
264        
265        Sstr = '';
266        for istr = 1:max(nDims,length(dims)),
267          if     istr == it,  Sstr = [Sstr,'dii,'];
268          elseif istr == 1,   Sstr = [Sstr,'i0+(1:ni),'];
269          elseif istr == 2,   Sstr = [Sstr,'j0+(1:nj),'];
270          elseif istr == 3,   Sstr = [Sstr,'(1:nk),'];
271          elseif istr == 4,   Sstr = [Sstr,'(1:nm),'];
272          elseif istr == 5,   Sstr = [Sstr,'(1:nn),'];
273          elseif istr == 6,   Sstr = [Sstr,'(1:no),'];
274          elseif istr == 7,   Sstr = [Sstr,'(1:np),'];
275          else, error('Can''t handle this many dimensions!');
276          end
277        end
278        eval(['S.(cvar)(',Sstr(1:end-1),')=tmpdata;'])
279        %S.(cvar)(i0+(1:ni),j0+(1:nj),(1:nk),(1:nm),(1:nn),(1:no),(1:np))=tmpdata;
280        if dBug > 1, fprintf(' %i',size(S.(cvar))); fprintf('\n'); end
281        
282        S.attributes.(cvar)=read_att(nc{cvar});
283        % replace missing or FillValues with NaN
284        if isfield(S.attributes.(cvar),'missing_value');
285          misval = S.attributes.(cvar).missing_value;
286          S.(cvar)(S.(cvar) == misval) = NaN;
287        end
288        if isfield(S.attributes.(cvar),'FillValue_');
289          misval = S.attributes.(cvar).FillValue_;
290          S.(cvar)(S.(cvar) == misval) = NaN;
291        end
292        
293      end % for ivar
294        
295      if isempty(S)
296        error('Something didn''t work!!!');
297      end
298      
299      return
300    
301    function [S] = rdmnc_local_matlabAPI(fname,varlist,iters,S,dBug)
302    
303      fiter = ncgetvar(fname,'iter');                     % File iterations present
304      if isempty(fiter), fiter = ncgetvar(fname,'T'); end
305      if isinf(iters); iters = fiter(end); end
306      if isnan(iters); iters = fiter; end
307      [fii,dii] = ismember(fiter,iters);  fii = find(fii); % File iteration index
308      dii = dii(find(dii ~= 0));                           % Data interation index
309      if dBug > 0,
310        fprintf(' ; fii='); fprintf(' %i',fii);
311        fprintf(' ; dii='); fprintf(' %i',dii); fprintf(' \n');
312      end
313    
314      % now open the file for reading
315      nc = netcdf.open(fname,'NC_NOWRITE');
316      % get basic information about netcdf file
317      [ndims nvars natts recdim] = netcdf.inq(nc);
318    
319      % No variables specified? Default to all
320      if isempty(varlist),
321        for k=0:nvars-1
322          varlist{k+1} = netcdf.inqVar(nc,k);
323        end
324      end
325      
326      % Attributes for structure
327      if iters>0; S.iters_from_file=iters; end
328      S.attributes.global=ncgetatt(nc,'global');
329      [pstr,netcdf_fname,ext] = fileparts(fname);
330      if strcmp(netcdf_fname(end-3:end),'glob')
331        % assume it is a global file produced by gluemnc and change some
332        % attributes
333        S.attributes.global.sNx = S.attributes.global.Nx;
334        S.attributes.global.sNy = S.attributes.global.Ny;
335        S.attributes.global.nPx = 1;
336        S.attributes.global.nSx = 1;
337        S.attributes.global.nPy = 1;
338        S.attributes.global.nSy = 1;
339        S.attributes.global.tile_number = 1;
340        S.attributes.global.nco_openmp_thread_number = 1;
341      end
342        
343      % Read variable data
344      for ivar=1:size(varlist,2)
345        
346        cvar=char(varlist{ivar});
347        varid=ncfindvarid(nc,cvar);
348        if isempty(varid)
349          disp(['No such variable ''',cvar,''' in MNC file ',fname]);
350          continue
351        end
352        % code by Bruno Deremble: if you do not want to read all tiles these
353        % modifications make the output field smaller, let us see, if it is
354        % robust
355        if (isfield(S,cvar) == 0); firstiter = 1; else firstiter = 0; end
356        % end code by Bruno Deremble
357        
358        [varname,xtype,dimids,natts] = netcdf.inqVar(nc,varid);
359        % does this variable contain a record (unlimited) dimension?
360        [isrecvar,recpos] = ismember(recdim,dimids);
361        
362        % Dimensions
363        clear sizVar dims
364        for k=1:length(dimids)
365          [dims{k}, sizVar(k)] = netcdf.inqDim(nc,dimids(k));
366        end
367        nDims=length(sizVar);
368        if isrecvar
369          if isempty(find(fii)), error('Iters not found'); end
370          it = length(dims);
371          if length(dimids) == 1
372            % just a time or iteration variable, this will result in a vector
373            % and requires special treatment
374            icount=1;
375            tmpdata = zeros(length(fii),1);
376            for k=1:length(fii)
377              istart = fii(k)-1;
378              tmpdata(k) = netcdf.getVar(nc,varid,istart,icount,'double');
379          end          end
380                  else
381          tmpdata=squeeze(permute(tmpdata,[9:-1:1]));          % from now on we assume that the record dimension is always the
382          [ni nj nk nm nn no np]=size(tmpdata);          % last dimension. This may not always be the case
383                    if recpos ~= nDims
384          [i0,j0,fn]=findTileOffset(S);            error(sprintf('%s\n%s\n%s%s%i%s%i', ...
385          cdim=dims{end}; if cdim(1)~='X'; i0=0; end                          ['The current code assumes that the record ' ...
386          cdim=dims{end}; if cdim(1)=='Y'; i0=j0; j0=0; end                           'dimension is the last dimension,'], ...
387          if length(dims)>1;                          'this is not the case for variable', cvar, ...
388              cdim=dims{end-1}; if cdim(1)~='Y'; j0=0; end                          ': nDims = ', nDims,  ...
389          else                          ', position of recDim = ', recpos))
             j0=0;  
390          end          end
391                    istart = zeros(1,it); % indexing starts a 0
392          Sstr = '';          icount = sizVar;
393          for istr = 1:7          % we always want to get only on time slice at a time
394              if     istr == it,  Sstr = [Sstr,'dii,'];          icount(recpos) = 1;
395              elseif istr == 1,   Sstr = [Sstr,'i0+(1:ni),'];          % make your life simpler by putting the time dimension first
396              elseif istr == 2,   Sstr = [Sstr,'j0+(1:nj),'];          tmpdata = zeros([length(fii) sizVar(1:end-1)]);
397              elseif istr == 3,   Sstr = [Sstr,'(1:nk),'];          for k=1:length(fii)
398              elseif istr == 4,   Sstr = [Sstr,'(1:nm),'];            istart(recpos) = fii(k)-1; % indexing starts at 0
399              elseif istr == 5,   Sstr = [Sstr,'(1:nn),'];            tmp = netcdf.getVar(nc,varid,istart,icount,'double');
400              elseif istr == 6,   Sstr = [Sstr,'(1:no),'];            tmpdata(k,:) = tmp(:);
401              elseif istr == 7,   Sstr = [Sstr,'(1:np),'];          end
402              else, error('Can''t handle this many dimensions!');          % move time dimension to the end ...
403              end          tmpdata = shiftdim(tmpdata,1);
404          end          % ... and restore original shape
405                    tmpdata = reshape(tmpdata,[sizVar(1:end-1) length(fii)]);
406          eval(['S.(cvar)(',Sstr(1:end-1),')=tmpdata;'])        end
407          %S.(cvar)(i0+(1:ni),j0+(1:nj),(1:nk),(1:nm),(1:nn),(1:no),(1:np))=tmpdata;      else
408          S.attributes.(cvar)=read_att(nc{cvar});        it = 0;
409          end        tmpdata = netcdf.getVar(nc,varid,'double');
410        end
411        %
412        if dBug > 1,
413          fprintf(['  var:',cvar,': nDims=%i ('],nDims);fprintf(' %i',sizVar);
414          fprintf('):%i,nD=%i,it=%i ;',length(size(tmpdata)),length(dims),it);
415        end
416        [ni nj nk nm nn no np]=size(tmpdata);
417        %
418        [i0,j0,fn]=findTileOffset(S);
419        cdim=dims{1}; if cdim(1)~='X'; i0=0; end
420        cdim=dims{1}; if cdim(1)=='Y'; i0=j0; j0=0; end
421        if length(dims)>1;
422          cdim=dims{2}; if cdim(1)~='Y'; j0=0; end
423        else
424          j0=0;
425        end
426        if dBug > 1,
427          fprintf(' i0,ni= %i %i; j,nj= %i %i; nk=%i :',i0,ni,j0,nj,nk);
428        end
429        % code by Bruno Deremble: if you do not want to read all tiles these
430        % modifications make the output field smaller, let us see, if it is
431        % robust
432        if (firstiter)
433          S.i_first.(cvar) = i0;
434          S.j_first.(cvar) = j0;
435        end
436        i0 = i0 - S.i_first.(cvar);
437        j0 = j0 - S.j_first.(cvar);
438        % end code by Bruno Deremble
439    
440        Sstr = '';
441        for istr = 1:max(nDims,length(dims)),
442          if     istr == it,  Sstr = [Sstr,'dii,'];
443          elseif istr == 1,   Sstr = [Sstr,'i0+(1:ni),'];
444          elseif istr == 2,   Sstr = [Sstr,'j0+(1:nj),'];
445          elseif istr == 3,   Sstr = [Sstr,'(1:nk),'];
446          elseif istr == 4,   Sstr = [Sstr,'(1:nm),'];
447          elseif istr == 5,   Sstr = [Sstr,'(1:nn),'];
448          elseif istr == 6,   Sstr = [Sstr,'(1:no),'];
449          elseif istr == 7,   Sstr = [Sstr,'(1:np),'];
450          else, error('Can''t handle this many dimensions!');
451          end
452        end
453        eval(['S.(cvar)(',Sstr(1:end-1),')=tmpdata;'])
454        %S.(cvar)(i0+(1:ni),j0+(1:nj),(1:nk),(1:nm),(1:nn),(1:no),(1:np))=tmpdata;
455        if dBug > 1, fprintf(' %i',size(S.(cvar))); fprintf('\n'); end
456        %
457        S.attributes.(cvar)=ncgetatt(nc,cvar);
458        % replace missing or FillValues with NaN
459        if isfield(S.attributes.(cvar),'missing_value');
460          misval = S.attributes.(cvar).missing_value;
461          S.(cvar)(S.(cvar) == misval) = NaN;
462        end
463        if isfield(S.attributes.(cvar),'FillValue_');
464          misval = S.attributes.(cvar).FillValue_;
465          S.(cvar)(S.(cvar) == misval) = NaN;
466        end
467        
468      end % for ivar
469    
470  if isempty(S)    % close the file
471      netcdf.close(nc);
472      
473      if isempty(S)
474      error('Something didn''t work!!!');      error('Something didn''t work!!!');
475  end    end
476      
477      return
478    
479    function vf = ncgetvar(fname,varname)
480    % read a netcdf variable
481      
482      nc=netcdf.open(fname,'NC_NOWRITE');
483      % find out basics about the files
484      [ndims nvars natts dimm] = netcdf.inq(nc);
485      vf = [];
486      varid = [];
487      for k=0:nvars-1
488        if strcmp(netcdf.inqVar(nc,k),varname)
489          varid = netcdf.inqVarID(nc,varname);
490        end
491      end
492      if ~isempty(varid);
493        [varn,xtype,dimids,natts] = netcdf.inqVar(nc,varid);
494        % get data
495        vf = double(netcdf.getVar(nc,varid));
496      else
497        % do nothing
498      end
499      netcdf.close(nc);
500      
501      return
502    
503    function misval = ncgetmisval(nc,varid)
504    
505      [varname,xtype,dimids,natts] = netcdf.inqVar(nc,varid);
506      misval = [];
507      for k=0:natts-1
508        attn = netcdf.inqAttName(nc,varid,k);
509        if strcmp(attn,'missing_value') | strcmp(attn,'_FillValue')
510          misval = double(netcdf.getAtt(nc,varid,attname));
511        end
512      end
513      
514    function A = ncgetatt(nc,varname)
515    % get all attributes and put them into a struct
516      
517    % 1. get global properties of file
518      [ndims nvars natts dimm] = netcdf.inq(nc);
519    
520      % get variable ID and properties
521      if strcmp('global',varname)
522        % "variable" is global
523        varid = netcdf.getConstant('NC_GLOBAL');
524      else
525        % find variable ID and properties
526        varid = ncfindvarid(nc,varname);
527        if ~isempty(varid)
528          [varn,xtype,dimids,natts] = netcdf.inqVar(nc,varid);
529        else
530          warning(sprintf('variable %s not found',varname))
531        end
532      end
533    
534      if natts > 1
535        for k=0:natts-1
536          attn = netcdf.inqAttName(nc,varid,k);
537          [xtype attlen]=netcdf.inqAtt(nc,varid,attn);
538          attval = netcdf.getAtt(nc,varid,attn);
539          if ~ischar(attval)
540            attval = double(attval);
541          end
542          if strcmp(attn,'_FillValue')
543            % matlab does not allow variable names to begin with an
544            % underscore ("_"), so we have to do change the name of this
545            % obsolete attribute.
546            A.FillValue_=attval;
547          else
548            A.(char(attn))=attval;
549          end
550        end
551      else
552          A = 'none';
553      end
554      
555      return
556    
557      
558    function varid = ncfindvarid(nc,varname)
559    
560      [ndims nvars natts dimm] = netcdf.inq(nc);
561      varid=[];
562      for k=0:nvars-1
563        if strcmp(netcdf.inqVar(nc,k),varname);
564          varid = netcdf.inqVarID(nc,varname);
565          break
566        end
567      end
568    
569      return
Legend:
Removed from v.1.7  
changed lines
  Added in v.1.25
  ViewVC Help
Powered by ViewVC 1.1.22