/[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.20 - (show annotations) (download)
Fri May 14 11:29:16 2010 UTC (14 years, 1 month ago) by mlosch
Branch: MAIN
Changes since 1.19: +239 -8 lines
adapt to be able to use the generic matlab netcdf API

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.19 2009/11/24 13:09:16 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 else
107 nc=netcdf.open(char(eachfile),'NC_NOWRITE');
108 nciters = ncgetvar(nc,'iter');
109 if isempty(nciters), nciters = ncgetvar(nc,'T'); end
110 end
111 iters = [iters,nciters'];
112 if usemexcdf
113 close(nc);
114 else
115 netcdf.close(nc);
116 end
117 end
118 iters = unique(iters');
119 end
120
121 % Cycle through files for data extraction.
122 S.attributes=[];
123 for ieachfile=1:length(files)
124 eachfile = [filepaths{ieachfile},files{ieachfile}];
125 if dBug > 0, fprintf([' open: ',eachfile]); end
126 if usemexcdf
127 nc=netcdf(char(eachfile),'read');
128 S=rdmnc_local(nc,varlist,iters,S,dBug);
129 close(nc);
130 else
131 nc=netcdf.open(char(eachfile),'NC_NOWRITE');
132 S=rdmnc_local_matlabAPI(char(eachfile),nc,varlist,iters,S,dBug);
133 netcdf.close(nc);
134 end
135 end
136
137 return
138
139 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
140 % Local functions %
141 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
142
143 function [A] = read_att(nc);
144 allatt=ncnames(att(nc));
145 if ~isempty(allatt)
146 for attr=allatt;
147 A.(char(attr))=nc.(char(attr))(:);
148 end
149 else
150 A = 'none';
151 end
152
153 function [i0,j0,fn] = findTileOffset(S);
154 fn=0;
155 if isfield(S,'attributes') & isfield(S.attributes,'global')
156 G=S.attributes.global;
157 tn=G.tile_number;
158 snx=G.sNx; sny=G.sNy; nsx=G.nSx; nsy=G.nSy; npx=G.nPx; npy=G.nPy;
159 ntx=nsx*npx;nty=nsy*npy;
160 gbi=mod(tn-1,ntx); gbj=(tn-gbi-1)/ntx;
161 i0=snx*gbi; j0=sny*gbj;
162 if isfield(S.attributes.global,'exch2_myFace')
163 fn=G.exch2_myFace;
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
212 dims = ncnames(dim(nc{cvar})); % Dimensions
213 sizVar = size(nc{cvar}); nDims=length(sizVar);
214 if dims{1} == 'T'
215 if isempty(find(fii)), error('Iters not found'); end
216 it = length(dims);
217 tmpdata = nc{cvar}(fii,:);
218 % leading unity dimensions get lost; add them back:
219 tmpdata=reshape(tmpdata,[length(fii) sizVar(2:end)]);
220 else
221 it = 0;
222 tmpdata = nc{cvar}(:);
223 % leading unity dimensions get lost; add them back:
224 tmpdata=reshape(tmpdata,sizVar);
225 end
226
227 if dBug > 1,
228 fprintf([' var:',cvar,': nDims=%i ('],nDims);fprintf(' %i',size(nc{cvar}));
229 fprintf('):%i,nD=%i,it=%i ;',length(size(tmpdata)),length(dims),it);
230 end
231 if length(dims) > 1,
232 tmpdata=permute(tmpdata,[nDims:-1:1]);
233 end
234 if dBug > 1,
235 % fprintf('(tmpdata:');fprintf(' %i',size(tmpdata)); fprintf(')');
236 end
237 [ni nj nk nm nn no np]=size(tmpdata);
238
239 [i0,j0,fn]=findTileOffset(S);
240 cdim=dims{end}; if cdim(1)~='X'; i0=0; end
241 cdim=dims{end}; if cdim(1)=='Y'; i0=j0; j0=0; end
242 if length(dims)>1;
243 cdim=dims{end-1}; if cdim(1)~='Y'; j0=0; end
244 else
245 j0=0;
246 end
247 if dBug > 1,
248 fprintf(' i0,ni= %i %i; j,nj= %i %i; nk=%i :',i0,ni,j0,nj,nk);
249 end
250
251 Sstr = '';
252 for istr = 1:max(nDims,length(dims)),
253 if istr == it, Sstr = [Sstr,'dii,'];
254 elseif istr == 1, Sstr = [Sstr,'i0+(1:ni),'];
255 elseif istr == 2, Sstr = [Sstr,'j0+(1:nj),'];
256 elseif istr == 3, Sstr = [Sstr,'(1:nk),'];
257 elseif istr == 4, Sstr = [Sstr,'(1:nm),'];
258 elseif istr == 5, Sstr = [Sstr,'(1:nn),'];
259 elseif istr == 6, Sstr = [Sstr,'(1:no),'];
260 elseif istr == 7, Sstr = [Sstr,'(1:np),'];
261 else, error('Can''t handle this many dimensions!');
262 end
263 end
264 eval(['S.(cvar)(',Sstr(1:end-1),')=tmpdata;'])
265 %S.(cvar)(i0+(1:ni),j0+(1:nj),(1:nk),(1:nm),(1:nn),(1:no),(1:np))=tmpdata;
266 if dBug > 1, fprintf(' %i',size(S.(cvar))); fprintf('\n'); end
267
268 S.attributes.(cvar)=read_att(nc{cvar});
269 % replace missing or FillValues with NaN
270 if isfield(S.attributes.(cvar),'missing_value');
271 misval = S.attributes.(cvar).missing_value;
272 S.(cvar)(S.(cvar) == misval) = NaN;
273 end
274 if isfield(S.attributes.(cvar),'FillValue_');
275 misval = S.attributes.(cvar).FillValue_;
276 S.(cvar)(S.(cvar) == misval) = NaN;
277 end
278
279 end % for ivar
280
281 if isempty(S)
282 error('Something didn''t work!!!');
283 end
284
285 return
286
287 function [S] = rdmnc_local_matlabAPI(fname,nc,varlist,iters,S,dBug)
288
289 fiter = ncgetvar(nc,'iter'); % File iterations present
290 if isempty(fiter), fiter = ncgetvar(nc,'T'); end
291 if isinf(iters); iters = fiter(end); end
292 if isnan(iters); iters = fiter; end
293 [fii,dii] = ismember(fiter,iters); fii = find(fii); % File iteration index
294 dii = dii(find(dii ~= 0)); % Data interation index
295 if dBug > 0,
296 fprintf(' ; fii='); fprintf(' %i',fii);
297 fprintf(' ; dii='); fprintf(' %i',dii); fprintf(' \n');
298 end
299
300 % get basic information about netcdf file
301 [ndims nvars natts dimm] = netcdf.inq(nc);
302
303 % No variables specified? Default to all
304 if isempty(varlist),
305 for k=0:nvars-1
306 varlist{k+1} = netcdf.inqVar(nc,k);
307 end
308 end
309
310 % Attributes for structure
311 if iters>0; S.iters_from_file=iters; end
312 S.attributes.global=ncgetatt(nc,'global');
313 [pstr,netcdf_fname,ext] = fileparts(fname);
314 if strcmp(netcdf_fname(end-3:end),'glob')
315 % assume it is a global file produced by gluemnc and change some
316 % attributes
317 S.attributes.global.sNx = S.attributes.global.Nx;
318 S.attributes.global.sNy = S.attributes.global.Ny;
319 S.attributes.global.nPx = 1;
320 S.attributes.global.nSx = 1;
321 S.attributes.global.nPy = 1;
322 S.attributes.global.nSy = 1;
323 S.attributes.global.tile_number = 1;
324 S.attributes.global.nco_openmp_thread_number = 1;
325 end
326
327 % Read variable data
328 for ivar=1:size(varlist,2)
329
330 cvar=char(varlist{ivar});
331 varid=ncfindvarid(nc,cvar);
332 if isempty(varid)
333 disp(['No such variable ''',cvar,''' in MNC file ',fname]);
334 continue
335 end
336
337 [varname,xtype,dimids,natts] = netcdf.inqVar(nc,varid);
338
339 % Dimensions
340 clear sizVar dims
341 for k=1:length(dimids)
342 [dims{k}, sizVar(k)] = netcdf.inqDim(nc,dimids(k));
343 end
344 if length(sizVar)==1; sizVar(2) = 1; end;
345 nDims=length(sizVar);
346 if dims{1} == 'T'
347 if isempty(find(fii)), error('Iters not found'); end
348 it = length(dims);
349 tmpdata = netcdf.getVar(nc,varid,'double');
350 tmpdata = tmpdata(fii,:);
351 % leading unity dimensions get lost; add them back:
352 tmpdata=reshape(tmpdata,[length(fii) sizVar(2:end)]);
353 else
354 it = 0;
355 tmpdata = netcdf.getVar(nc,varid,'double');
356 % leading unity dimensions get lost; add them back:
357 tmpdata=reshape(tmpdata,sizVar);
358 end
359
360 if dBug > 1,
361 fprintf([' var:',cvar,': nDims=%i ('],nDims);fprintf(' %i',sizVar);
362 fprintf('):%i,nD=%i,it=%i ;',length(size(tmpdata)),length(dims),it);
363 end
364 [ni nj nk nm nn no np]=size(tmpdata);
365
366 [i0,j0,fn]=findTileOffset(S);
367 cdim=dims{1}; if cdim(1)~='X'; i0=0; end
368 cdim=dims{1}; if cdim(1)=='Y'; i0=j0; j0=0; end
369 if length(dims)>1;
370 cdim=dims{2}; if cdim(1)~='Y'; j0=0; end
371 else
372 j0=0;
373 end
374 if dBug > 1,
375 fprintf(' i0,ni= %i %i; j,nj= %i %i; nk=%i :',i0,ni,j0,nj,nk);
376 end
377
378 Sstr = '';
379 for istr = 1:max(nDims,length(dims)),
380 if istr == it, Sstr = [Sstr,'dii,'];
381 elseif istr == 1, Sstr = [Sstr,'i0+(1:ni),'];
382 elseif istr == 2, Sstr = [Sstr,'j0+(1:nj),'];
383 elseif istr == 3, Sstr = [Sstr,'(1:nk),'];
384 elseif istr == 4, Sstr = [Sstr,'(1:nm),'];
385 elseif istr == 5, Sstr = [Sstr,'(1:nn),'];
386 elseif istr == 6, Sstr = [Sstr,'(1:no),'];
387 elseif istr == 7, Sstr = [Sstr,'(1:np),'];
388 else, error('Can''t handle this many dimensions!');
389 end
390 end
391 eval(['S.(cvar)(',Sstr(1:end-1),')=tmpdata;'])
392 %S.(cvar)(i0+(1:ni),j0+(1:nj),(1:nk),(1:nm),(1:nn),(1:no),(1:np))=tmpdata;
393 if dBug > 1, fprintf(' %i',size(S.(cvar))); fprintf('\n'); end
394
395 S.attributes.(cvar)=ncgetatt(nc,cvar);
396 % replace missing or FillValues with NaN
397 if isfield(S.attributes.(cvar),'missing_value');
398 misval = S.attributes.(cvar).missing_value;
399 S.(cvar)(S.(cvar) == misval) = NaN;
400 end
401 if isfield(S.attributes.(cvar),'FillValue_');
402 misval = S.attributes.(cvar).FillValue_;
403 S.(cvar)(S.(cvar) == misval) = NaN;
404 end
405
406 end % for ivar
407
408 if isempty(S)
409 error('Something didn''t work!!!');
410 end
411
412 return
413
414 function vf = ncgetvar(nc,varname)
415 % read a netcdf variable
416
417 % find out basics about the files
418 [ndims nvars natts dimm] = netcdf.inq(nc);
419 vf = [];
420 varid = [];
421 for k=0:nvars-1
422 if strcmp(netcdf.inqVar(nc,k),varname)
423 varid = netcdf.inqVarId(nc,varname);
424 end
425 end
426 if ~isempty(varid);
427 [varn,xtype,dimids,natts] = netcdf.inqVar(nc,varid);
428 % get data
429 vf = double(netcdf.getVar(nc,varid));
430 else
431 % do nothing
432 end
433
434 return
435
436 function misval = ncgetmisval(nc,varid)
437
438 [varname,xtype,dimids,natts] = netcdf.inqVar(nc,varid);
439 misval = [];
440 for k=0:natts-1
441 attn = netcdf.inqAttName(nc,varid,k);
442 if strcmp(attn,'missing_value') | strcmp(attn,'_FillValue')
443 misval = double(netcdf.getAtt(nc,varid,attname));
444 end
445 end
446
447 function A = ncgetatt(nc,varname)
448 % get all attributes and put them into a struct
449
450 % 1. get global properties of file
451 [ndims nvars natts dimm] = netcdf.inq(nc);
452
453 % get variable ID and properties
454 if strcmp('global',varname)
455 % "variable" is global
456 varid = netcdf.getConstant('NC_GLOBAL');
457 else
458 % find variable ID and properties
459 varid = ncfindvarid(nc,varname);
460 if ~isempty(varid)
461 [varn,xtype,dimids,natts] = netcdf.inqVar(nc,varid);
462 else
463 warning(sprintf('variable %s not found',varname))
464 end
465 end
466
467 if natts > 1
468 for k=0:natts-1
469 attn = netcdf.inqAttName(nc,varid,k);
470 [xtype attlen]=netcdf.inqAtt(nc,varid,attn);
471 attval = netcdf.getAtt(nc,varid,attn);
472 if ~ischar(attval)
473 attval = double(attval);
474 end
475 A.(char(attn))=attval;
476 end
477 else
478 A = 'none';
479 end
480
481 return
482
483
484 function varid = ncfindvarid(nc,varname)
485
486 [ndims nvars natts dimm] = netcdf.inq(nc);
487 varid=[];
488 for k=0:nvars-1
489 if strcmp(netcdf.inqVar(nc,k),varname);
490 varid = netcdf.inqVarId(nc,varname);
491 break
492 end
493 end
494
495 return

  ViewVC Help
Powered by ViewVC 1.1.22