/[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.26 - (show annotations) (download)
Tue Jun 28 09:36:42 2011 UTC (12 years, 10 months ago) by mlosch
Branch: MAIN
CVS Tags: checkpoint64a, checkpoint64c, checkpoint64b, checkpoint63p, checkpoint63q, checkpoint63r, checkpoint63s, checkpoint63l, checkpoint63m, checkpoint63n, checkpoint63o, checkpoint63h, checkpoint63i, checkpoint63j, checkpoint63k, checkpoint63d, checkpoint63e, checkpoint63f, checkpoint63g, checkpoint63a, checkpoint63b, checkpoint63c, checkpoint64, checkpoint63
Changes since 1.25: +9 -9 lines
move attributes i_first and j_first into the "attritubes" attribute,
so that my scripts do not break

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.25 2011/06/23 08:40:09 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 % 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.attributes.i_first.(cvar) = i0;
259 S.attributes.j_first.(cvar) = j0;
260 end
261 i0 = i0 - S.attributes.i_first.(cvar);
262 j0 = j0 - S.attributes.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
380 else
381 % from now on we assume that the record dimension is always the
382 % last dimension. This may not always be the case
383 if recpos ~= nDims
384 error(sprintf('%s\n%s\n%s%s%i%s%i', ...
385 ['The current code assumes that the record ' ...
386 'dimension is the last dimension,'], ...
387 'this is not the case for variable', cvar, ...
388 ': nDims = ', nDims, ...
389 ', position of recDim = ', recpos))
390 end
391 istart = zeros(1,it); % indexing starts a 0
392 icount = sizVar;
393 % we always want to get only on time slice at a time
394 icount(recpos) = 1;
395 % make your life simpler by putting the time dimension first
396 tmpdata = zeros([length(fii) sizVar(1:end-1)]);
397 for k=1:length(fii)
398 istart(recpos) = fii(k)-1; % indexing starts at 0
399 tmp = netcdf.getVar(nc,varid,istart,icount,'double');
400 tmpdata(k,:) = tmp(:);
401 end
402 % move time dimension to the end ...
403 tmpdata = shiftdim(tmpdata,1);
404 % ... and restore original shape
405 tmpdata = reshape(tmpdata,[sizVar(1:end-1) length(fii)]);
406 end
407 else
408 it = 0;
409 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.attributes.i_first.(cvar) = i0;
434 S.attributes.j_first.(cvar) = j0;
435 end
436 i0 = i0 - S.attributes.i_first.(cvar);
437 j0 = j0 - S.attributes.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 % close the file
471 netcdf.close(nc);
472
473 if isempty(S)
474 error('Something didn''t work!!!');
475 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

  ViewVC Help
Powered by ViewVC 1.1.22