/[MITgcm]/MITgcm_contrib/enderton/Diagnostics/DiagLoadGradsData.m
ViewVC logotype

Contents of /MITgcm_contrib/enderton/Diagnostics/DiagLoadGradsData.m

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.5 - (show annotations) (download)
Mon Sep 5 18:48:27 2005 UTC (19 years, 10 months ago) by enderton
Branch: MAIN
CVS Tags: HEAD
Changes since 1.4: +126 -41 lines
The long-awaited modification to DiagLoadGradsData.m to allow for reading in data spanning multiple files for Andrea.

1 function [data,xax,yax,zax,months,time,dim] = ...
2 DiagLoadGradsData(Grads,dad,fln,DiagDebug)
3
4
5 % Initial assumptions about data.
6 ShiftData = 0;
7 MultFiles = 0;
8 format = 'NONSEQUENTIAL';
9
10
11 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
12 % Parse table file %
13 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
14 tablfile = [dad,'/',Grads];
15 file = textread(tablfile,'%s','delimiter','\n','whitespace','');
16
17 % 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 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 % 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 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 % 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 if min(xax) >= 0 && max(xax) > 180
57 ShiftData = 1;
58 end
59 nx = length(xax);
60 end
61
62 % Read in y-axis information.
63 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 % Read in z-axis information.
74 if isequal(tokens{1},'ZDEF')
75 if isequal(tokens{2},'1')
76 zax = [0]; zax_found = 1; dim = 2;
77 else
78 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 end
92 end
93 nz = length(zax);
94 end
95
96 % Read in t-axis (time) information.
97 if isequal(tokens{1},'TDEF'), ini=tokens{4};
98 if ~isequal(tokens{3},'LINEAR')
99 error('Currently unable to deal with nonlinear grads t-axis data!');
100 end
101 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 error('Currently unable to deal with non monthly mean grads data!');
106 end
107 % 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 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 elseif isequal(monchar,'DEC'), inimonth = 12;
134 else, error('Can''t determine initial month in GRADS data!'); end
135 num = str2num(tokens{2});
136 % 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 months=[inimonth:num+inimonth-1];
150 time=months/12; nt = length(months);
151 end
152
153 % Find where the desired variable is located.
154 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 ivar = iline - VARSline;
161 end
162
163 % Find the values defonting undefined data.
164 if isequal(tokens{1},'UNDEF')
165 undef = str2num(tokens{2});
166 end
167
168 % Determine data format type.
169 if isequal(tokens{1},'FORMAT')
170 if isequal(tokens{2},'SEQUENTIAL') | isequal(tokens{2},'sequential')
171 format = 'SEQUENTIAL';
172 else
173 disp(['Unrecognized grads FORMAT: ',tokens{2}]);
174 end
175 end
176 end
177
178
179 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
180 % Verification and debugging %
181 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
182
183 % Verify that everything we need from the table file has been read in.
184 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 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
196 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
197 % Read data file %
198 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
199
200 % 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([nx*ny*nz*nv*nt+2*nv*nt,1]);
207 index([1:nx*ny*nz+2:end])=false;
208 index([2:nx*ny*nz+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 tol = 1e-5;
239 data( abs((data-undef)/undef) < tol ) = NaN;
240
241 % Shift data such that it is on a -180 to 180 longitude axis.
242 if ShiftData
243 indexWestHemi = xax>=180;
244 indexEastHemi = xax<180;
245 data = cat(1,data(indexWestHemi,:,:,:),data(indexEastHemi,:,:,:));
246 xax = cat(2,xax(indexWestHemi)-360,xax(indexEastHemi));
247 end

  ViewVC Help
Powered by ViewVC 1.1.22