/[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.2 by mlosch, Wed Feb 2 10:31:22 2005 UTC revision 1.7 by enderton, Fri May 12 16:52:09 2006 UTC
# Line 1  Line 1 
1  function [S] = rdmnc(varargin)  function [S] = rdmnc(varargin)
2  % S=RDMNC(FILE1,FILE2,...)  
3  % S=RDMNC(FILE1,...,ITER)  % Usage:
4  % S=RDMNC(FILE1,...,'VAR1','VAR2',...)  %   S=RDMNC(FILE1,FILE2,...)
5  % S=RDMNC(FILE1,...,'VAR1','VAR2',...,ITER)  %   S=RDMNC(FILE1,...,ITER)
6    %   S=RDMNC(FILE1,...,'VAR1','VAR2',...)
7    %   S=RDMNC(FILE1,...,'VAR1','VAR2',...,ITER)
8    %
9    % Input:
10    %   FILE1   Either a single file name (e.g. 'state.nc') or a wild-card
11    %           strings expanding to a group of file names (e.g. 'state.*.nc').
12    %           There are no assumptions about suffices (e.g. 'state.*' works).
13    %   VAR1    Model variable names as written in the MNC/netcdf file.
14    %   ITER    Vector of iterations in the MNC/netcdf files, not model time.
15  %  %
16  %  FILE1 ... is either a single file name (e.g. 'state.nc') or a wild-card  % Output:
17  %        strings expanding to a group of file names (e.g. 'state.*.nc').  %   S       Structure with fields corresponding to 'VAR1', 'VAR2', ...
 %        There are no assumptions about suffices (e.g. 'state.*' works)  
 %  VAR1 ... are model variable names as written in the MNC/netcdf file  
 %  ITER  is a vector of time levels (consecutive indices referring to snap-shot  
 %        in the MNC/netcdf file i.e. 1,2,3,... and not model time)  
 %  S     is a structure with mutliple fields  
18  %  %
19  % rdmnc() is a rudimentary wrapper for joining and reading netcdf files  % Description:
20  % produced by MITgcm.  It does not give the same flexibility as opening the  %   This function is a rudimentary wrapper for joining and reading netcdf
21  % netcdf files directly using netcdf(). It is useful for quick loading of  %   files produced by MITgcm.  It does not give the same flexibility as
22  % entire model fields which are distributed in multiple netcdf files.  %   opening the netcdf files directly using netcdf(), but is useful for
23  % rdmnc() also reads non-MITgcm generated netcdf files.  %   quick loading of entire model fields which are distributed in multiple
24    %   netcdf files.
25  %  %
26  % e.g. To plot the surface temperature in last time-level written to file  % Example:
27  % >> S=rdmnd('state.*','XC','YC','RC','T',Inf);  %   >> S=rdmnd('state.*','XC','YC','T');
28  % >> imagesc( S.XC, S.YC, S.T(:,:,1)' );  %   >> imagesc( S.XC, S.YC, S.T(:,:,1)' );
29  %  %
30  % $Header$  %  Author:  Alistair Adcroft
31    %  Modifications:  Daniel Enderton
32    
33  %Defaults  % Initializations
 global verbose  
34  file={};  file={};
35    filepaths={};
36  files={};  files={};
 verbose=0;  
37  varlist={};  varlist={};
38  timelevels=[];  iters=[];
39    
40  % Process function arguments  % Process function arguments
41  for iarg=1:nargin;  for iarg=1:nargin;
42   arg=varargin{iarg};      arg=varargin{iarg};
43   if ischar(arg)      if ischar(arg)
44     switch char(arg)          if isempty(dir(char(arg)))
45       case 'verbose',              varlist{end+1}=arg;
46         verbose=1;          else
47       otherwise,              file{end+1}=arg;
48         if isempty( dir(char(arg)) )          end
49          varlist{end+1}={arg};      else
50         else          if isempty(iters)
51          file{end+1}={arg};              iters=arg;
52         end          else
53     end              error(['The only allowed numeric argument is iterations',...
54   else                     ' to read in as a vector for the last arguement']);
55     if isempty(timelevels)          end
56       timelevels=arg;      end
    else  
      error('Specify the time levels in a vector and not as multiple numeric arguments')  
    end  
  end  
57  end  end
58    
 if verbose; disp('Verbose mode enabled'); end  
   
59  % Create list of filenames  % Create list of filenames
60  for eachfile=file  for eachfile=file
61   expandedEachFile=dir(char(eachfile{1}));          filepathtemp=eachfile{:};
62   for i=1:size(expandedEachFile,1);          indecies = find(filepathtemp=='/');
63    if expandedEachFile(i).isdir==0; files{end+1}=expandedEachFile(i).name; end          if ~isempty(indecies)
64   end          filepathtemp = filepathtemp(1:indecies(end));
65            else
66            filepathtemp = '';
67            end
68        expandedEachFile=dir(char(eachfile{1}));
69        for i=1:size(expandedEachFile,1);
70            if expandedEachFile(i).isdir==0
71                files{end+1}=expandedEachFile(i).name;
72                filepaths{end+1}=filepathtemp;
73            end
74        end
75  end  end
76    
 % Open each file  
 S.attributes=[];  
 for eachfile={files{1:end}}  
  nc=netcdf(char(eachfile),'read');  
  if ismncfile(nc)  
   S=rdmnc_local(nc,varlist,timelevels,S);  
  else  
   S=rdnetcdf_local(nc,varlist,S);  
  end  
 end  
77    
78  % -----------------------------------------------------------------------------  % If iterations unspecified, open all the files and make list of all the
79  function [result] = ismncfile(nc);  % iterations that appear, use this iterations list for data extraction.
80  result=~isempty(nc.MITgcm_mnc_ver);  if isempty(iters)
81  %MLresult=~isempty(nc.('MITgcm_mnc_ver'));      iters = [];
82  % -----------------------------------------------------------------------------      for ieachfile=1:length(files)
83  function [A] = read_att(nc);          eachfile = [filepaths{ieachfile},files{ieachfile}];
84  global verbose          nc=netcdf(char(eachfile),'read');
85  allatt=ncnames(att(nc)); if verbose; allatt, end          nciters = nc{'iter'}(:);
86  A='none';          if isempty(nciters), nciters = nc{'T'}(:); end
87  for attr=allatt;          iters = [iters,nciters];
88   tmp = char(attr);          close(nc);
89   eval(['A.' tmp '= nc.' tmp '(:);'])      end
90  %ML A.(char(attr))=nc.(char(attr))(:);      iters = unique(iters);
 end  
 % -----------------------------------------------------------------------------  
 function [i0,j0,fn] = findTileOffset(S);  
 global verbose  
 fn=0;  
 if isfield(S,'attributes') & isfield(S.attributes,'global')  
  G=S.attributes.global;  
  tn=G.tile_number;  
  snx=G.sNx; sny=G.sNy; nsx=G.nSx; nsy=G.nSy; npx=G.nPx; npy=G.nPy;  
  ntx=nsx*npx;nty=nsy*npy;  
  gbi=mod(tn-1,ntx); gbj=(tn-gbi-1)/ntx;  
  i0=snx*gbi; j0=sny*gbj;  
  if isfield(S.attributes.global,'exch2_myFace')  
   fn=G.exch2_myFace;  
  end  
 else  
  i0=0;j0=0;  
 end  
 % -----------------------------------------------------------------------------  
 function [S] = rdnetcdf_local(nc,varlist,S)  
 % Read all attributes and variable data from a netcdf file  
 % with special operations for MNC style data  
 global verbose  
   
 % No variables specified? Default to all  
 if isempty(varlist)  
  varlist=ncnames(var(nc)); if verbose; varlist, end  
91  end  end
92    
93  % Attributes for structure  % Cycle through files for data extraction.
94  S.attributes=read_att(nc);  S.attributes=[];
95    for ieachfile=1:length(files)
96  % Read variable data      eachfile = [filepaths{ieachfile},files{ieachfile}];
97  for ivar=1:size(varlist,2);      nc=netcdf(char(eachfile),'read');
98   cvar=char(varlist{ivar});      S=rdmnc_local(nc,varlist,iters,S);
99   tmpdata=nc{cvar}(:);      close(nc);
  if isempty(tmpdata)  
   disp(['No such variable ''' cvar ''' in netcdf file' name(nc)])  
  else  
   tmpdata=squeeze(permute(tmpdata,[9:-1:1]));  
   eval(['S.' cvar '=tmpdata;'])  
   eval(['S.attributes.' cvar '=read_att(nc{' cvar '});'])  
 %ML  S.(cvar)=tmpdata;  
 %ML  S.attributes.(cvar)=read_att(nc{cvar});  
  end  
 end  
 % -----------------------------------------------------------------------------  
 function [S] = rdmnc_local(nc,varlist,timelevels,S)  
 % Read all attributes and variable data from a netcdf file  
 % with special operations for MNC style data  
 global verbose  
   
 nt=size( nc('T'), 1); if verbose; nt, end  
   
 % No variables specified? Default to all  
 if isempty(varlist)  
  varlist=ncnames(var(nc)); if verbose; varlist, end  
100  end  end
101    
 % No iterations specified? Default to all  
 if isempty(timelevels) | isnan(timelevels)  
  timelevels=1:nt;  
 elseif timelevels == Inf  
  timelevels=nt;  
 end  
102    
103  % Sanity checks  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
104  if max( find(timelevels > nt) )  %                             Local functions                             %
105   error('Requested time level is beyond time dimension in netcdf file')  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 end  
106    
107  % Attributes for structure  function [A] = read_att(nc);
108  if timelevels>0; S.timelevels_read_from_file=timelevels; end      allatt=ncnames(att(nc));
109  S.attributes.global=read_att(nc);      if ~isempty(allatt)
110            for attr=allatt;
111  % Read variable data              A.(char(attr))=nc.(char(attr))(:);
112  for ivar=1:size(varlist,2);          end
113   cvar=char(varlist{ivar}); if verbose; cvar, end      else
114   if isempty(nc{cvar})          A = 'none';
115    disp(['No such variable ''' cvar ''' in MNC file ' name(nc)])      end
116    continue  
117   end  function [i0,j0,fn] = findTileOffset(S);
118   dims=ncnames(dim(nc{cvar}));      fn=0;
119   if dims{1}=='T'      if isfield(S,'attributes') & isfield(S.attributes,'global')
120    if verbose; disp(['Reading variable ''' cvar ''' with record dimension']); end          G=S.attributes.global;
121    tmpdata=nc{cvar}(timelevels,:);          tn=G.tile_number;
122   else          snx=G.sNx; sny=G.sNy; nsx=G.nSx; nsy=G.nSy; npx=G.nPx; npy=G.nPy;
123    if verbose; disp(['Reading variable ''' cvar '''']); end          ntx=nsx*npx;nty=nsy*npy;
124    tmpdata=nc{cvar}(:);          gbi=mod(tn-1,ntx); gbj=(tn-gbi-1)/ntx;
125   end          i0=snx*gbi; j0=sny*gbj;
126   if isempty(tmpdata)          if isfield(S.attributes.global,'exch2_myFace')
127    error(['No such variable ''' cvar ''' in MNC file ' name(nc)])              fn=G.exch2_myFace;
128   else          end
129    tmpdata=squeeze(permute(tmpdata,[9:-1:1]));      else
130    [ni nj nk nm nn no np]=size(tmpdata);          i0=0;j0=0;
131    if np~=1; error('Wow! This is a very high dimension variable...'); end      end
132    [i0,j0,fn]=findTileOffset(S);      %[snx,sny,nsx,nsy,npx,npy,ntx,nty,i0,j0,fn];
133    cdim=dims{end}; if cdim(1)~='X'; i0=0; end  
134    cdim=dims{end}; if cdim(1)=='Y'; i0=j0; j0=0; end  function [S] = rdmnc_local(nc,varlist,iters,S)
135    if length(dims)>1;  
136     cdim=dims{end-1}; if cdim(1)~='Y'; j0=0; end      fiter = nc{'iter'}(:);                               % File iterations present
137    else      if isempty(fiter), fiter = nc{'T'}(:); end
138     j0=0;      [fii,dii] = ismember(fiter,iters);  fii = find(fii); % File iteration index
139    end      dii = dii(find(dii ~= 0));                           % Data interation index
140    eval(['S.' cvar ...      
141          '(i0+(1:ni),j0+(1:nj),(1:nk),(1:nm),(1:nn),(1:no),(1:np))=tmpdata;'])      % No variables specified? Default to all
142    eval(['S.attributes.' cvar ' =read_att(nc{''' cvar '''});'])      if isempty(varlist), varlist=ncnames(var(nc)); end
143  %ML  S.(cvar)(i0+(1:ni),j0+(1:nj),(1:nk),(1:nm),(1:nn),(1:no),(1:np))=tmpdata;      
144  %ML  S.attributes.(cvar)=read_att(nc{cvar});      % Attributes for structure
145   end      if iters>0; S.iters_read_from_file=iters; end
146  end      S.attributes.global=read_att(nc);
147        
148            % Read variable data
149            for ivar=1:size(varlist,2)
150            
151            cvar=char(varlist{ivar});
152            if isempty(nc{cvar})
153                disp(['No such variable ''',cvar,''' in MNC file ',name(nc)]);
154                continue
155            end
156            
157            dims = ncnames(dim(nc{cvar}));        % Dimensions
158            adj = 0;
159            if dims{1} == 'T'
160                if isempty(find(fii)), return, end
161                
162                tmpdata = nc{cvar}(fii,:);
163                if ismember('Zd000001' ,dims), adj = adj - 1; end
164                if ismember('Zmd000001',dims), adj = adj - 1; end
165                it = length(dims) + adj;
166            else
167                tmpdata = nc{cvar}(:);
168                it = 0;
169            end
170            
171            tmpdata=squeeze(permute(tmpdata,[9:-1:1]));
172            [ni nj nk nm nn no np]=size(tmpdata);
173            
174            [i0,j0,fn]=findTileOffset(S);
175            cdim=dims{end}; if cdim(1)~='X'; i0=0; end
176            cdim=dims{end}; if cdim(1)=='Y'; i0=j0; j0=0; end
177            if length(dims)>1;
178                cdim=dims{end-1}; if cdim(1)~='Y'; j0=0; end
179            else
180                j0=0;
181            end
182            
183            Sstr = '';
184            for istr = 1:7
185                if     istr == it,  Sstr = [Sstr,'dii,'];
186                elseif istr == 1,   Sstr = [Sstr,'i0+(1:ni),'];
187                elseif istr == 2,   Sstr = [Sstr,'j0+(1:nj),'];
188                elseif istr == 3,   Sstr = [Sstr,'(1:nk),'];
189                elseif istr == 4,   Sstr = [Sstr,'(1:nm),'];
190                elseif istr == 5,   Sstr = [Sstr,'(1:nn),'];
191                elseif istr == 6,   Sstr = [Sstr,'(1:no),'];
192                elseif istr == 7,   Sstr = [Sstr,'(1:np),'];
193                else, error('Can''t handle this many dimensions!');
194                end
195            end
196            
197            eval(['S.(cvar)(',Sstr(1:end-1),')=tmpdata;'])
198            %S.(cvar)(i0+(1:ni),j0+(1:nj),(1:nk),(1:nm),(1:nn),(1:no),(1:np))=tmpdata;
199            S.attributes.(cvar)=read_att(nc{cvar});
200            end
201    
202  if isempty(S)  if isempty(S)
203   error('Something didn''t work!!!')      error('Something didn''t work!!!');
204  end  end
 % -----------------------------------------------------------------------------  
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.7
  ViewVC Help
Powered by ViewVC 1.1.22