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

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

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


Revision 1.22 - (show 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 function [S] = rdmnc(varargin)
2
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 %
16 % Output:
17 % S Structure with fields corresponding to 'VAR1', 'VAR2', ...
18 %
19 % 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 %
26 % Example:
27 % >> S=rdmnd('state.*','XC','YC','T');
28 % >> imagesc( S.XC, S.YC, S.T(:,:,1)' );
29 %
30 % Author: Alistair Adcroft
31 % Modifications: Daniel Enderton
32
33 % $Header: /u/gcmpack/MITgcm/utils/matlab/rdmnc.m,v 1.21 2010/05/16 08:56:43 mlosch Exp $
34 % $Name: $
35
36 % Initializations
37 dBug=0;
38 file={};
39 filepaths={};
40 files={};
41 varlist={};
42 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
49 for iarg=1:nargin;
50 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 ' to read in as a vector for the last argument']);
63 end
64 end
65 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
78 for eachfile=file
79 filepathtemp=eachfile{:};
80 indecies = find(filepathtemp==filesep);
81 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 end
94
95
96 % 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 if usemexcdf
103 nc=netcdf(char(eachfile),'read');
104 nciters = nc{'iter'}(:);
105 if isempty(nciters), nciters = nc{'T'}(:); end
106 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
116 iters = unique(iters');
117 end
118
119 % Cycle through files for data extraction.
120 S.attributes=[];
121 for ieachfile=1:length(files)
122 eachfile = [filepaths{ieachfile},files{ieachfile}];
123 if dBug > 0, fprintf([' open: ',eachfile]); end
124 if usemexcdf
125 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
135
136 return
137
138 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
139 % Local functions %
140 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
141
142 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
152 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 function [S] = rdmnc_local(nc,varlist,iters,S,dBug)
170
171 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 % fprintf('(tmpdata:');fprintf(' %i',size(tmpdata)); fprintf(')');
235 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 error('Something didn''t work!!!');
282 end
283
284 return
285
286 function [S] = rdmnc_local_matlabAPI(fname,varlist,iters,S,dBug)
287
288 fiter = ncgetvar(fname,'iter'); % File iterations present
289 if isempty(fiter), fiter = ncgetvar(fname,'T'); end
290 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 % now open the file for reading
300 nc = netcdf.open(fname,'NC_NOWRITE');
301 % get basic information about netcdf file
302 [ndims nvars natts recdim] = netcdf.inq(nc);
303
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 % does this variable contain a record (unlimited) dimension?
340 [isrecvar,recpos] = ismember(recdim,dimids);
341
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 if isrecvar
349 if isempty(find(fii)), error('Iters not found'); end
350 it = length(dims);
351 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 else
388 it = 0;
389 tmpdata = netcdf.getVar(nc,varid,'double');
390 end
391 %
392 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 %
398 [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 %
410 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 %
427 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
440 % close the file
441 netcdf.close(nc);
442
443 if isempty(S)
444 error('Something didn''t work!!!');
445 end
446
447 return
448
449 function vf = ncgetvar(fname,varname)
450 % read a netcdf variable
451
452 nc=netcdf.open(fname,'NC_NOWRITE');
453 % find out basics about the files
454 [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 netcdf.close(nc);
470
471 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