enderton/Diagnostics/DiagLoadGradsData.m

function [data,xax,yax,zax,months,time,dim] = ...
    DiagLoadGradsData(Grads,dad,fln,DiagDebug)


% Initial assumptions about data.
ShiftData = 0;
MultFiles = 0;
format = 'NONSEQUENTIAL';


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                          Parse table file                               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tablfile = [dad,'/',Grads];
file = textread(tablfile,'%s','delimiter','\n','whitespace','');

% Cycle though lines of the GRADS table file, pasring the file line-by-line
% to properly read in and format the data.  The first line of every file
% has an identifier denoting the information on that line.  The following
% blocks of code parse the information according to different identifiers.
for iline = 1:length(file)
    rem = file{iline}; tokens = {}; itoken = 0;
    while isstr(rem)
        [token,rem] = strtok(rem); itoken = itoken + 1;
        if ~isempty(token), tokens{itoken} = token; end
    end
    
    % Determine file(s) to read in.
    if isequal(tokens{1},'DSET')
        if isempty(findstr(tokens{2},'%'))
            datafile = [dad,tokens{2}(2:end)];
        else
            str = tokens{2}(2:end);
            alldatafile = strrep(str,'%y2','*');
            alldatafile = strrep(alldatafile,'%m2','*');
            files = ls([dad,alldatafile]);
            breaks = [0,find(isspace(files))];
            for ibreak = 1:length(breaks)-1
                datafile{ibreak} = files([breaks(ibreak)+1:breaks(ibreak+1)-1]);
            end
            MultFiles = 1;
        end
    end  
    
    % Read in x-axis information.
    if isequal(tokens{1},'XDEF')
        if ~isequal(tokens{3},'LINEAR')
            error('Unable to deal with nonlinear grads x-axis data!');
        end
        num = str2num(tokens{2});
        ini = str2num(tokens{4});
        inc = str2num(tokens{5});
        xax = [ini:inc:ini+(num-1)*inc];
        % If the x-axis is from 0 to 360, the axis information and the data
        % will later have to be shifted to a -180 to 180 axis.
        if min(xax) >= 0 && max(xax) > 180
            ShiftData = 1;
        end
        nx = length(xax);
    end
    
    % Read in y-axis information.
    if isequal(tokens{1},'YDEF')
        if ~isequal(tokens{3},'LINEAR')
            error('Unable to deal with nonlinear grads y-axis data!');
        end
        num = str2num(tokens{2});
        ini = str2num(tokens{4});
        inc = str2num(tokens{5});
        yax = [ini:inc:ini+(num-1)*inc]; ny = length(yax);
    end
    
    % Read in z-axis information.
    if isequal(tokens{1},'ZDEF')
        if isequal(tokens{2},'1')
            zax = [0]; zax_found = 1; dim = 2;
        else
            dim = 3;
            if isequal(tokens{3},'LINEAR')
                num = str2num(tokens{2});
                ini = str2num(tokens{4});
                inc = str2num(tokens{5});
                zax = [ini:inc:ini+(num-1)*inc];
            elseif isequal(tokens{3},'LEVELS')
                num = str2num(tokens{2});
                for inum = 1:num
                    zax(inum) = 100.*str2num(tokens{inum+3});
                end
            else
                error('Unknown information in z-axis data!');
            end
        end
        nz = length(zax);
    end
    
    % Read in t-axis (time) information.
    if isequal(tokens{1},'TDEF'), ini=tokens{4};
        if ~isequal(tokens{3},'LINEAR')
            error('Currently unable to deal with nonlinear grads t-axis data!');
        end
        if ~isequal(tokens{5},'1mo') & ~isequal(tokens{5},'01mo')
            % Note that monthly mean assumptions are also made when data
            % spanning multiple files are read in later on and that this
            % part must be updated if non-monthly mean data is allowed.
            error('Currently unable to deal with non monthly mean grads data!');
        end
        % Determine initial month and initial year.  This is done in an ad
        % hoc manner in that depending on the length start time string, the
        % month is pluck out.  There would be problems is there are start
        % time strings of different lengths have months in different
        % locations.  In this case the month will not be found and an error
        % will be thrown.  The year is only required for data spanning
        % multiple files, which thus fas has a start time string length of
        % 12.  If it is not found, and error of undefined yrchar will be
        % cast and the code will need to be appropriately updated.
        if     length(ini) == 13, monchar=ini(9:11);
        elseif length(ini) == 12, monchar=ini(6:8);  yrchar = ini(9:12);
        elseif length(ini) == 10, monchar=ini(6:8);
        elseif length(ini) == 7,  monchar=ini(1:3);
        else, error('Cannot parse TDEF correctly'); end
        monchar = upper(monchar);
        if     isequal(monchar,'JAN'), inimonth = 1;
        elseif isequal(monchar,'FEB'), inimonth = 2;
        elseif isequal(monchar,'MAR'), inimonth = 3;
        elseif isequal(monchar,'APR'), inimonth = 4;
        elseif isequal(monchar,'MAY'), inimonth = 5;
        elseif isequal(monchar,'JUN'), inimonth = 6;
        elseif isequal(monchar,'JUL'), inimonth = 7;
        elseif isequal(monchar,'AUG'), inimonth = 8;
        elseif isequal(monchar,'SEP'), inimonth = 9;
        elseif isequal(monchar,'OCT'), inimonth = 10;
        elseif isequal(monchar,'NOV'), inimonth = 11;
        elseif isequal(monchar,'DEC'), inimonth = 12;
        else, error('Can''t determine initial month in GRADS data!'); end
        num = str2num(tokens{2});
        % Determine num and inimonth parameters accounting for the unknown
        % number of data files present when dealing with multiple data
        % files.  Assumption of monthly mean data is made.
        if MultFiles;
            initFileFound = 0;
            for ifile = 1:length(datafile)
                initFileFound = ~isempty(findstr(yrchar,datafile{ifile})) ...
                                | initFileFound;
            end
            if ~initFileFound, inimonth = 1; end
            num = 12.*length(datafile);
        end
        % Assumption of monthly mean data is made.
        months=[inimonth:num+inimonth-1];
        time=months/12; nt = length(months);
    end
    
    % Find where the desired variable is located.
    if isequal(tokens{1},'VARS')
        nv = str2num(tokens{2});
        VARSline = iline;
    end
    if isequal(tokens{1},fln)
        FIELDline = iline;
        ivar = iline - VARSline;
    end
    
    % Find the values defonting undefined data.
    if isequal(tokens{1},'UNDEF')
        undef = str2num(tokens{2});
    end
    
    % Determine data format type.
    if isequal(tokens{1},'FORMAT')
        if isequal(tokens{2},'SEQUENTIAL') | isequal(tokens{2},'sequential')
            format = 'SEQUENTIAL';
        else
            disp(['Unrecognized grads FORMAT:  ',tokens{2}]);
        end
    end 
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                       Verification and debugging                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Verify that everything we need from the table file has been read in.
if DiagDebug, disp(['  Debug -- grads data format:  ',format]); end
parameters = {'xax' ,'yax' ,'zax' ,'time','nv'  ,'undef','ivar'};
GRADSnames = {'XDEF','YDEF','ZDEF','TDEF','VARS','UNDEF',fln   };
for ii = 1:length(parameters)
    try
        eval([parameters{ii},';']);
    catch
        error(['GRADS information not found:  ',GRADSnames{ii}]);
    end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                              Read data file                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Read in data from a single file.
if ~MultFiles
    fid=fopen(datafile,'r','b');
    data = fread(fid,'real*4');
    fclose(fid);
        if isequal(format,'SEQUENTIAL') | isequal(format,'sequential')
        index=true([nx*ny*nz*nv*nt+2*nv*nt,1]);
        index([1:nx*ny*nz+2:end])=false;
        index([2:nx*ny*nz+2:end])=false;
        data = data(index);
        end
    data = reshape(data,[nx,ny,nz,nv,nt]);
    data = data(:,:,:,ivar,:);
    datatmp(:,:,:,:,inimonth:num+inimonth-1) = data;
    datatmp(:,:,:,:,1:inimonth-1) = NaN;
    data=datatmp;
% Read in data from multiple files.
else
    for ifile = 1:length(datafile)
        if DiagDebug, disp(['  Debug -- Reading grads ',...
                            'data file:  ',datafile{ifile}]); end
        fid=fopen(datafile{ifile},'r','b');
        datatemp = fread(fid,'real*4');
        fclose(fid);
        % Assumption of monthly mean data is made.
        if MultFiles, nt = 12; end
        if initFileFound & (inimonth ~= 1)
            error(['Currently not able to handle Grads data spanning ',...
                    'multiple files not starting in January']);
        end
        datatemp = reshape(datatemp,[nx,ny,nz,nv,nt]);
        datatemp = squeeze(datatemp(:,:,:,ivar,:));
        data(1:nx,1:ny,1:nz,nt.*(ifile-1)+1:nt.*ifile) = datatemp;
    end
end
data = squeeze(data);

% Turn undefined values into NaN.
tol = 1e-5;
data( abs((data-undef)/undef) < tol ) = NaN;

% Shift data such that it is on a -180 to 180 longitude axis.
if ShiftData
    indexWestHemi = xax>=180;
    indexEastHemi = xax<180;
    data = cat(1,data(indexWestHemi,:,:,:),data(indexEastHemi,:,:,:));
    xax = cat(2,xax(indexWestHemi)-360,xax(indexEastHemi));
end
1	enderton	1.1	function [data,xax,yax,zax,months,time,dim] = ...
2	enderton	1.2	DiagLoadGradsData(Grads,dad,fln,DiagDebug)
3
4
5	enderton	1.5	% Initial assumptions about data.
6			ShiftData = 0;
7			MultFiles = 0;
8			format = 'NONSEQUENTIAL';
9
10
11	enderton	1.2	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
12			% Parse table file %
13			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
14	molod	1.4	tablfile = [dad,'/',Grads];
15	enderton	1.1	file = textread(tablfile,'%s','delimiter','\n','whitespace','');
16	enderton	1.2
17	enderton	1.5	% Cycle though lines of the GRADS table file, pasring the file line-by-line
18			% to properly read in and format the data. The first line of every file
19			% has an identifier denoting the information on that line. The following
20			% blocks of code parse the information according to different identifiers.
21	enderton	1.1	for iline = 1:length(file)
22			rem = file{iline}; tokens = {}; itoken = 0;
23			while isstr(rem)
24			[token,rem] = strtok(rem); itoken = itoken + 1;
25			if ~isempty(token), tokens{itoken} = token; end
26			end
27
28	enderton	1.5	% Determine file(s) to read in.
29			if isequal(tokens{1},'DSET')
30			if isempty(findstr(tokens{2},'%'))
31			datafile = [dad,tokens{2}(2:end)];
32			else
33			str = tokens{2}(2:end);
34			alldatafile = strrep(str,'%y2','*');
35			alldatafile = strrep(alldatafile,'%m2','*');
36			files = ls([dad,alldatafile]);
37			breaks = [0,find(isspace(files))];
38			for ibreak = 1:length(breaks)-1
39			datafile{ibreak} = files([breaks(ibreak)+1:breaks(ibreak+1)-1]);
40			end
41			MultFiles = 1;
42			end
43			end
44
45			% Read in x-axis information.
46	enderton	1.1	if isequal(tokens{1},'XDEF')
47			if ~isequal(tokens{3},'LINEAR')
48			error('Unable to deal with nonlinear grads x-axis data!');
49			end
50			num = str2num(tokens{2});
51			ini = str2num(tokens{4});
52			inc = str2num(tokens{5});
53			xax = [ini:inc:ini+(num-1)*inc];
54	enderton	1.5	% If the x-axis is from 0 to 360, the axis information and the data
55			% will later have to be shifted to a -180 to 180 axis.
56	enderton	1.1	if min(xax) >= 0 && max(xax) > 180
57	molod	1.3	ShiftData = 1;
58	enderton	1.1	end
59			nx = length(xax);
60			end
61
62	enderton	1.5	% Read in y-axis information.
63	enderton	1.1	if isequal(tokens{1},'YDEF')
64			if ~isequal(tokens{3},'LINEAR')
65			error('Unable to deal with nonlinear grads y-axis data!');
66			end
67			num = str2num(tokens{2});
68			ini = str2num(tokens{4});
69			inc = str2num(tokens{5});
70			yax = [ini:inc:ini+(num-1)*inc]; ny = length(yax);
71			end
72
73	enderton	1.5	% Read in z-axis information.
74	enderton	1.1	if isequal(tokens{1},'ZDEF')
75			if isequal(tokens{2},'1')
76			zax = [0]; zax_found = 1; dim = 2;
77			else
78	enderton	1.5	dim = 3;
79			if isequal(tokens{3},'LINEAR')
80			num = str2num(tokens{2});
81			ini = str2num(tokens{4});
82			inc = str2num(tokens{5});
83			zax = [ini:inc:ini+(num-1)*inc];
84			elseif isequal(tokens{3},'LEVELS')
85			num = str2num(tokens{2});
86			for inum = 1:num
87			zax(inum) = 100.*str2num(tokens{inum+3});
88			end
89			else
90			error('Unknown information in z-axis data!');
91	enderton	1.1	end
92			end
93			nz = length(zax);
94			end
95
96	enderton	1.5	% Read in t-axis (time) information.
97	enderton	1.1	if isequal(tokens{1},'TDEF'), ini=tokens{4};
98			if ~isequal(tokens{3},'LINEAR')
99	enderton	1.5	error('Currently unable to deal with nonlinear grads t-axis data!');
100	enderton	1.1	end
101	enderton	1.5	if ~isequal(tokens{5},'1mo') & ~isequal(tokens{5},'01mo')
102			% Note that monthly mean assumptions are also made when data
103			% spanning multiple files are read in later on and that this
104			% part must be updated if non-monthly mean data is allowed.
105	enderton	1.1	error('Currently unable to deal with non monthly mean grads data!');
106			end
107	enderton	1.5	% Determine initial month and initial year. This is done in an ad
108			% hoc manner in that depending on the length start time string, the
109			% month is pluck out. There would be problems is there are start
110			% time strings of different lengths have months in different
111			% locations. In this case the month will not be found and an error
112			% will be thrown. The year is only required for data spanning
113			% multiple files, which thus fas has a start time string length of
114			% 12. If it is not found, and error of undefined yrchar will be
115			% cast and the code will need to be appropriately updated.
116			if length(ini) == 13, monchar=ini(9:11);
117			elseif length(ini) == 12, monchar=ini(6:8); yrchar = ini(9:12);
118			elseif length(ini) == 10, monchar=ini(6:8);
119			elseif length(ini) == 7, monchar=ini(1:3);
120			else, error('Cannot parse TDEF correctly'); end
121			monchar = upper(monchar);
122	enderton	1.1	if isequal(monchar,'JAN'), inimonth = 1;
123			elseif isequal(monchar,'FEB'), inimonth = 2;
124			elseif isequal(monchar,'MAR'), inimonth = 3;
125			elseif isequal(monchar,'APR'), inimonth = 4;
126			elseif isequal(monchar,'MAY'), inimonth = 5;
127			elseif isequal(monchar,'JUN'), inimonth = 6;
128			elseif isequal(monchar,'JUL'), inimonth = 7;
129			elseif isequal(monchar,'AUG'), inimonth = 8;
130			elseif isequal(monchar,'SEP'), inimonth = 9;
131			elseif isequal(monchar,'OCT'), inimonth = 10;
132			elseif isequal(monchar,'NOV'), inimonth = 11;
133	enderton	1.5	elseif isequal(monchar,'DEC'), inimonth = 12;
134			else, error('Can''t determine initial month in GRADS data!'); end
135	enderton	1.1	num = str2num(tokens{2});
136	enderton	1.5	% Determine num and inimonth parameters accounting for the unknown
137			% number of data files present when dealing with multiple data
138			% files. Assumption of monthly mean data is made.
139			if MultFiles;
140			initFileFound = 0;
141			for ifile = 1:length(datafile)
142			initFileFound = ~isempty(findstr(yrchar,datafile{ifile})) ...
143			\| initFileFound;
144			end
145			if ~initFileFound, inimonth = 1; end
146			num = 12.*length(datafile);
147			end
148			% Assumption of monthly mean data is made.
149	enderton	1.2	months=[inimonth:num+inimonth-1];
150			time=months/12; nt = length(months);
151	enderton	1.1	end
152
153	enderton	1.5	% Find where the desired variable is located.
154	enderton	1.1	if isequal(tokens{1},'VARS')
155			nv = str2num(tokens{2});
156			VARSline = iline;
157			end
158			if isequal(tokens{1},fln)
159			FIELDline = iline;
160	enderton	1.2	ivar = iline - VARSline;
161	enderton	1.1	end
162
163	enderton	1.5	% Find the values defonting undefined data.
164	enderton	1.1	if isequal(tokens{1},'UNDEF')
165			undef = str2num(tokens{2});
166			end
167
168	enderton	1.5	% Determine data format type.
169	enderton	1.2	if isequal(tokens{1},'FORMAT')
170	molod	1.4	if isequal(tokens{2},'SEQUENTIAL') \| isequal(tokens{2},'sequential')
171	enderton	1.2	format = 'SEQUENTIAL';
172			else
173			disp(['Unrecognized grads FORMAT: ',tokens{2}]);
174			end
175			end
176	enderton	1.1	end
177
178	enderton	1.2
179			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
180			% Verification and debugging %
181			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
182
183	enderton	1.5	% Verify that everything we need from the table file has been read in.
184	enderton	1.2	if DiagDebug, disp([' Debug -- grads data format: ',format]); end
185			parameters = {'xax' ,'yax' ,'zax' ,'time','nv' ,'undef','ivar'};
186			GRADSnames = {'XDEF','YDEF','ZDEF','TDEF','VARS','UNDEF',fln };
187	enderton	1.1	for ii = 1:length(parameters)
188			try
189			eval([parameters{ii},';']);
190			catch
191			error(['GRADS information not found: ',GRADSnames{ii}]);
192			end
193			end
194
195	enderton	1.2
196			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
197			% Read data file %
198			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
199
200	enderton	1.5	% Read in data from a single file.
201			if ~MultFiles
202			fid=fopen(datafile,'r','b');
203			data = fread(fid,'real*4');
204			fclose(fid);
205			if isequal(format,'SEQUENTIAL') \| isequal(format,'sequential')
206			index=true([nxnynznvnt+2nvnt,1]);
207			index([1:nxnynz+2:end])=false;
208			index([2:nxnynz+2:end])=false;
209			data = data(index);
210			end
211			data = reshape(data,[nx,ny,nz,nv,nt]);
212			data = data(:,:,:,ivar,:);
213			datatmp(:,:,:,:,inimonth:num+inimonth-1) = data;
214			datatmp(:,:,:,:,1:inimonth-1) = NaN;
215			data=datatmp;
216			% Read in data from multiple files.
217			else
218			for ifile = 1:length(datafile)
219			if DiagDebug, disp([' Debug -- Reading grads ',...
220			'data file: ',datafile{ifile}]); end
221			fid=fopen(datafile{ifile},'r','b');
222			datatemp = fread(fid,'real*4');
223			fclose(fid);
224			% Assumption of monthly mean data is made.
225			if MultFiles, nt = 12; end
226			if initFileFound & (inimonth ~= 1)
227			error(['Currently not able to handle Grads data spanning ',...
228			'multiple files not starting in January']);
229			end
230			datatemp = reshape(datatemp,[nx,ny,nz,nv,nt]);
231			datatemp = squeeze(datatemp(:,:,:,ivar,:));
232			data(1:nx,1:ny,1:nz,nt.(ifile-1)+1:nt.ifile) = datatemp;
233			end
234			end
235			data = squeeze(data);
236
237			% Turn undefined values into NaN.
238	enderton	1.1	tol = 1e-5;
239	molod	1.3	data( abs((data-undef)/undef) < tol ) = NaN;
240	enderton	1.5
241			% Shift data such that it is on a -180 to 180 longitude axis.
242	molod	1.3	if ShiftData
243			indexWestHemi = xax>=180;
244			indexEastHemi = xax<180;
245	enderton	1.5	data = cat(1,data(indexWestHemi,:,:,:),data(indexEastHemi,:,:,:));
246	molod	1.3	xax = cat(2,xax(indexWestHemi)-360,xax(indexEastHemi));
247	molod	1.4	end