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