/[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.27 - (show annotations) (download)
Tue Feb 12 21:54:01 2013 UTC (11 years, 4 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint64y, checkpoint64x, checkpoint64z, checkpoint64q, checkpoint64p, checkpoint64s, checkpoint64r, checkpoint64u, checkpoint64t, checkpoint64w, checkpoint64v, checkpoint64i, checkpoint64h, checkpoint64k, checkpoint64j, checkpoint64m, checkpoint64l, checkpoint64o, checkpoint64n, checkpoint64e, checkpoint64d, checkpoint64g, checkpoint64f, checkpoint65, checkpoint66g, checkpoint66f, checkpoint66e, checkpoint66d, checkpoint66c, checkpoint66b, checkpoint66a, checkpoint66o, checkpoint66n, checkpoint66m, checkpoint66l, checkpoint66k, checkpoint66j, checkpoint66i, checkpoint66h, checkpoint65z, checkpoint65x, checkpoint65y, checkpoint65r, checkpoint65s, checkpoint65p, checkpoint65q, checkpoint65v, checkpoint65w, checkpoint65t, checkpoint65u, checkpoint65j, checkpoint65k, checkpoint65h, checkpoint65i, checkpoint65n, checkpoint65o, checkpoint65l, checkpoint65m, checkpoint65b, checkpoint65c, checkpoint65a, checkpoint65f, checkpoint65g, checkpoint65d, checkpoint65e, HEAD
Changes since 1.26: +47 -46 lines
fix tile-offset when output was produced using EXCH2

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

  ViewVC Help
Powered by ViewVC 1.1.22