--- MITgcm_contrib/enderton/Diagnostics/DiagLoadGradsData.m 2005/01/31 15:43:26 1.1 +++ MITgcm_contrib/enderton/Diagnostics/DiagLoadGradsData.m 2005/09/05 18:48:27 1.5 @@ -1,10 +1,23 @@ function [data,xax,yax,zax,months,time,dim] = ... - DiagLoadGradsData(Grads,dad,fln) + DiagLoadGradsData(Grads,dad,fln,DiagDebug) -% Read table file. -tablfile = [dad,'/',Grads,'.tabl']; -datafile = [dad,'/',Grads,'.data']; + +% 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) @@ -12,6 +25,24 @@ 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!'); @@ -20,12 +51,15 @@ 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 - xax = 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!'); @@ -36,31 +70,55 @@ 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 - if ~isequal(tokens{3},'LINEAR') - error('Unable to deal with nonlinear grads z-axis data!'); + 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 - num = str2num(tokens{2}); - ini = str2num(tokens{4}); - inc = str2num(tokens{5}); - zax = [ini:inc:ini+(num-1)*inc]; dim = 3; 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 z-axis data!'); + error('Currently unable to deal with nonlinear grads t-axis data!'); end - if ~isequal(tokens{5},'1mo') + 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 - if length(ini) == 13 - monchar=ini(9:11); - 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; @@ -72,31 +130,60 @@ elseif isequal(monchar,'SEP'), inimonth = 9; elseif isequal(monchar,'OCT'), inimonth = 10; elseif isequal(monchar,'NOV'), inimonth = 11; - elseif isequal(monchar,'DEC'), inimonth = 12; end + elseif isequal(monchar,'DEC'), inimonth = 12; + else, error('Can''t determine initial month in GRADS data!'); end num = str2num(tokens{2}); - months=[inimonth:num]; time=months/12; nt = length(months); - + % 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; - index = iline - VARSline; + 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 -% Verify that things are there. -parameters = {'xax' ,'yax' ,'zax' ,'time','nv' ,'undef','index'}; -GRADSnames = {'XDEF','YDEF','ZDEF','TDEF','VARS','UNDEF',fln }; + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% 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},';']); @@ -105,11 +192,56 @@ end end -% Read in grads data using axis infomation. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% 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; -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; \ No newline at end of file +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