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

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

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

revision 1.1 by enderton, Mon Jan 31 15:43:26 2005 UTC revision 1.5 by enderton, Mon Sep 5 18:48:27 2005 UTC
# Line 1  Line 1 
1  function [data,xax,yax,zax,months,time,dim] = ...  function [data,xax,yax,zax,months,time,dim] = ...
2      DiagLoadGradsData(Grads,dad,fln)      DiagLoadGradsData(Grads,dad,fln,DiagDebug)
3    
4  % Read table file.  
5  tablfile = [dad,'/',Grads,'.tabl'];  % Initial assumptions about data.
6  datafile = [dad,'/',Grads,'.data'];  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','');  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)  for iline = 1:length(file)
22      rem = file{iline}; tokens = {}; itoken = 0;      rem = file{iline}; tokens = {}; itoken = 0;
23      while isstr(rem)      while isstr(rem)
# Line 12  for iline = 1:length(file) Line 25  for iline = 1:length(file)
25          if ~isempty(token), tokens{itoken} = token; end          if ~isempty(token), tokens{itoken} = token; end
26      end      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')      if isequal(tokens{1},'XDEF')
47          if ~isequal(tokens{3},'LINEAR')          if ~isequal(tokens{3},'LINEAR')
48              error('Unable to deal with nonlinear grads x-axis data!');              error('Unable to deal with nonlinear grads x-axis data!');
# Line 20  for iline = 1:length(file) Line 51  for iline = 1:length(file)
51          ini = str2num(tokens{4});          ini = str2num(tokens{4});
52          inc = str2num(tokens{5});          inc = str2num(tokens{5});
53          xax = [ini:inc:ini+(num-1)*inc];          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          if min(xax) >= 0 && max(xax) > 180
57              xax = xax - 180;              ShiftData = 1;
58          end          end
59          nx = length(xax);          nx = length(xax);
60      end      end
61            
62        % Read in y-axis information.
63      if isequal(tokens{1},'YDEF')      if isequal(tokens{1},'YDEF')
64          if ~isequal(tokens{3},'LINEAR')          if ~isequal(tokens{3},'LINEAR')
65              error('Unable to deal with nonlinear grads y-axis data!');              error('Unable to deal with nonlinear grads y-axis data!');
# Line 36  for iline = 1:length(file) Line 70  for iline = 1:length(file)
70          yax = [ini:inc:ini+(num-1)*inc]; ny = length(yax);          yax = [ini:inc:ini+(num-1)*inc]; ny = length(yax);
71      end      end
72            
73        % Read in z-axis information.
74      if isequal(tokens{1},'ZDEF')      if isequal(tokens{1},'ZDEF')
75          if isequal(tokens{2},'1')          if isequal(tokens{2},'1')
76              zax = [0]; zax_found = 1; dim = 2;              zax = [0]; zax_found = 1; dim = 2;
77          else          else
78              if ~isequal(tokens{3},'LINEAR')              dim = 3;
79                  error('Unable to deal with nonlinear grads z-axis data!');              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              end
             num = str2num(tokens{2});  
             ini = str2num(tokens{4});  
             inc = str2num(tokens{5});  
             zax = [ini:inc:ini+(num-1)*inc]; dim = 3;  
92          end          end
93          nz = length(zax);          nz = length(zax);
94      end      end
95            
96        % Read in t-axis (time) information.
97      if isequal(tokens{1},'TDEF'), ini=tokens{4};      if isequal(tokens{1},'TDEF'), ini=tokens{4};
98          if ~isequal(tokens{3},'LINEAR')          if ~isequal(tokens{3},'LINEAR')
99              error('Currently unable to deal with nonlinear grads z-axis data!');              error('Currently unable to deal with nonlinear grads t-axis data!');
100          end          end
101          if ~isequal(tokens{5},'1mo')          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!');              error('Currently unable to deal with non monthly mean grads data!');
106          end          end
107          if length(ini) == 13          % Determine initial month and initial year.  This is done in an ad
108              monchar=ini(9:11);          % hoc manner in that depending on the length start time string, the
109          end          % 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;          if     isequal(monchar,'JAN'), inimonth = 1;
123          elseif isequal(monchar,'FEB'), inimonth = 2;          elseif isequal(monchar,'FEB'), inimonth = 2;
124          elseif isequal(monchar,'MAR'), inimonth = 3;          elseif isequal(monchar,'MAR'), inimonth = 3;
# Line 72  for iline = 1:length(file) Line 130  for iline = 1:length(file)
130          elseif isequal(monchar,'SEP'), inimonth = 9;          elseif isequal(monchar,'SEP'), inimonth = 9;
131          elseif isequal(monchar,'OCT'), inimonth = 10;          elseif isequal(monchar,'OCT'), inimonth = 10;
132          elseif isequal(monchar,'NOV'), inimonth = 11;          elseif isequal(monchar,'NOV'), inimonth = 11;
133          elseif isequal(monchar,'DEC'), inimonth = 12; end          elseif isequal(monchar,'DEC'), inimonth = 12;
134            else, error('Can''t determine initial month in GRADS data!'); end
135          num = str2num(tokens{2});          num = str2num(tokens{2});
136          months=[inimonth:num]; time=months/12; nt = length(months);          % 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      end
152            
153        % Find where the desired variable is located.
154      if isequal(tokens{1},'VARS')      if isequal(tokens{1},'VARS')
155          nv = str2num(tokens{2});          nv = str2num(tokens{2});
156          VARSline = iline;          VARSline = iline;
157      end      end
       
158      if isequal(tokens{1},fln)      if isequal(tokens{1},fln)
159          FIELDline = iline;          FIELDline = iline;
160          index = iline - VARSline;          ivar = iline - VARSline;
161      end      end
162            
163        % Find the values defonting undefined data.
164      if isequal(tokens{1},'UNDEF')      if isequal(tokens{1},'UNDEF')
165          undef = str2num(tokens{2});          undef = str2num(tokens{2});
166      end      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  end
177    
178  % Verify that things are there.  
179  parameters = {'xax' ,'yax' ,'zax' ,'time','nv'  ,'undef','index'};  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
180  GRADSnames = {'XDEF','YDEF','ZDEF','TDEF','VARS','UNDEF',fln    };  %                       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)  for ii = 1:length(parameters)
188      try      try
189          eval([parameters{ii},';']);          eval([parameters{ii},';']);
# Line 105  for ii = 1:length(parameters) Line 192  for ii = 1:length(parameters)
192      end      end
193  end  end
194    
195  % Read in grads data using axis infomation.  
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;  tol = 1e-5;
 fid=fopen(datafile,'r','b');  
 data = fread(fid,'real*4');  
 fclose(fid);  
 data = reshape(data,[nx,ny,nz,nv,nt]);  
 data = squeeze(data(:,:,:,index,:));  
 data( abs((data-undef)/undef) < tol ) = NaN;  
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
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.5
  ViewVC Help
Powered by ViewVC 1.1.22