Warning, /utils/matlab/Graphix/GraphixLoadGradsData.m is written in an unsupported language. File is not indexed.
view on githubraw file Latest commit c578fd6b on 2005-10-20 18:14:57 UTC
5bbd0d07c9 Dani*0001 function [data,xax,yax,zax,months,time,dim] = ...
0002 DiagLoadGradsData(Grads,dad,fln,DiagDebug)
0003
0004
0005 % Initial assumptions about data.
0006 ShiftData = 0;
0007 MultFiles = 0;
0008 format = 'NONSEQUENTIAL';
0009
0010
0011 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0012 % Parse table file %
0013 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0014 tablfile = [dad,'/',Grads];
0015 file = textread(tablfile,'%s','delimiter','\n','whitespace','');
0016
0017 % Cycle though lines of the GRADS table file, pasring the file line-by-line
0018 % to properly read in and format the data. The first line of every file
0019 % has an identifier denoting the information on that line. The following
0020 % blocks of code parse the information according to different identifiers.
0021 for iline = 1:length(file)
0022 rem = file{iline}; tokens = {}; itoken = 0;
c578fd6b1e Dani*0023 while ~isempty(rem)
5bbd0d07c9 Dani*0024 [token,rem] = strtok(rem); itoken = itoken + 1;
0025 if ~isempty(token), tokens{itoken} = token; end
0026 end
0027
0028 % Determine file(s) to read in.
0029 if isequal(tokens{1},'DSET')
0030 if isempty(findstr(tokens{2},'%'))
0031 datafile = [dad,tokens{2}(2:end)];
0032 else
0033 str = tokens{2}(2:end);
0034 alldatafile = strrep(str,'%y2','*');
0035 alldatafile = strrep(alldatafile,'%m2','*');
0036 files = ls([dad,alldatafile]);
0037 breaks = [0,find(isspace(files))];
0038 for ibreak = 1:length(breaks)-1
0039 datafile{ibreak} = files([breaks(ibreak)+1:breaks(ibreak+1)-1]);
0040 end
0041 MultFiles = 1;
0042 end
0043 end
0044
0045 % Read in x-axis information.
0046 if isequal(tokens{1},'XDEF')
0047 if ~isequal(tokens{3},'LINEAR')
0048 error('Unable to deal with nonlinear grads x-axis data!');
0049 end
0050 num = str2num(tokens{2});
0051 ini = str2num(tokens{4});
0052 inc = str2num(tokens{5});
0053 xax = [ini:inc:ini+(num-1)*inc];
0054 % If the x-axis is from 0 to 360, the axis information and the data
0055 % will later have to be shifted to a -180 to 180 axis.
0056 if min(xax) >= 0 && max(xax) > 180
0057 ShiftData = 1;
0058 end
0059 nx = length(xax);
0060 end
0061
0062 % Read in y-axis information.
0063 if isequal(tokens{1},'YDEF')
0064 if ~isequal(tokens{3},'LINEAR')
0065 error('Unable to deal with nonlinear grads y-axis data!');
0066 end
0067 num = str2num(tokens{2});
0068 ini = str2num(tokens{4});
0069 inc = str2num(tokens{5});
0070 yax = [ini:inc:ini+(num-1)*inc]; ny = length(yax);
0071 end
0072
0073 % Read in z-axis information.
0074 if isequal(tokens{1},'ZDEF')
0075 if isequal(tokens{2},'1')
0076 zax = [0]; zax_found = 1; dim = 2;
0077 else
0078 dim = 3;
0079 if isequal(tokens{3},'LINEAR')
0080 num = str2num(tokens{2});
0081 ini = str2num(tokens{4});
0082 inc = str2num(tokens{5});
0083 zax = [ini:inc:ini+(num-1)*inc];
0084 elseif isequal(tokens{3},'LEVELS')
0085 num = str2num(tokens{2});
0086 for inum = 1:num
0087 zax(inum) = 100.*str2num(tokens{inum+3});
0088 end
0089 else
0090 error('Unknown information in z-axis data!');
0091 end
0092 end
0093 nz = length(zax);
0094 end
0095
0096 % Read in t-axis (time) information.
0097 if isequal(tokens{1},'TDEF'), ini=tokens{4};
0098 if ~isequal(tokens{3},'LINEAR')
0099 error('Currently unable to deal with nonlinear grads t-axis data!');
0100 end
0101 if ~isequal(tokens{5},'1mo') & ~isequal(tokens{5},'01mo')
0102 % Note that monthly mean assumptions are also made when data
0103 % spanning multiple files are read in later on and that this
0104 % part must be updated if non-monthly mean data is allowed.
0105 error('Currently unable to deal with non monthly mean grads data!');
0106 end
0107 % Determine initial month and initial year. This is done in an ad
0108 % hoc manner in that depending on the length start time string, the
0109 % month is pluck out. There would be problems is there are start
0110 % time strings of different lengths have months in different
0111 % locations. In this case the month will not be found and an error
0112 % will be thrown. The year is only required for data spanning
0113 % multiple files, which thus fas has a start time string length of
0114 % 12. If it is not found, and error of undefined yrchar will be
0115 % cast and the code will need to be appropriately updated.
0116 if length(ini) == 13, monchar=ini(9:11);
0117 elseif length(ini) == 12, monchar=ini(6:8); yrchar = ini(9:12);
0118 elseif length(ini) == 10, monchar=ini(6:8);
0119 elseif length(ini) == 7, monchar=ini(1:3);
0120 else, error('Cannot parse TDEF correctly'); end
0121 monchar = upper(monchar);
0122 if isequal(monchar,'JAN'), inimonth = 1;
0123 elseif isequal(monchar,'FEB'), inimonth = 2;
0124 elseif isequal(monchar,'MAR'), inimonth = 3;
0125 elseif isequal(monchar,'APR'), inimonth = 4;
0126 elseif isequal(monchar,'MAY'), inimonth = 5;
0127 elseif isequal(monchar,'JUN'), inimonth = 6;
0128 elseif isequal(monchar,'JUL'), inimonth = 7;
0129 elseif isequal(monchar,'AUG'), inimonth = 8;
0130 elseif isequal(monchar,'SEP'), inimonth = 9;
0131 elseif isequal(monchar,'OCT'), inimonth = 10;
0132 elseif isequal(monchar,'NOV'), inimonth = 11;
0133 elseif isequal(monchar,'DEC'), inimonth = 12;
0134 else, error('Can''t determine initial month in GRADS data!'); end
0135 num = str2num(tokens{2});
0136 % Determine num and inimonth parameters accounting for the unknown
0137 % number of data files present when dealing with multiple data
0138 % files. Assumption of monthly mean data is made.
0139 if MultFiles;
0140 initFileFound = 0;
0141 for ifile = 1:length(datafile)
0142 initFileFound = ~isempty(findstr(yrchar,datafile{ifile})) ...
0143 | initFileFound;
0144 end
0145 if ~initFileFound, inimonth = 1; end
0146 num = 12.*length(datafile);
0147 end
0148 % Assumption of monthly mean data is made.
0149 months=[inimonth:num+inimonth-1];
0150 time=months/12; nt = length(months);
0151 end
0152
0153 % Find where the desired variable is located.
0154 if isequal(tokens{1},'VARS')
0155 nv = str2num(tokens{2});
0156 VARSline = iline;
0157 end
0158 if isequal(tokens{1},fln)
0159 FIELDline = iline;
0160 ivar = iline - VARSline;
0161 end
0162
0163 % Find the values defonting undefined data.
0164 if isequal(tokens{1},'UNDEF')
0165 undef = str2num(tokens{2});
0166 end
0167
0168 % Determine data format type.
0169 if isequal(tokens{1},'FORMAT')
0170 if isequal(tokens{2},'SEQUENTIAL') | isequal(tokens{2},'sequential')
0171 format = 'SEQUENTIAL';
0172 else
0173 disp(['Unrecognized grads FORMAT: ',tokens{2}]);
0174 end
0175 end
0176 end
0177
0178
0179 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0180 % Verification and debugging %
0181 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0182
0183 % Verify that everything we need from the table file has been read in.
0184 if DiagDebug, disp([' Debug -- grads data format: ',format]); end
0185 parameters = {'xax' ,'yax' ,'zax' ,'time','nv' ,'undef','ivar'};
0186 GRADSnames = {'XDEF','YDEF','ZDEF','TDEF','VARS','UNDEF',fln };
0187 for ii = 1:length(parameters)
0188 try
0189 eval([parameters{ii},';']);
0190 catch
0191 error(['GRADS information not found: ',GRADSnames{ii}]);
0192 end
0193 end
0194
0195
0196 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0197 % Read data file %
0198 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0199
0200 % Read in data from a single file.
0201 if ~MultFiles
0202 fid=fopen(datafile,'r','b');
0203 data = fread(fid,'real*4');
0204 fclose(fid);
0205 if isequal(format,'SEQUENTIAL') | isequal(format,'sequential')
0206 index=true([nx*ny*nz*nv*nt+2*nv*nt,1]);
0207 index([1:nx*ny*nz+2:end])=false;
0208 index([2:nx*ny*nz+2:end])=false;
0209 data = data(index);
0210 end
0211 data = reshape(data,[nx,ny,nz,nv,nt]);
0212 data = data(:,:,:,ivar,:);
0213 datatmp(:,:,:,:,inimonth:num+inimonth-1) = data;
0214 datatmp(:,:,:,:,1:inimonth-1) = NaN;
0215 data=datatmp;
0216 % Read in data from multiple files.
0217 else
0218 for ifile = 1:length(datafile)
0219 if DiagDebug, disp([' Debug -- Reading grads ',...
0220 'data file: ',datafile{ifile}]); end
0221 fid=fopen(datafile{ifile},'r','b');
0222 datatemp = fread(fid,'real*4');
0223 fclose(fid);
0224 % Assumption of monthly mean data is made.
0225 if MultFiles, nt = 12; end
0226 if initFileFound & (inimonth ~= 1)
0227 error(['Currently not able to handle Grads data spanning ',...
0228 'multiple files not starting in January']);
0229 end
0230 datatemp = reshape(datatemp,[nx,ny,nz,nv,nt]);
0231 datatemp = squeeze(datatemp(:,:,:,ivar,:));
0232 data(1:nx,1:ny,1:nz,nt.*(ifile-1)+1:nt.*ifile) = datatemp;
0233 end
0234 end
0235 data = squeeze(data);
0236
0237 % Turn undefined values into NaN.
0238 tol = 1e-5;
0239 data( abs((data-undef)/undef) < tol ) = NaN;
0240
0241 % Shift data such that it is on a -180 to 180 longitude axis.
0242 if ShiftData
0243 indexWestHemi = xax>=180;
0244 indexEastHemi = xax<180;
0245 data = cat(1,data(indexWestHemi,:,:,:),data(indexEastHemi,:,:,:));
0246 xax = cat(2,xax(indexWestHemi)-360,xax(indexEastHemi));
0247 end
