enderton/Diagnostics/DiagSlice.m

function [data,xax,yax,pltslc] = ...
    DiagSlice(data,fln,trl,dat,dad,grd,itr,tst,flu,ddf,gdf,...
              avg,slc,pst,LoadGridData,GridSuffix,ZcordFile);

% Function: DiagSlice
% Author:   Daniel Enderton
%
% Input Fields:
%
%   Field        Type        (Brief) Description
%   -----------------------------------------------------------------------
%   data         array       Field name
%   fln          string      Field name
%   dad          string      Data directory.
%   grd          string      Grid data directory.
%   avg          string      Averaging scheme ('Ann','DJF', 'JJA',...)
%   slc          string      Slice type ('Zon','k=#',...)
%   pst          string      Plot style ('Con','Sur',...)
%   flu          string      Fluid ('A','O')
%   dfm          string      Data format ('MDS' or 'MNC')
%   LoadGridData 0/1         Optionally load grid data
%
% Output Fields:
%
%   Field       Type        (Brief) Description
%   -----------------------------------------------------------------------
%   data        array       Sliced data.
%
% Descripton:
%   This function slices the data according to SliceType 'slc'.  Depending
%   on the SliceType and PlotStyle, the data may be converted to a lat-lon
%   grid.  At the end there is also a range check.  If the value is out of
%   the user specified range given in 'DaigFieldParamA'.

% Load diagnostics parameters, grid data, mark data size
DiagGenParam;
DiagFieldParamA;
DiagFieldParamO;
DiagFieldParamC;
DiagFieldParamI;
[XC,XG,YC,YG,Ylat,ZC,ZG,RAC,drC,drF,HFacC,HFacW,HFacS,dxG,dyG,dxC,dyC] = ...
    DiagLoadGridData(LoadGridData,grd,gdf,flu,GridSuffix,ZcordFile);
datasize = size(data);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                               Surface plot                              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% All surface plots are direct cube sphere outputs, and therefore must have
% the dimensions defined in 'hres' in 'DiagLoadGridData'.  This is commonly
% [192,32].  If the 'pst' is 'Grd' or 'Int', everything is already all set
% since those handle the cubed sphere data.  If the 'pst' is 'Con' or 'Cnf',
% convert data to a lat-lon grid so that it can easily be contours.
if isequal(slc,'Sur')
    if ismember(fln,{'TX','TY','USTR','VSTR'})
        data = data'; xax = XL; yax = YL;
    elseif ~isequal(datasize,[faces*hres,hres]) 
        error('Incorrect dimensions for slc:  ',slc);
    else
        if ismember(pst,{'Grd','Int'})
            if     isequal(pst,'Grd'),
                xax = XG(1:faces*hres,1:hres);
                yax = YG(1:faces*hres,1:hres);
            elseif isequal(pst,'Int'),
                xax = XC(1:faces*hres,1:hres);
                yax = YC(1:faces*hres,1:hres);
            end
        elseif ismember(pst,{'Con','Cnf'})
            data = cube2latlon(XC,YC,data,XL,YL)';
            xax = XL; yax = YL;
        else
            error(['slc = ''Sur'' can not use pst:  ',pst]);
        end
    end
    pltslc='lonlat';
    
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                            Zonal average plot                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% When computing a zonal average, there are many different options for the
% initial data.  If the data is cube-sphere, it could be of the size
% [hres,z-axis length] or just [hres].  In either case, 'calc_ZonAv_CS' (a
% nifty zonally averaging script from JMC) is called to compute the zonal
% average for raw cube sphere data.  For horizontal velocities, the data
% could be of the size [length(XL),length(YL),z-axis length], in which case
% you just need to take the mean over the longitude axis (always 1).
elseif isequal(slc,'Zon')
    if isequal(datasize(1:2),[faces*hres,hres])
        if isequal(flu,'O'), nBas = 0; end
            [data,dump1,dump2] = ...
                calc_ZonAv_CS(data,kpr,kwr,nBas,XC,YC,XG,YG,RAC,dad,HFacC);
        if isequal(avg,'Tse')
            data=data(:,:,1);  xax=years;  yax=Ylat; pltslc='timlat';
        elseif isequal(avg,'Tyr')
            data=data(:,:,1);  xax=[0:12]; yax=Ylat; pltslc='timlat';
        elseif isequal(pst,'Lin')
            data=data(:,:,1)'; xax=Ylat;   yax=NaN;  pltslc='latfld';
        else
            data=data(:,:,1)'; xax=Ylat;   yax=ZC;   pltslc='lathgt';
        end
    elseif isequal(datasize(1:2),[length(XL),length(YL)])
        if ~isequal(pst,'Lin')
            data = squeeze(mean(data,1))'; xax=YL; yax=ZC;  pltslc='lathgt';
        else
            data = squeeze(mean(data,1))'; xax=YL; yax=NaN; pltslc='latfld';
        end
    else
        error('Incorrect dimensions for slc = ''Zon''');
    end

    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                   i,j,k=#                               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% These slicing settings just take a slice on one of the axes.  If the
% slice option is 'i=' or 'j=' and the data is cube-sphere, it is converted
% to lat-lon first and indexed then (thus the index number should be based
% off of the XL and YL interpolated grids defined in 'DiagGenParam'.  It
% would be nice to add functions that can slice by longitude, latitude, and
% height values directly.  Note that there is transposing of all the data
% to prepare it for contour plots (harmless is line plots).
elseif isequal(slc(1:2),'i=')
    ii = str2num(slc(3:end));
    data = cube2latlon(XC,YC,data,XL,YL);
    if ismember(fln,fields2D)
        data=squeeze(data(ii,:));
        xax=YL; yax=NaN; pltslc='latfld';
    elseif ismember(fln,fields3D)
        data=squeeze(data(ii,:,:))';
        xax=YL; yax=ZC;  pltslc='lathgt';
    end
    
elseif isequal(slc(1:2),'j=')
    jj = str2num(slc(3:end));
    data = cube2latlon(XC,YC,data,XL,YL);
    if ismember(fln,fields2D)
        data = squeeze(data(:,jj));
        xax=XL; yax=NaN; pltslc='lonfld';
    elseif ismember(fln,fields3D)
        data = squeeze(data(:,jj,:))';
        xax=XL; yax=ZC;  pltslc='lonhgt';
    end
    
elseif isequal(slc(1:2),'k=')
    kk = str2num(slc(3:end));
    data = squeeze(data(:,:,kk));
    if ismember(fln,{'U','V','uVel','vVel','fizhi_U','fizhi_V'})
        data = data'; xax = XL; yax = YL;
    elseif ismember(pst,{'Grd','Int'})
        if isequal(pst,'Grd'),
            xax = XG(1:faces*hres,1:hres);
            yax = YG(1:faces*hres,1:hres);
        elseif isequal(pst,'Int'),
            xax = XC(1:faces*hres,1:hres);
            yax = YC(1:faces*hres,1:hres);
        end
    elseif ismember(pst,{'Con','Cnf'})
        data = cube2latlon(XC,YC,data,XL,YL)'; xax = XL; yax = YL;
    end
    pltslc='lonlat';
    
else
    error(['Unrecognized SliceType:  ',slc]);
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                               Check range                               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Load fixed and data ranges, throw a warning is data out of range.
datarange = [min(data(:)),max(data(:))];
try
    eval(['fixedrange = ',fln,'range',flu,';']);
    if datarange(1) < fixedrange(1) | ...
       datarange(2) > fixedrange(2)
        disp(['***Warning***  Value out of range for ',fln]);
        disp(['               Data range:  ',mat2str(datarange)]);
        disp(['               Fixed range:  ',mat2str(fixedrange)]);
    end
catch
    disp(['***Warning***  No range information found for ',fln]);
    disp(['               Data range:  ',mat2str(datarange)]);
end
1	enderton	1.1	function [data,xax,yax,pltslc] = ...
2	enderton	1.2	DiagSlice(data,fln,trl,dat,dad,grd,itr,tst,flu,ddf,gdf,...
3			avg,slc,pst,LoadGridData,GridSuffix,ZcordFile);
4	enderton	1.1
5			% Function: DiagSlice
6			% Author: Daniel Enderton
7			%
8			% Input Fields:
9			%
10			% Field Type (Brief) Description
11			% -----------------------------------------------------------------------
12			% data array Field name
13			% fln string Field name
14			% dad string Data directory.
15			% grd string Grid data directory.
16			% avg string Averaging scheme ('Ann','DJF', 'JJA',...)
17			% slc string Slice type ('Zon','k=#',...)
18			% pst string Plot style ('Con','Sur',...)
19			% flu string Fluid ('A','O')
20			% dfm string Data format ('MDS' or 'MNC')
21			% LoadGridData 0/1 Optionally load grid data
22			%
23			% Output Fields:
24			%
25			% Field Type (Brief) Description
26			% -----------------------------------------------------------------------
27			% data array Sliced data.
28			%
29			% Descripton:
30			% This function slices the data according to SliceType 'slc'. Depending
31			% on the SliceType and PlotStyle, the data may be converted to a lat-lon
32			% grid. At the end there is also a range check. If the value is out of
33			% the user specified range given in 'DaigFieldParamA'.
34
35			% Load diagnostics parameters, grid data, mark data size
36			DiagGenParam;
37			DiagFieldParamA;
38			DiagFieldParamO;
39			DiagFieldParamC;
40			DiagFieldParamI;
41			[XC,XG,YC,YG,Ylat,ZC,ZG,RAC,drC,drF,HFacC,HFacW,HFacS,dxG,dyG,dxC,dyC] = ...
42			DiagLoadGridData(LoadGridData,grd,gdf,flu,GridSuffix,ZcordFile);
43			datasize = size(data);
44
45
46			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
47			% Surface plot %
48			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
49
50			% All surface plots are direct cube sphere outputs, and therefore must have
51			% the dimensions defined in 'hres' in 'DiagLoadGridData'. This is commonly
52			% [192,32]. If the 'pst' is 'Grd' or 'Int', everything is already all set
53			% since those handle the cubed sphere data. If the 'pst' is 'Con' or 'Cnf',
54			% convert data to a lat-lon grid so that it can easily be contours.
55			if isequal(slc,'Sur')
56			if ismember(fln,{'TX','TY','USTR','VSTR'})
57			data = data'; xax = XL; yax = YL;
58			elseif ~isequal(datasize,[faces*hres,hres])
59			error('Incorrect dimensions for slc: ',slc);
60			else
61			if ismember(pst,{'Grd','Int'})
62			if isequal(pst,'Grd'),
63			xax = XG(1:faces*hres,1:hres);
64			yax = YG(1:faces*hres,1:hres);
65			elseif isequal(pst,'Int'),
66			xax = XC(1:faces*hres,1:hres);
67			yax = YC(1:faces*hres,1:hres);
68			end
69			elseif ismember(pst,{'Con','Cnf'})
70			data = cube2latlon(XC,YC,data,XL,YL)';
71			xax = XL; yax = YL;
72			else
73			error(['slc = ''Sur'' can not use pst: ',pst]);
74			end
75			end
76			pltslc='lonlat';
77
78
79			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
80			% Zonal average plot %
81			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
82
83			% When computing a zonal average, there are many different options for the
84			% initial data. If the data is cube-sphere, it could be of the size
85			% [hres,z-axis length] or just [hres]. In either case, 'calc_ZonAv_CS' (a
86			% nifty zonally averaging script from JMC) is called to compute the zonal
87			% average for raw cube sphere data. For horizontal velocities, the data
88			% could be of the size [length(XL),length(YL),z-axis length], in which case
89			% you just need to take the mean over the longitude axis (always 1).
90			elseif isequal(slc,'Zon')
91			if isequal(datasize(1:2),[faces*hres,hres])
92			if isequal(flu,'O'), nBas = 0; end
93			[data,dump1,dump2] = ...
94			calc_ZonAv_CS(data,kpr,kwr,nBas,XC,YC,XG,YG,RAC,dad,HFacC);
95			if isequal(avg,'Tse')
96			data=data(:,:,1); xax=years; yax=Ylat; pltslc='timlat';
97			elseif isequal(avg,'Tyr')
98			data=data(:,:,1); xax=[0:12]; yax=Ylat; pltslc='timlat';
99			elseif isequal(pst,'Lin')
100			data=data(:,:,1)'; xax=Ylat; yax=NaN; pltslc='latfld';
101			else
102			data=data(:,:,1)'; xax=Ylat; yax=ZC; pltslc='lathgt';
103			end
104			elseif isequal(datasize(1:2),[length(XL),length(YL)])
105	molod	1.4	if ~isequal(pst,'Lin')
106	enderton	1.1	data = squeeze(mean(data,1))'; xax=YL; yax=ZC; pltslc='lathgt';
107	molod	1.4	else
108	enderton	1.1	data = squeeze(mean(data,1))'; xax=YL; yax=NaN; pltslc='latfld';
109			end
110			else
111			error('Incorrect dimensions for slc = ''Zon''');
112			end
113
114
115			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
116			% i,j,k=# %
117			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
118
119			% These slicing settings just take a slice on one of the axes. If the
120			% slice option is 'i=' or 'j=' and the data is cube-sphere, it is converted
121			% to lat-lon first and indexed then (thus the index number should be based
122			% off of the XL and YL interpolated grids defined in 'DiagGenParam'. It
123			% would be nice to add functions that can slice by longitude, latitude, and
124			% height values directly. Note that there is transposing of all the data
125			% to prepare it for contour plots (harmless is line plots).
126			elseif isequal(slc(1:2),'i=')
127			ii = str2num(slc(3:end));
128			data = cube2latlon(XC,YC,data,XL,YL);
129			if ismember(fln,fields2D)
130			data=squeeze(data(ii,:));
131			xax=YL; yax=NaN; pltslc='latfld';
132			elseif ismember(fln,fields3D)
133			data=squeeze(data(ii,:,:))';
134			xax=YL; yax=ZC; pltslc='lathgt';
135			end
136
137			elseif isequal(slc(1:2),'j=')
138			jj = str2num(slc(3:end));
139			data = cube2latlon(XC,YC,data,XL,YL);
140			if ismember(fln,fields2D)
141			data = squeeze(data(:,jj));
142			xax=XL; yax=NaN; pltslc='lonfld';
143			elseif ismember(fln,fields3D)
144			data = squeeze(data(:,jj,:))';
145			xax=XL; yax=ZC; pltslc='lonhgt';
146			end
147
148			elseif isequal(slc(1:2),'k=')
149			kk = str2num(slc(3:end));
150			data = squeeze(data(:,:,kk));
151	molod	1.3	if ismember(fln,{'U','V','uVel','vVel','fizhi_U','fizhi_V'})
152	enderton	1.1	data = data'; xax = XL; yax = YL;
153			elseif ismember(pst,{'Grd','Int'})
154			if isequal(pst,'Grd'),
155			xax = XG(1:faces*hres,1:hres);
156			yax = YG(1:faces*hres,1:hres);
157			elseif isequal(pst,'Int'),
158			xax = XC(1:faces*hres,1:hres);
159			yax = YC(1:faces*hres,1:hres);
160			end
161			elseif ismember(pst,{'Con','Cnf'})
162			data = cube2latlon(XC,YC,data,XL,YL)'; xax = XL; yax = YL;
163			end
164			pltslc='lonlat';
165
166			else
167			error(['Unrecognized SliceType: ',slc]);
168			end
169
170
171			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
172			% Check range %
173			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
174
175			% Load fixed and data ranges, throw a warning is data out of range.
176			datarange = [min(data(:)),max(data(:))];
177			try
178			eval(['fixedrange = ',fln,'range',flu,';']);
179			if datarange(1) < fixedrange(1) \| ...
180			datarange(2) > fixedrange(2)
181			disp(['*Warning* Value out of range for ',fln]);
182			disp([' Data range: ',mat2str(datarange)]);
183			disp([' Fixed range: ',mat2str(fixedrange)]);
184			end
185			catch
186			disp(['*Warning* No range information found for ',fln]);
187			disp([' Data range: ',mat2str(datarange)]);
188	molod	1.3	end