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

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

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


Revision 1.22 - (hide annotations) (download)
Fri May 28 07:16:43 2010 UTC (14 years ago) by mlosch
Branch: MAIN
CVS Tags: checkpoint62k, checkpoint62j, checkpoint62i, checkpoint62h, checkpoint62o, checkpoint62n, checkpoint62m, checkpoint62l, checkpoint62r, checkpoint62q, checkpoint62p
Changes since 1.21: +45 -14 lines
fix a problem with new matlab/netcdf version. Now really read only the
"iters" that are specified in the parameter list

1 enderton 1.6 function [S] = rdmnc(varargin)
2 enderton 1.3
3     % Usage:
4     % S=RDMNC(FILE1,FILE2,...)
5     % 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 adcroft 1.1 %
16 enderton 1.3 % Output:
17     % S Structure with fields corresponding to 'VAR1', 'VAR2', ...
18 adcroft 1.1 %
19 enderton 1.3 % Description:
20     % This function is a rudimentary wrapper for joining and reading netcdf
21     % files produced by MITgcm. It does not give the same flexibility as
22     % opening the netcdf files directly using netcdf(), but is useful for
23     % quick loading of entire model fields which are distributed in multiple
24     % netcdf files.
25 adcroft 1.1 %
26 enderton 1.3 % Example:
27     % >> S=rdmnd('state.*','XC','YC','T');
28     % >> imagesc( S.XC, S.YC, S.T(:,:,1)' );
29 adcroft 1.1 %
30 enderton 1.4 % Author: Alistair Adcroft
31 enderton 1.3 % Modifications: Daniel Enderton
32 adcroft 1.1
33 mlosch 1.22 % $Header: /u/gcmpack/MITgcm/utils/matlab/rdmnc.m,v 1.21 2010/05/16 08:56:43 mlosch Exp $
34 jmc 1.12 % $Name: $
35    
36 enderton 1.3 % Initializations
37 jmc 1.13 dBug=0;
38 adcroft 1.1 file={};
39 enderton 1.3 filepaths={};
40 adcroft 1.1 files={};
41     varlist={};
42 enderton 1.3 iters=[];
43 adcroft 1.1
44 mlosch 1.20 % 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 adcroft 1.1 % Process function arguments
49     for iarg=1:nargin;
50 enderton 1.3 arg=varargin{iarg};
51     if ischar(arg)
52     if isempty(dir(char(arg)))
53     varlist{end+1}=arg;
54     else
55     file{end+1}=arg;
56     end
57     else
58     if isempty(iters)
59     iters=arg;
60     else
61     error(['The only allowed numeric argument is iterations',...
62 jmc 1.14 ' to read in as a vector for the last argument']);
63 enderton 1.3 end
64     end
65 adcroft 1.1 end
66 jmc 1.14 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 adcroft 1.1
77     % Create list of filenames
78     for eachfile=file
79 enderton 1.3 filepathtemp=eachfile{:};
80 mlosch 1.17 indecies = find(filepathtemp==filesep);
81 enderton 1.3 if ~isempty(indecies)
82     filepathtemp = filepathtemp(1:indecies(end));
83     else
84     filepathtemp = '';
85     end
86     expandedEachFile=dir(char(eachfile{1}));
87     for i=1:size(expandedEachFile,1);
88     if expandedEachFile(i).isdir==0
89     files{end+1}=expandedEachFile(i).name;
90     filepaths{end+1}=filepathtemp;
91     end
92     end
93 adcroft 1.1 end
94    
95    
96 enderton 1.3 % If iterations unspecified, open all the files and make list of all the
97     % iterations that appear, use this iterations list for data extraction.
98     if isempty(iters)
99     iters = [];
100     for ieachfile=1:length(files)
101     eachfile = [filepaths{ieachfile},files{ieachfile}];
102 mlosch 1.20 if usemexcdf
103     nc=netcdf(char(eachfile),'read');
104     nciters = nc{'iter'}(:);
105     if isempty(nciters), nciters = nc{'T'}(:); end
106 mlosch 1.21 close(nc);
107 mlosch 1.20 else
108 mlosch 1.21 % 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 mlosch 1.20 end
114 jmc 1.15 iters = [iters,nciters'];
115 enderton 1.3 end
116 jmc 1.15 iters = unique(iters');
117 adcroft 1.1 end
118    
119 enderton 1.3 % Cycle through files for data extraction.
120     S.attributes=[];
121     for ieachfile=1:length(files)
122     eachfile = [filepaths{ieachfile},files{ieachfile}];
123 jmc 1.13 if dBug > 0, fprintf([' open: ',eachfile]); end
124 mlosch 1.20 if usemexcdf
125     nc=netcdf(char(eachfile),'read');
126     S=rdmnc_local(nc,varlist,iters,S,dBug);
127     close(nc);
128     else
129 mlosch 1.21 % 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 mlosch 1.20 end
134 adcroft 1.1 end
135    
136 mlosch 1.20 return
137 enderton 1.5
138 enderton 1.3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
139     % Local functions %
140     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
141 adcroft 1.1
142 enderton 1.3 function [A] = read_att(nc);
143     allatt=ncnames(att(nc));
144     if ~isempty(allatt)
145     for attr=allatt;
146     A.(char(attr))=nc.(char(attr))(:);
147     end
148     else
149     A = 'none';
150     end
151 adcroft 1.1
152 enderton 1.3 function [i0,j0,fn] = findTileOffset(S);
153     fn=0;
154     if isfield(S,'attributes') & isfield(S.attributes,'global')
155     G=S.attributes.global;
156     tn=G.tile_number;
157     snx=G.sNx; sny=G.sNy; nsx=G.nSx; nsy=G.nSy; npx=G.nPx; npy=G.nPy;
158     ntx=nsx*npx;nty=nsy*npy;
159     gbi=mod(tn-1,ntx); gbj=(tn-gbi-1)/ntx;
160     i0=snx*gbi; j0=sny*gbj;
161     if isfield(S.attributes.global,'exch2_myFace')
162     fn=G.exch2_myFace;
163     end
164     else
165     i0=0;j0=0;
166     end
167     %[snx,sny,nsx,nsy,npx,npy,ntx,nty,i0,j0,fn];
168    
169 jmc 1.13 function [S] = rdmnc_local(nc,varlist,iters,S,dBug)
170 enderton 1.3
171 mlosch 1.19 fiter = nc{'iter'}(:); % File iterations present
172     if isempty(fiter), fiter = nc{'T'}(:); end
173     if isinf(iters); iters = fiter(end); end
174     if isnan(iters); iters = fiter; end
175     [fii,dii] = ismember(fiter,iters); fii = find(fii); % File iteration index
176     dii = dii(find(dii ~= 0)); % Data interation index
177     if dBug > 0,
178     fprintf(' ; fii='); fprintf(' %i',fii);
179     fprintf(' ; dii='); fprintf(' %i',dii); fprintf(' \n');
180     end
181    
182     % No variables specified? Default to all
183     if isempty(varlist), varlist=ncnames(var(nc)); end
184    
185     % Attributes for structure
186     if iters>0; S.iters_from_file=iters; end
187     S.attributes.global=read_att(nc);
188     [pstr,netcdf_fname,ext] = fileparts(name(nc));
189     if strcmp(netcdf_fname(end-3:end),'glob')
190     % assume it is a global file produced by gluemnc and change some
191     % attributes
192     S.attributes.global.sNx = S.attributes.global.Nx;
193     S.attributes.global.sNy = S.attributes.global.Ny;
194     S.attributes.global.nPx = 1;
195     S.attributes.global.nSx = 1;
196     S.attributes.global.nPy = 1;
197     S.attributes.global.nSy = 1;
198     S.attributes.global.tile_number = 1;
199     S.attributes.global.nco_openmp_thread_number = 1;
200     end
201    
202     % Read variable data
203     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    
211     dims = ncnames(dim(nc{cvar})); % Dimensions
212     sizVar = size(nc{cvar}); nDims=length(sizVar);
213     if dims{1} == 'T'
214     if isempty(find(fii)), error('Iters not found'); end
215     it = length(dims);
216     tmpdata = nc{cvar}(fii,:);
217     % leading unity dimensions get lost; add them back:
218     tmpdata=reshape(tmpdata,[length(fii) sizVar(2:end)]);
219     else
220     it = 0;
221     tmpdata = nc{cvar}(:);
222     % leading unity dimensions get lost; add them back:
223     tmpdata=reshape(tmpdata,sizVar);
224     end
225    
226     if dBug > 1,
227     fprintf([' var:',cvar,': nDims=%i ('],nDims);fprintf(' %i',size(nc{cvar}));
228     fprintf('):%i,nD=%i,it=%i ;',length(size(tmpdata)),length(dims),it);
229     end
230     if length(dims) > 1,
231     tmpdata=permute(tmpdata,[nDims:-1:1]);
232     end
233     if dBug > 1,
234 jmc 1.13 % fprintf('(tmpdata:');fprintf(' %i',size(tmpdata)); fprintf(')');
235 mlosch 1.19 end
236     [ni nj nk nm nn no np]=size(tmpdata);
237    
238     [i0,j0,fn]=findTileOffset(S);
239     cdim=dims{end}; if cdim(1)~='X'; i0=0; end
240     cdim=dims{end}; if cdim(1)=='Y'; i0=j0; j0=0; end
241     if length(dims)>1;
242     cdim=dims{end-1}; if cdim(1)~='Y'; j0=0; end
243     else
244     j0=0;
245     end
246     if dBug > 1,
247     fprintf(' i0,ni= %i %i; j,nj= %i %i; nk=%i :',i0,ni,j0,nj,nk);
248     end
249    
250     Sstr = '';
251     for istr = 1:max(nDims,length(dims)),
252     if istr == it, Sstr = [Sstr,'dii,'];
253     elseif istr == 1, Sstr = [Sstr,'i0+(1:ni),'];
254     elseif istr == 2, Sstr = [Sstr,'j0+(1:nj),'];
255     elseif istr == 3, Sstr = [Sstr,'(1:nk),'];
256     elseif istr == 4, Sstr = [Sstr,'(1:nm),'];
257     elseif istr == 5, Sstr = [Sstr,'(1:nn),'];
258     elseif istr == 6, Sstr = [Sstr,'(1:no),'];
259     elseif istr == 7, Sstr = [Sstr,'(1:np),'];
260     else, error('Can''t handle this many dimensions!');
261     end
262     end
263     eval(['S.(cvar)(',Sstr(1:end-1),')=tmpdata;'])
264     %S.(cvar)(i0+(1:ni),j0+(1:nj),(1:nk),(1:nm),(1:nn),(1:no),(1:np))=tmpdata;
265     if dBug > 1, fprintf(' %i',size(S.(cvar))); fprintf('\n'); end
266    
267     S.attributes.(cvar)=read_att(nc{cvar});
268     % replace missing or FillValues with NaN
269     if isfield(S.attributes.(cvar),'missing_value');
270     misval = S.attributes.(cvar).missing_value;
271     S.(cvar)(S.(cvar) == misval) = NaN;
272     end
273     if isfield(S.attributes.(cvar),'FillValue_');
274     misval = S.attributes.(cvar).FillValue_;
275     S.(cvar)(S.(cvar) == misval) = NaN;
276     end
277    
278     end % for ivar
279    
280     if isempty(S)
281 enderton 1.3 error('Something didn''t work!!!');
282 mlosch 1.19 end
283    
284     return
285 mlosch 1.20
286 mlosch 1.21 function [S] = rdmnc_local_matlabAPI(fname,varlist,iters,S,dBug)
287 mlosch 1.20
288 mlosch 1.21 fiter = ncgetvar(fname,'iter'); % File iterations present
289     if isempty(fiter), fiter = ncgetvar(fname,'T'); end
290 mlosch 1.20 if isinf(iters); iters = fiter(end); end
291     if isnan(iters); iters = fiter; end
292     [fii,dii] = ismember(fiter,iters); fii = find(fii); % File iteration index
293     dii = dii(find(dii ~= 0)); % Data interation index
294     if dBug > 0,
295     fprintf(' ; fii='); fprintf(' %i',fii);
296     fprintf(' ; dii='); fprintf(' %i',dii); fprintf(' \n');
297     end
298    
299 mlosch 1.21 % now open the file for reading
300     nc = netcdf.open(fname,'NC_NOWRITE');
301 mlosch 1.20 % get basic information about netcdf file
302 mlosch 1.22 [ndims nvars natts recdim] = netcdf.inq(nc);
303 mlosch 1.20
304     % No variables specified? Default to all
305     if isempty(varlist),
306     for k=0:nvars-1
307     varlist{k+1} = netcdf.inqVar(nc,k);
308     end
309     end
310    
311     % Attributes for structure
312     if iters>0; S.iters_from_file=iters; end
313     S.attributes.global=ncgetatt(nc,'global');
314     [pstr,netcdf_fname,ext] = fileparts(fname);
315     if strcmp(netcdf_fname(end-3:end),'glob')
316     % assume it is a global file produced by gluemnc and change some
317     % attributes
318     S.attributes.global.sNx = S.attributes.global.Nx;
319     S.attributes.global.sNy = S.attributes.global.Ny;
320     S.attributes.global.nPx = 1;
321     S.attributes.global.nSx = 1;
322     S.attributes.global.nPy = 1;
323     S.attributes.global.nSy = 1;
324     S.attributes.global.tile_number = 1;
325     S.attributes.global.nco_openmp_thread_number = 1;
326     end
327    
328     % Read variable data
329     for ivar=1:size(varlist,2)
330    
331     cvar=char(varlist{ivar});
332     varid=ncfindvarid(nc,cvar);
333     if isempty(varid)
334     disp(['No such variable ''',cvar,''' in MNC file ',fname]);
335     continue
336     end
337    
338     [varname,xtype,dimids,natts] = netcdf.inqVar(nc,varid);
339 mlosch 1.22 % does this variable contain a record (unlimited) dimension?
340     [isrecvar,recpos] = ismember(recdim,dimids);
341 mlosch 1.20
342     % Dimensions
343     clear sizVar dims
344     for k=1:length(dimids)
345     [dims{k}, sizVar(k)] = netcdf.inqDim(nc,dimids(k));
346     end
347     nDims=length(sizVar);
348 mlosch 1.22 if isrecvar
349 mlosch 1.20 if isempty(find(fii)), error('Iters not found'); end
350     it = length(dims);
351 mlosch 1.22 if length(dimids) == 1
352     % just a time or iteration variable, this will result in a vector
353     % and requires special treatment
354     icount=1;
355     tmpdata = zeros(length(fii),1);
356     for k=1:length(fii)
357     istart = fii(k)-1;
358     tmpdata(k) = netcdf.getVar(nc,varid,istart,icount,'double');
359     end
360     else
361     % from now on we assume that the record dimension is always the
362     % last dimension. This may not always be the case
363     if recpos ~= nDims
364     error(sprintf('%s\n%s\n%s%s%i%s%i', ...
365     ['The current code assumes that the record ' ...
366     'dimension is the last dimension,'], ...
367     'this is not the case for variable', cvar, ...
368     ': nDims = ', nDims, ...
369     ', position of recDim = ', recpos))
370     end
371     istart = zeros(1,it); % indexing starts a 0
372     icount = sizVar;
373     % we always want to get only on time slice at a time
374     icount(recpos) = 1;
375     % make your life simpler by putting the time dimension first
376     tmpdata = zeros([length(fii) sizVar(1:end-1)]);
377     for k=1:length(fii)
378     istart(recpos) = fii(k)-1; % indexing starts at 0
379     tmp = netcdf.getVar(nc,varid,istart,icount,'double');
380     tmpdata(k,:) = tmp(:);
381     end
382     % move time dimension to the end ...
383     tmpdata = shiftdim(tmpdata,1);
384     % ... and restore original shape
385     tmpdata = reshape(tmpdata,[sizVar(1:end-1) length(fii)]);
386     end
387 mlosch 1.20 else
388     it = 0;
389     tmpdata = netcdf.getVar(nc,varid,'double');
390     end
391 mlosch 1.22 %
392 mlosch 1.20 if dBug > 1,
393     fprintf([' var:',cvar,': nDims=%i ('],nDims);fprintf(' %i',sizVar);
394     fprintf('):%i,nD=%i,it=%i ;',length(size(tmpdata)),length(dims),it);
395     end
396     [ni nj nk nm nn no np]=size(tmpdata);
397 mlosch 1.22 %
398 mlosch 1.20 [i0,j0,fn]=findTileOffset(S);
399     cdim=dims{1}; if cdim(1)~='X'; i0=0; end
400     cdim=dims{1}; if cdim(1)=='Y'; i0=j0; j0=0; end
401     if length(dims)>1;
402     cdim=dims{2}; if cdim(1)~='Y'; j0=0; end
403     else
404     j0=0;
405     end
406     if dBug > 1,
407     fprintf(' i0,ni= %i %i; j,nj= %i %i; nk=%i :',i0,ni,j0,nj,nk);
408     end
409 mlosch 1.22 %
410 mlosch 1.20 Sstr = '';
411     for istr = 1:max(nDims,length(dims)),
412     if istr == it, Sstr = [Sstr,'dii,'];
413     elseif istr == 1, Sstr = [Sstr,'i0+(1:ni),'];
414     elseif istr == 2, Sstr = [Sstr,'j0+(1:nj),'];
415     elseif istr == 3, Sstr = [Sstr,'(1:nk),'];
416     elseif istr == 4, Sstr = [Sstr,'(1:nm),'];
417     elseif istr == 5, Sstr = [Sstr,'(1:nn),'];
418     elseif istr == 6, Sstr = [Sstr,'(1:no),'];
419     elseif istr == 7, Sstr = [Sstr,'(1:np),'];
420     else, error('Can''t handle this many dimensions!');
421     end
422     end
423     eval(['S.(cvar)(',Sstr(1:end-1),')=tmpdata;'])
424     %S.(cvar)(i0+(1:ni),j0+(1:nj),(1:nk),(1:nm),(1:nn),(1:no),(1:np))=tmpdata;
425     if dBug > 1, fprintf(' %i',size(S.(cvar))); fprintf('\n'); end
426 mlosch 1.22 %
427 mlosch 1.20 S.attributes.(cvar)=ncgetatt(nc,cvar);
428     % replace missing or FillValues with NaN
429     if isfield(S.attributes.(cvar),'missing_value');
430     misval = S.attributes.(cvar).missing_value;
431     S.(cvar)(S.(cvar) == misval) = NaN;
432     end
433     if isfield(S.attributes.(cvar),'FillValue_');
434     misval = S.attributes.(cvar).FillValue_;
435     S.(cvar)(S.(cvar) == misval) = NaN;
436     end
437    
438     end % for ivar
439 mlosch 1.21
440     % close the file
441     netcdf.close(nc);
442    
443 mlosch 1.20 if isempty(S)
444     error('Something didn''t work!!!');
445     end
446    
447     return
448    
449 mlosch 1.21 function vf = ncgetvar(fname,varname)
450 mlosch 1.20 % read a netcdf variable
451    
452 mlosch 1.21 nc=netcdf.open(fname,'NC_NOWRITE');
453     % find out basics about the files
454 mlosch 1.20 [ndims nvars natts dimm] = netcdf.inq(nc);
455     vf = [];
456     varid = [];
457     for k=0:nvars-1
458     if strcmp(netcdf.inqVar(nc,k),varname)
459     varid = netcdf.inqVarId(nc,varname);
460     end
461     end
462     if ~isempty(varid);
463     [varn,xtype,dimids,natts] = netcdf.inqVar(nc,varid);
464     % get data
465     vf = double(netcdf.getVar(nc,varid));
466     else
467     % do nothing
468     end
469 mlosch 1.21 netcdf.close(nc);
470    
471 mlosch 1.20 return
472    
473     function misval = ncgetmisval(nc,varid)
474    
475     [varname,xtype,dimids,natts] = netcdf.inqVar(nc,varid);
476     misval = [];
477     for k=0:natts-1
478     attn = netcdf.inqAttName(nc,varid,k);
479     if strcmp(attn,'missing_value') | strcmp(attn,'_FillValue')
480     misval = double(netcdf.getAtt(nc,varid,attname));
481     end
482     end
483    
484     function A = ncgetatt(nc,varname)
485     % get all attributes and put them into a struct
486    
487     % 1. get global properties of file
488     [ndims nvars natts dimm] = netcdf.inq(nc);
489    
490     % get variable ID and properties
491     if strcmp('global',varname)
492     % "variable" is global
493     varid = netcdf.getConstant('NC_GLOBAL');
494     else
495     % find variable ID and properties
496     varid = ncfindvarid(nc,varname);
497     if ~isempty(varid)
498     [varn,xtype,dimids,natts] = netcdf.inqVar(nc,varid);
499     else
500     warning(sprintf('variable %s not found',varname))
501     end
502     end
503    
504     if natts > 1
505     for k=0:natts-1
506     attn = netcdf.inqAttName(nc,varid,k);
507     [xtype attlen]=netcdf.inqAtt(nc,varid,attn);
508     attval = netcdf.getAtt(nc,varid,attn);
509     if ~ischar(attval)
510     attval = double(attval);
511     end
512     A.(char(attn))=attval;
513     end
514     else
515     A = 'none';
516     end
517    
518     return
519    
520    
521     function varid = ncfindvarid(nc,varname)
522    
523     [ndims nvars natts dimm] = netcdf.inq(nc);
524     varid=[];
525     for k=0:nvars-1
526     if strcmp(netcdf.inqVar(nc,k),varname);
527     varid = netcdf.inqVarId(nc,varname);
528     break
529     end
530     end
531    
532     return

  ViewVC Help
Powered by ViewVC 1.1.22