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:  $
% $Name:  $

% Initializations
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 arguement']);
        end
    end
end

% Create list of filenames
for eachfile=file
        filepathtemp=eachfile{:};
        indecies = find(filepathtemp=='/');
        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}];
    nc=netcdf(char(eachfile),'read');
    S=rdmnc_local(nc,varlist,iters,S);
    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)

    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
    
    % No variables specified? Default to all
    if isempty(varlist), varlist=ncnames(var(nc)); end
    
    % Attributes for structure
    if iters>0; S.iters_read_from_file=iters; end
    S.attributes.global=read_att(nc);
    
        % 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
        adj = 0;
        if dims{1} == 'T'
            
            if isempty(find(fii)), disp('Iters not found'); return, end
            
            tmpdata = nc{cvar}(fii,:);
            if ismember('Zd000001' ,dims), adj = adj - 1; end
            if ismember('Zmd000001',dims), adj = adj - 1; end
            it = length(dims) + adj;
        else
            tmpdata = nc{cvar}(:);
            it = 0;
        end
        
        tmpdata=shiftdim(permute(tmpdata,[9:-1:1]));
        [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
        
        Sstr = '';
        for istr = 1:7
            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;
        S.attributes.(cvar)=read_att(nc{cvar});
        end

if isempty(S)
    error('Something didn''t work!!!');
end
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: $
34	% $Name: $
35
36	% Initializations
37	file={};
38	filepaths={};
39	files={};
40	varlist={};
41	iters=[];
42
43	% Process function arguments
44	for iarg=1:nargin;
45	arg=varargin{iarg};
46	if ischar(arg)
47	if isempty(dir(char(arg)))
48	varlist{end+1}=arg;
49	else
50	file{end+1}=arg;
51	end
52	else
53	if isempty(iters)
54	iters=arg;
55	else
56	error(['The only allowed numeric argument is iterations',...
57	' to read in as a vector for the last arguement']);
58	end
59	end
60	end
61
62	% Create list of filenames
63	for eachfile=file
64	filepathtemp=eachfile{:};
65	indecies = find(filepathtemp=='/');
66	if ~isempty(indecies)
67	filepathtemp = filepathtemp(1:indecies(end));
68	else
69	filepathtemp = '';
70	end
71	expandedEachFile=dir(char(eachfile{1}));
72	for i=1:size(expandedEachFile,1);
73	if expandedEachFile(i).isdir==0
74	files{end+1}=expandedEachFile(i).name;
75	filepaths{end+1}=filepathtemp;
76	end
77	end
78	end
79
80
81	% If iterations unspecified, open all the files and make list of all the
82	% iterations that appear, use this iterations list for data extraction.
83	if isempty(iters)
84	iters = [];
85	for ieachfile=1:length(files)
86	eachfile = [filepaths{ieachfile},files{ieachfile}];
87	nc=netcdf(char(eachfile),'read');
88	nciters = nc{'iter'}(:);
89	if isempty(nciters), nciters = nc{'T'}(:); end
90	iters = [iters,nciters];
91	close(nc);
92	end
93	iters = unique(iters);
94	end
95
96	% Cycle through files for data extraction.
97	S.attributes=[];
98	for ieachfile=1:length(files)
99	eachfile = [filepaths{ieachfile},files{ieachfile}];
100	nc=netcdf(char(eachfile),'read');
101	S=rdmnc_local(nc,varlist,iters,S);
102	close(nc);
103	end
104
105
106	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
107	% Local functions %
108	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
109
110	function [A] = read_att(nc);
111	allatt=ncnames(att(nc));
112	if ~isempty(allatt)
113	for attr=allatt;
114	A.(char(attr))=nc.(char(attr))(:);
115	end
116	else
117	A = 'none';
118	end
119
120	function [i0,j0,fn] = findTileOffset(S);
121	fn=0;
122	if isfield(S,'attributes') & isfield(S.attributes,'global')
123	G=S.attributes.global;
124	tn=G.tile_number;
125	snx=G.sNx; sny=G.sNy; nsx=G.nSx; nsy=G.nSy; npx=G.nPx; npy=G.nPy;
126	ntx=nsxnpx;nty=nsynpy;
127	gbi=mod(tn-1,ntx); gbj=(tn-gbi-1)/ntx;
128	i0=snxgbi; j0=snygbj;
129	if isfield(S.attributes.global,'exch2_myFace')
130	fn=G.exch2_myFace;
131	end
132	else
133	i0=0;j0=0;
134	end
135	%[snx,sny,nsx,nsy,npx,npy,ntx,nty,i0,j0,fn];
136
137	function [S] = rdmnc_local(nc,varlist,iters,S)
138
139	fiter = nc{'iter'}(:); % File iterations present
140	if isempty(fiter), fiter = nc{'T'}(:); end
141	if isinf(iters); iters = fiter(end); end
142	if isnan(iters); iters = fiter; end
143	[fii,dii] = ismember(fiter,iters); fii = find(fii); % File iteration index
144	dii = dii(find(dii ~= 0)); % Data interation index
145
146	% No variables specified? Default to all
147	if isempty(varlist), varlist=ncnames(var(nc)); end
148
149	% Attributes for structure
150	if iters>0; S.iters_read_from_file=iters; end
151	S.attributes.global=read_att(nc);
152
153	% Read variable data
154	for ivar=1:size(varlist,2)
155
156	cvar=char(varlist{ivar});
157	if isempty(nc{cvar})
158	disp(['No such variable ''',cvar,''' in MNC file ',name(nc)]);
159	continue
160	end
161
162	dims = ncnames(dim(nc{cvar})); % Dimensions
163	adj = 0;
164	if dims{1} == 'T'
165
166	if isempty(find(fii)), disp('Iters not found'); return, end
167
168	tmpdata = nc{cvar}(fii,:);
169	if ismember('Zd000001' ,dims), adj = adj - 1; end
170	if ismember('Zmd000001',dims), adj = adj - 1; end
171	it = length(dims) + adj;
172	else
173	tmpdata = nc{cvar}(:);
174	it = 0;
175	end
176
177	tmpdata=shiftdim(permute(tmpdata,[9:-1:1]));
178	[ni nj nk nm nn no np]=size(tmpdata);
179
180	[i0,j0,fn]=findTileOffset(S);
181	cdim=dims{end}; if cdim(1)~='X'; i0=0; end
182	cdim=dims{end}; if cdim(1)=='Y'; i0=j0; j0=0; end
183	if length(dims)>1;
184	cdim=dims{end-1}; if cdim(1)~='Y'; j0=0; end
185	else
186	j0=0;
187	end
188
189	Sstr = '';
190	for istr = 1:7
191	if istr == it, Sstr = [Sstr,'dii,'];
192	elseif istr == 1, Sstr = [Sstr,'i0+(1:ni),'];
193	elseif istr == 2, Sstr = [Sstr,'j0+(1:nj),'];
194	elseif istr == 3, Sstr = [Sstr,'(1:nk),'];
195	elseif istr == 4, Sstr = [Sstr,'(1:nm),'];
196	elseif istr == 5, Sstr = [Sstr,'(1:nn),'];
197	elseif istr == 6, Sstr = [Sstr,'(1:no),'];
198	elseif istr == 7, Sstr = [Sstr,'(1:np),'];
199	else, error('Can''t handle this many dimensions!');
200	end
201	end
202
203	eval(['S.(cvar)(',Sstr(1:end-1),')=tmpdata;'])
204	%S.(cvar)(i0+(1:ni),j0+(1:nj),(1:nk),(1:nm),(1:nn),(1:no),(1:np))=tmpdata;
205	S.attributes.(cvar)=read_att(nc{cvar});
206	end
207
208	if isempty(S)
209	error('Something didn''t work!!!');
210	end