% ct_ts.m  An mfile to plot the timeseries of temperature, salinity,
%          conductivity, and sigma-theta data from South Carolina CT's.  

%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Use of this program is described in:
%
% Sullivan, C.M., Warner, J.C., Martini, M.A., Voulgaris, G., 
% Work, P.A., Haas, K.A., and Hanes, D.H. (2006) 
% South Carolina Coastal Erosion Study Data Report for Observations
% October 2003 - April 2004., USGS Open-File Report 2005-1429.
%
% Program written in Matlab v7.1.0 SP3
% Program ran on PC with Windows XP Professional OS.
%
% "Although this program has been used by the USGS, no warranty, 
% expressed or implied, is made by the USGS or the United States 
% Government as to the accuracy and functioning of the program 
% and related program material nor shall the fact of distribution 
% constitute any such warranty, and no responsibility is assumed 
% by the USGS in connection therewith."
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% C. Sullivan 09/21/04
%             02/09/05, Interpolate data onto a standard time axis.  Name
%                       pdf's ct_site#_surface(bottom).pdf.
%             02/24/05, Reformatting of labels and ticks
%             05/24/06, Update paths to files and file names.  Update to
%                       show Hs from site 2 in the top panel of the plot. 
%                       
% Dependencies:
%       the netcdf toolbox
%       julian.m
%       timeplt.m
%       dolandscape.m
%       fixpaper.m
%       suptitle.m

% TODO:  -ask Polly about site 6 metadata b/c it is wrong
%        -cleaned nc files *-b.nc need metadata copied from *-a.nc

clear
more off

%define data and output directories
dataDirectory='E:\CSULLIVAN\SOUTH_CAROLINA\DATA_DVD\DATAFILES';
plotDirectory='E:\CSULLIVAN\SOUTH_CAROLINA\DATA_DVD\PLOTS\CT';

% define data file names
dataFiles={'7231mc-a.nc','7431mc-a.nc',...
           '7222sc-a.nc','7422sc-a.nc',...
           '7251mc-a.nc','7451mc-a.nc',...
           '7242sc-a.nc','7442sc-a.nc',...
           '7265sc-a.nc',['7464Asc-a.nc';'7464Bsc-a.nc'],...
           '7291mc-a.nc','7491mc-a.nc'};
nFiles=length(dataFiles);

% define parameters for plots
startTime=julian(2003,10,25,00);        %plot start stime
stopTime=julian(2004,04,25,00);         %plot stop time  
col=get(0,'defaultaxescolororder');     %plot color order
ylab={ {'Sigma-theta','(kg/m^3)'},...   %plot ylabels
       {'Salinity','(psu)'},...
       {'Conductivity','(S/m)'},...
       {'Temperature','(\circC)'},...
       {'Site 2','Significant Wave','Height (m)'} };
timediff=stopTime-startTime;                                %plot xlimits extend 5 percent
xlims=[(startTime-timediff*0.05) (stopTime+timediff*0.05)]; %past start time and stop time
ylims=[ 22 28;...   %plot ylimits
        32 36;...
        3 5.5;...
        5  28;...
        0   2];
xt=[];
for yy=2003:2004
    for mm=1:12
        tic=julian(yy,mm,01,00); 
        if tic>=xlims(1) & tic<=xlims(2)
            xt=[xt; tic];   % place xticks on the first of every month
        end
    end
end

% 'jd_int' is the standard time array onto which
% data is interpolated. data gaps larger than
% 'maxGap' hours are not interpolated over
jd_int = [startTime:1/24:stopTime]';
maxGap = 1;

% load site 2 wave data for the top panel
% of the plot
nc=netcdf(fullfile(dataDirectory,'7221adwp-cal.nc'), 'nowrite');
jd_waves{1} = nc{'time'}(:) + (nc{'time2'}(:)/3600/1000/24);
waveDelta_t=gmean(diff(jd_waves{1}));
hsig{1}=nc{'wh_4061'}(:);
theFillValue=nc{'wh_4061'}.FillValue_(:);
badInd=find(hsig{1} >= theFillValue);
hsig{1}(badInd)=nan;
nc=close(nc);

nc=netcdf(fullfile(dataDirectory,'7421adwp-cal.nc'), 'nowrite');
jd_waves{2} = nc{'time'}(:) + (nc{'time2'}(:)/3600/1000/24);
hsig{2}=nc{'wh_4061'}(:);
theFillValue=nc{'wh_4061'}.FillValue_(:);
badInd=find(hsig{2} >= theFillValue);
hsig{2}(badInd)=nan;
nc=close(nc);

hsig_all_comb = smart_interp([jd_waves{1};jd_waves{2}], [hsig{1};hsig{2}],...
                              jd_int, 3);
                     
% initialize variables
jd=cell(nFiles,1);
T_28=cell(nFiles,1);
C_51=cell(nFiles,1);
S_40=cell(nFiles,1);
STH_71=cell(nFiles,1);



for f=1:length(dataFiles)
     
     % load data
     if min(size(dataFiles{f})) > 1
         % combine files for 2nd deployment at site 6
         % where there are two files
         nc1=netcdf(dataFiles{f}(1,:), 'nowrite');
         nc2=netcdf(dataFiles{f}(2,:), 'nowrite');
         jd{f}=[nc1{'time'}(:) + (nc1{'time2'}(:)/3600/1000/24);...
                nc2{'time'}(:) + (nc2{'time2'}(:)/3600/1000/24)];
         T_28{f}=[nc1{'T_28'}(:); nc2{'T_28'}(:)];       %Temperature, deg. C
         C_51{f}=[nc1{'C_51'}(:); nc2{'C_51'}(:)];       %Conductivity, S/m
         S_40{f}=[nc1{'S_40'}(:); nc2{'S_40'}(:)];       %Salinity, psu
         STH_71{f}=[nc1{'STH_71'}(:); nc2{'STH_71'}(:)]; %Sigma-theta, kg/m^3
         site='6';
         inst_depth=6.7; %this is wrong!
         nc1=close(nc1);
         nc2=close(nc2);
     else
         nc=netcdf(dataFiles{f}, 'nowrite');
         jd{f}=nc{'time'}(:) + (nc{'time2'}(:)/3600/1000/24);
         T_28{f}=nc{'T_28'}(:);
         C_51{f}=nc{'C_51'}(:);
         S_40{f}=nc{'S_40'}(:);
         STH_71{f}=nc{'STH_71'}(:);
         site=nc.DATA_CMNT(6);                          %site number
         inst_depth=nc.inst_depth(:);                   %instrument depth, m
         theFillValue=nc{'T_28'}.FillValue_(:);         %netcdf fill value
         delta_t = gmean(diff(jd{f}));                  %data time step, s
         nc=close(nc);
     end

     if isequal(f,8)
         nc=netcdf('E:\CSULLIVAN\SOUTH_CAROLINA\DATA_DVD\DATAFILES\7442sc-b.nc', 'nowrite');
         site='3';
         inst_depth=9.6;
         nc=close(nc);
     end
     
     if mod(f,2)==0
         
            % combine data from deployment 1 and
            % deployment 2, and replace netcdf
            % fill values with nan's.
            jd_all=[jd{f-1}; jd{f}];
            T_28_all=[T_28{f-1}; T_28{f}];
            T_28_all(T_28_all>=1e35)=nan;
            C_51_all=[C_51{f-1}; C_51{f}];
            C_51_all(C_51_all>=1e35)=nan;
            S_40_all=[S_40{f-1}; S_40{f}];
            S_40_all(S_40_all>=1e35)=nan;
            STH_71_all=[STH_71{f-1}; STH_71{f}];
            STH_71_all(STH_71_all>=1e35)=nan;
            
            % interpolate data onto 'jd_int'
            T_28_all_int = smart_interp(jd_all, T_28_all, jd_int, maxGap);
            C_51_all_int = smart_interp(jd_all, C_51_all, jd_int, maxGap);
            S_40_all_int = smart_interp(jd_all, S_40_all, jd_int, maxGap);
            STH_71_all_int = smart_interp(jd_all, STH_71_all, jd_int, maxGap);
            
            % call timeplt.m to create individual plots
            % (clearing global timeplt*, the ResizeFcn,
            % and the Tag are necessary to release
            % timeplt.m's control over your figure!)            xdata=jd_int;
            xdata=jd_int;
            ydata=[STH_71_all_int S_40_all_int C_51_all_int T_28_all_int hsig_all_comb];
            panelOrd=[1:min(size(ydata))];            
            figure
            han=timeplt(xdata,ydata,panelOrd);
            clear global timeplt*
            set(gcf,'ResizeFcn',[])
            set(gcf,'Tag','')
            dolandscape
            fixpaper

            % reformat xlims, ylims, ylabels, and line color
            for h=1:length(han)
                axes(han(h))
                set(gca,'xlim',xlims,'xtick',xt)
                set(gca,'ylim',ylims(h,:))
                ylabhan=get(gca,'Ylabel');
                set(ylabhan,'String',ylab{h})
                pos=get(ylabhan,'position');
                set(get(gca,'Ylabel'),'position',[julian(2003,10,8) pos(2) pos(3)]) %align ylabels nicely
                if h~=length(han)
                    linehan=findobj(gca,'Type','line','color','b');
                    set(linehan,'color',col(f/2,:));
                end
            end

            % re-do xlabels
            axes(han(1))
            oldhan=findobj(gca,'Tag','xtick labels');
            delete(oldhan);
            for tic=1:length(xt)
                xtg=gregorian(xt(tic));
                xtd=datenum(xtg);
                xtl={[datestr(xtd,3),' 1'];datestr(xtd,10)}; %2 line label
                text(xt(tic),21.7,xtl,'horizontalalignment','center',...
                    'verticalalignment','top')
            end

            % give it a title and export to .pdf
            set(gcf,'PaperPositionMode','Auto')
            suptitle({'Salinity and Temperature Observations',['Site ',num2str(site),],...
                [num2str(inst_depth),' mbs']})
            if inst_depth <= 1.0
                print(fullfile(plotDirectory,['ct_site',site,'_surface']),'-dpdf');
            else
                print(fullfile(plotDirectory,['ct_site',site,'_bottom']),'-dpdf');
            end

            % get data ready for combined plot
            STH_71_all_int_comb(:,f/2)=STH_71_all_int;
            S_40_all_int_comb(:,f/2)=S_40_all_int;
            C_51_all_int_comb(:,f/2)=C_51_all_int;
            T_28_all_int_comb(:,f/2)=T_28_all_int;

     end

end
          
% now call timeplt.m to create a combined plot. 
xdata=jd_int;
ydata=[STH_71_all_int_comb S_40_all_int_comb,...
       C_51_all_int_comb T_28_all_int_comb hsig_all_comb];
panelOrd=[repmat(1,1,nFiles/2),...
          repmat(2,1,nFiles/2),...
          repmat(3,1,nFiles/2),...
          repmat(4,1,nFiles/2),...
          repmat(5,1,1)];
figure
han=timeplt(xdata,ydata,panelOrd);
clear global timeplt*
set(gcf,'ResizeFcn',[])
set(gcf,'Tag','')
LEG=legend(han(4),'Site 2 (1 mbs)','Site 2 (8.1 mbs)',...
                  'Site 3 (1 mbs)','Site 3 (9.6 mbs)',...
                  'Site 6 (6.7 mbs)',...
                  'Site 7 (1 mbs)');
set(LEG,'Location','NorthOutside')   
set(LEG,'Orientation','Horizontal')
figCol=get(gcf,'color');
set(LEG,'color',figCol)
set(LEG,'edgecolor',figCol)
dolandscape
fixpaper
pos=get(LEG,'position');
set(LEG,'position',[pos(1) .95 pos(3) pos(4)])

% reformat xlims, ylims, ylabels
for h=1:length(han)
    axes(han(h))
    set(gca,'xlim',xlims,'xtick',xt)
    set(gca,'ylim',ylims(h,:))
    ylabhan=get(gca,'Ylabel');
    set(ylabhan,'String',ylab{h})
    pos=get(ylabhan,'position');
    set(get(gca,'Ylabel'),'position',[julian(2003,10,8) pos(2) pos(3)]) %align ylabels nicely
end

% re-do xlabels
axes(han(1))
oldhan=findobj(gca,'Tag','xtick labels');
delete(oldhan);
for tic=1:length(xt)
    xtg=gregorian(xt(tic));
    xtd=datenum(xtg);
    xtl={[datestr(xtd,3),' 1'];datestr(xtd,10)}; %2 line label
    text(xt(tic),21.7,xtl,'horizontalalignment','center',...
        'verticalalignment','top')
end

% give it a title and export to .pdf
set(gcf,'PaperPositionMode','Auto')
suptitle({['Salinity and Temperature Observations'],...
          ['All Sites']})
print(fullfile(plotDirectory,['ct_all_sites']),'-dpdf');          
