utils/matlab/rdmnc.m

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

% $Header: /u/gcmpack/MITgcm/utils/matlab/rdmnc.m,v 1.16 2007/03/07 09:59:53 mlosch Exp $
% $Name:  $

% Initializations
dBug=0;
file={};
filepaths={};
files={};
varlist={};
iters=[];

% 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}];
        nc=netcdf(char(eachfile),'read');
        nciters = nc{'iter'}(:);
        if isempty(nciters), nciters = nc{'T'}(:); end
        iters = [iters,nciters'];
        close(nc);
    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
    nc=netcdf(char(eachfile),'read');
    S=rdmnc_local(nc,varlist,iters,S,dBug);
    close(nc);
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                             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;
        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
        
        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
        
        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});
        end

if isempty(S)
    error('Something didn''t work!!!');
end
1	enderton	1.6	function [S] = rdmnc(varargin)
2	enderton	1.3
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	adcroft	1.1	%
16	enderton	1.3	% Output:
17			% S Structure with fields corresponding to 'VAR1', 'VAR2', ...
18	adcroft	1.1	%
19	enderton	1.3	% 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	adcroft	1.1	%
26	enderton	1.3	% Example:
27			% >> S=rdmnd('state.*','XC','YC','T');
28			% >> imagesc( S.XC, S.YC, S.T(:,:,1)' );
29	adcroft	1.1	%
30	enderton	1.4	% Author: Alistair Adcroft
31	enderton	1.3	% Modifications: Daniel Enderton
32	adcroft	1.1
33	mlosch	1.17	% $Header: /u/gcmpack/MITgcm/utils/matlab/rdmnc.m,v 1.16 2007/03/07 09:59:53 mlosch Exp $
34	jmc	1.12	% $Name: $
35
36	enderton	1.3	% Initializations
37	jmc	1.13	dBug=0;
38	adcroft	1.1	file={};
39	enderton	1.3	filepaths={};
40	adcroft	1.1	files={};
41			varlist={};
42	enderton	1.3	iters=[];
43	adcroft	1.1
44			% Process function arguments
45			for iarg=1:nargin;
46	enderton	1.3	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	jmc	1.14	' to read in as a vector for the last argument']);
59	enderton	1.3	end
60			end
61	adcroft	1.1	end
62	jmc	1.14	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	adcroft	1.1
73			% Create list of filenames
74			for eachfile=file
75	enderton	1.3	filepathtemp=eachfile{:};
76	mlosch	1.17	indecies = find(filepathtemp==filesep);
77	enderton	1.3	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	adcroft	1.1	end
90
91
92	enderton	1.3	% 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	enderton	1.4	nciters = nc{'iter'}(:);
100			if isempty(nciters), nciters = nc{'T'}(:); end
101	jmc	1.15	iters = [iters,nciters'];
102	enderton	1.3	close(nc);
103			end
104	jmc	1.15	iters = unique(iters');
105	adcroft	1.1	end
106
107	enderton	1.3	% Cycle through files for data extraction.
108			S.attributes=[];
109			for ieachfile=1:length(files)
110			eachfile = [filepaths{ieachfile},files{ieachfile}];
111	jmc	1.13	if dBug > 0, fprintf([' open: ',eachfile]); end
112	enderton	1.3	nc=netcdf(char(eachfile),'read');
113	jmc	1.13	S=rdmnc_local(nc,varlist,iters,S,dBug);
114	enderton	1.3	close(nc);
115	adcroft	1.1	end
116
117	enderton	1.5
118	enderton	1.3	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
119			% Local functions %
120			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
121	adcroft	1.1
122	enderton	1.3	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	adcroft	1.1
132	enderton	1.3	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=nsxnpx;nty=nsynpy;
139			gbi=mod(tn-1,ntx); gbj=(tn-gbi-1)/ntx;
140			i0=snxgbi; j0=snygbj;
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	jmc	1.13	function [S] = rdmnc_local(nc,varlist,iters,S,dBug)
150	enderton	1.3
151			fiter = nc{'iter'}(:); % File iterations present
152	enderton	1.4	if isempty(fiter), fiter = nc{'T'}(:); end
153	mlosch	1.10	if isinf(iters); iters = fiter(end); end
154			if isnan(iters); iters = fiter; end
155	enderton	1.3	[fii,dii] = ismember(fiter,iters); fii = find(fii); % File iteration index
156			dii = dii(find(dii ~= 0)); % Data interation index
157	jmc	1.13	if dBug > 0,
158			fprintf(' ; fii='); fprintf(' %i',fii);
159			fprintf(' ; dii='); fprintf(' %i',dii); fprintf(' \n');
160			end
161	enderton	1.4
162	enderton	1.3	% No variables specified? Default to all
163			if isempty(varlist), varlist=ncnames(var(nc)); end
164	enderton	1.4
165	enderton	1.3	% Attributes for structure
166	jmc	1.14	if iters>0; S.iters_from_file=iters; end
167	enderton	1.3	S.attributes.global=read_att(nc);
168	mlosch	1.16	[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	enderton	1.4
182	enderton	1.3	% Read variable data
183	enderton	1.4	for ivar=1:size(varlist,2)
184
185	enderton	1.3	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	enderton	1.6	dims = ncnames(dim(nc{cvar})); % Dimensions
192	jmc	1.14	sizVar = size(nc{cvar}); nDims=length(sizVar);
193	enderton	1.6	if dims{1} == 'T'
194	jmc	1.14	if isempty(find(fii)), error('Iters not found'); end
195			it = length(dims);
196	enderton	1.7	tmpdata = nc{cvar}(fii,:);
197	jmc	1.14	%- leading unity dimensions get lost; add them back:
198			tmpdata=reshape(tmpdata,[length(fii) sizVar(2:end)]);
199	enderton	1.3	else
200	jmc	1.14	it = 0;
201	enderton	1.6	tmpdata = nc{cvar}(:);
202	jmc	1.14	%- leading unity dimensions get lost; add them back:
203			tmpdata=reshape(tmpdata,sizVar);
204	enderton	1.3	end
205	enderton	1.4
206	jmc	1.13	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	enderton	1.3	[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	jmc	1.13	if dBug > 1,
227			fprintf(' i0,ni= %i %i; j,nj= %i %i; nk=%i :',i0,ni,j0,nj,nk);
228			end
229	enderton	1.3
230			Sstr = '';
231	jmc	1.13	for istr = 1:max(nDims,length(dims)),
232	enderton	1.3	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	jmc	1.13	if dBug > 1, fprintf(' %i',size(S.(cvar))); fprintf('\n'); end
246
247	enderton	1.3	S.attributes.(cvar)=read_att(nc{cvar});
248			end
249	adcroft	1.1
250			if isempty(S)
251	enderton	1.3	error('Something didn''t work!!!');
252	adcroft	1.1	end