/[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.18 - (show annotations) (download)
Fri Nov 13 07:57:16 2009 UTC (14 years, 8 months ago) by mlosch
Branch: MAIN
CVS Tags: checkpoint61y
Changes since 1.17: +10 -0 lines
let rdmnc_local find missing_value and replace it with NaN

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.16 2007/03/07 09:59:53 mlosch Exp $
34 % $Name: $
35
36 % Initializations
37 dBug=0;
38 file={};
39 filepaths={};
40 files={};
41 varlist={};
42 iters=[];
43
44 % Process function arguments
45 for iarg=1:nargin;
46 arg=varargin{iarg};
47 if ischar(arg)
48 if isempty(dir(char(arg)))
49 varlist{end+1}=arg;
50 else
51 file{end+1}=arg;
52 end
53 else
54 if isempty(iters)
55 iters=arg;
56 else
57 error(['The only allowed numeric argument is iterations',...
58 ' to read in as a vector for the last argument']);
59 end
60 end
61 end
62 if isempty(file)
63 if isempty(varlist),
64 fprintf( 'No file name in argument list\n');
65 else
66 fprintf(['No file in argument list:\n ==> ',char(varlist(1))]);
67 for i=2:size(varlist,2), fprintf([' , ',char(varlist(i))]); end
68 fprintf(' <==\n');
69 end
70 error(' check argument list !!!');
71 end
72
73 % Create list of filenames
74 for eachfile=file
75 filepathtemp=eachfile{:};
76 indecies = find(filepathtemp==filesep);
77 if ~isempty(indecies)
78 filepathtemp = filepathtemp(1:indecies(end));
79 else
80 filepathtemp = '';
81 end
82 expandedEachFile=dir(char(eachfile{1}));
83 for i=1:size(expandedEachFile,1);
84 if expandedEachFile(i).isdir==0
85 files{end+1}=expandedEachFile(i).name;
86 filepaths{end+1}=filepathtemp;
87 end
88 end
89 end
90
91
92 % If iterations unspecified, open all the files and make list of all the
93 % iterations that appear, use this iterations list for data extraction.
94 if isempty(iters)
95 iters = [];
96 for ieachfile=1:length(files)
97 eachfile = [filepaths{ieachfile},files{ieachfile}];
98 nc=netcdf(char(eachfile),'read');
99 nciters = nc{'iter'}(:);
100 if isempty(nciters), nciters = nc{'T'}(:); end
101 iters = [iters,nciters'];
102 close(nc);
103 end
104 iters = unique(iters');
105 end
106
107 % Cycle through files for data extraction.
108 S.attributes=[];
109 for ieachfile=1:length(files)
110 eachfile = [filepaths{ieachfile},files{ieachfile}];
111 if dBug > 0, fprintf([' open: ',eachfile]); end
112 nc=netcdf(char(eachfile),'read');
113 S=rdmnc_local(nc,varlist,iters,S,dBug);
114 close(nc);
115 end
116
117
118 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
119 % Local functions %
120 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
121
122 function [A] = read_att(nc);
123 allatt=ncnames(att(nc));
124 if ~isempty(allatt)
125 for attr=allatt;
126 A.(char(attr))=nc.(char(attr))(:);
127 end
128 else
129 A = 'none';
130 end
131
132 function [i0,j0,fn] = findTileOffset(S);
133 fn=0;
134 if isfield(S,'attributes') & isfield(S.attributes,'global')
135 G=S.attributes.global;
136 tn=G.tile_number;
137 snx=G.sNx; sny=G.sNy; nsx=G.nSx; nsy=G.nSy; npx=G.nPx; npy=G.nPy;
138 ntx=nsx*npx;nty=nsy*npy;
139 gbi=mod(tn-1,ntx); gbj=(tn-gbi-1)/ntx;
140 i0=snx*gbi; j0=sny*gbj;
141 if isfield(S.attributes.global,'exch2_myFace')
142 fn=G.exch2_myFace;
143 end
144 else
145 i0=0;j0=0;
146 end
147 %[snx,sny,nsx,nsy,npx,npy,ntx,nty,i0,j0,fn];
148
149 function [S] = rdmnc_local(nc,varlist,iters,S,dBug)
150
151 fiter = nc{'iter'}(:); % File iterations present
152 if isempty(fiter), fiter = nc{'T'}(:); end
153 if isinf(iters); iters = fiter(end); end
154 if isnan(iters); iters = fiter; end
155 [fii,dii] = ismember(fiter,iters); fii = find(fii); % File iteration index
156 dii = dii(find(dii ~= 0)); % Data interation index
157 if dBug > 0,
158 fprintf(' ; fii='); fprintf(' %i',fii);
159 fprintf(' ; dii='); fprintf(' %i',dii); fprintf(' \n');
160 end
161
162 % No variables specified? Default to all
163 if isempty(varlist), varlist=ncnames(var(nc)); end
164
165 % Attributes for structure
166 if iters>0; S.iters_from_file=iters; end
167 S.attributes.global=read_att(nc);
168 [pstr,netcdf_fname,ext] = fileparts(name(nc));
169 if strcmp(netcdf_fname(end-3:end),'glob')
170 % assume it is a global file produced by gluemnc and change some
171 % attributes
172 S.attributes.global.sNx = S.attributes.global.Nx;
173 S.attributes.global.sNy = S.attributes.global.Ny;
174 S.attributes.global.nPx = 1;
175 S.attributes.global.nSx = 1;
176 S.attributes.global.nPy = 1;
177 S.attributes.global.nSy = 1;
178 S.attributes.global.tile_number = 1;
179 S.attributes.global.nco_openmp_thread_number = 1;
180 end
181
182 % Read variable data
183 for ivar=1:size(varlist,2)
184
185 cvar=char(varlist{ivar});
186 if isempty(nc{cvar})
187 disp(['No such variable ''',cvar,''' in MNC file ',name(nc)]);
188 continue
189 end
190
191 dims = ncnames(dim(nc{cvar})); % Dimensions
192 sizVar = size(nc{cvar}); nDims=length(sizVar);
193 if dims{1} == 'T'
194 if isempty(find(fii)), error('Iters not found'); end
195 it = length(dims);
196 tmpdata = nc{cvar}(fii,:);
197 %- leading unity dimensions get lost; add them back:
198 tmpdata=reshape(tmpdata,[length(fii) sizVar(2:end)]);
199 else
200 it = 0;
201 tmpdata = nc{cvar}(:);
202 %- leading unity dimensions get lost; add them back:
203 tmpdata=reshape(tmpdata,sizVar);
204 end
205
206 if dBug > 1,
207 fprintf([' var:',cvar,': nDims=%i ('],nDims);fprintf(' %i',size(nc{cvar}));
208 fprintf('):%i,nD=%i,it=%i ;',length(size(tmpdata)),length(dims),it);
209 end
210 if length(dims) > 1,
211 tmpdata=permute(tmpdata,[nDims:-1:1]);
212 end
213 if dBug > 1,
214 % fprintf('(tmpdata:');fprintf(' %i',size(tmpdata)); fprintf(')');
215 end
216 [ni nj nk nm nn no np]=size(tmpdata);
217
218 [i0,j0,fn]=findTileOffset(S);
219 cdim=dims{end}; if cdim(1)~='X'; i0=0; end
220 cdim=dims{end}; if cdim(1)=='Y'; i0=j0; j0=0; end
221 if length(dims)>1;
222 cdim=dims{end-1}; if cdim(1)~='Y'; j0=0; end
223 else
224 j0=0;
225 end
226 if dBug > 1,
227 fprintf(' i0,ni= %i %i; j,nj= %i %i; nk=%i :',i0,ni,j0,nj,nk);
228 end
229
230 Sstr = '';
231 for istr = 1:max(nDims,length(dims)),
232 if istr == it, Sstr = [Sstr,'dii,'];
233 elseif istr == 1, Sstr = [Sstr,'i0+(1:ni),'];
234 elseif istr == 2, Sstr = [Sstr,'j0+(1:nj),'];
235 elseif istr == 3, Sstr = [Sstr,'(1:nk),'];
236 elseif istr == 4, Sstr = [Sstr,'(1:nm),'];
237 elseif istr == 5, Sstr = [Sstr,'(1:nn),'];
238 elseif istr == 6, Sstr = [Sstr,'(1:no),'];
239 elseif istr == 7, Sstr = [Sstr,'(1:np),'];
240 else, error('Can''t handle this many dimensions!');
241 end
242 end
243 eval(['S.(cvar)(',Sstr(1:end-1),')=tmpdata;'])
244 %S.(cvar)(i0+(1:ni),j0+(1:nj),(1:nk),(1:nm),(1:nn),(1:no),(1:np))=tmpdata;
245 if dBug > 1, fprintf(' %i',size(S.(cvar))); fprintf('\n'); end
246
247 S.attributes.(cvar)=read_att(nc{cvar});
248 % replace missing or FillValues with NaN
249 attnames=fieldnames(S.attributes.(cvar));
250 if ~isempty(attnames)
251 for k=1:length(attnames)
252 if strcmp(attnames{k},'missing_value') ...
253 | strcmp(attnames{k},'FillValue_')
254 S.(cvar)(S.(cvar) == S.attributes.(cvar).(attnames{k})) = NaN;
255 end
256 end
257 end
258 end
259
260 if isempty(S)
261 error('Something didn''t work!!!');
262 end

  ViewVC Help
Powered by ViewVC 1.1.22