% 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:
%
% Armstrong, B.N., Warner, J.C., Voulgaris, G., List, J.H., Thieler, E.R., 
% Martini, M.M., Montgomery, E.T., 2011, Carolinas Coastal Change Processes 
% Project Data Report for Observations near Diamond Shoals, NC from January-April 
% 2009: U.S. Geological Survey Open-File Report 2011-1156.
%
% 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."
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 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 = 'Z:\data\hatteras09\BBV_reviewed';
plotDirectory = 'D:\BNA\data\Carolinas\webpage\dataDVD\Hatteras2009DVD\mfiles';

% define data file names
dataFiles={'8531mc-a.nc','8532sc-a.nc','8533pt-a.nc',...
    '8541mc-a.nc',...
    '8543sc-a.nc','8554sc-a.nc',...
    '8562sc-a.nc','8563pt-a.nc',...
    '8571mc-a.nc','8573sc-a.nc',...
    '8582sc-a.nc'};
waveFiles={'8542p-cal.nc','8551p-cal.nc',...
    '8572p-cal.nc',''};

nFiles=length(dataFiles);

% define parameters for plots
startTime=julian(2009,1,8,00);
stopTime=julian(2009,05,15,00);
col=get(0,'defaultaxescolororder');     %plot color order
ylab={ {'Sigma-theta','(kg/m^3)'},...   %plot ylabels
    {'Salinity','(psu)'},...
    {'Conductivity','(S/m)'},...
    {'Temperature','(\circC)'},...
    {'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
    29 38;...
    3 5.5;...
    5  28;...
    0   4];
xt=[];
for yy=2009:2009
    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;

nc=netcdf(['Z:\data\hatteras09\candidate_final_files\',waveFiles{2},], '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(['Z:\data\hatteras09\candidate_final_files\',waveFiles{3},], 'nowrite');
jd_waves{2} = nc{'time'}(:) + (nc{'time2'}(:)/3600/1000/24);
waveDelta_t=gmean(diff(jd_waves{2}));
hsig{2}=nc{'wh_4061'}(:);
theFillValue=nc{'wh_4061'}.FillValue_(:);
badInd=find(hsig{2} >= theFillValue);
hsig{2}(badInd)=nan;
nc=close(nc);

hsig_north_comb = smart_interp([jd_waves{1}], [hsig{1}],...
    jd_int, 3);
hsig_south_comb = smart_interp([jd_waves{2}], [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
    nc=netcdf(fullfile(dataDirectory,dataFiles{f}), 'nowrite');
    jd_all=nc{'time'}(:) + (nc{'time2'}(:)/3600/1000/24);
    T_28_all=nc{'T_28'}(:);
    T_28_all(T_28_all>=1e35)=nan;
    C_51_all=nc{'C_51'}(:);
    C_51_all(C_51_all>=1e35)=nan;
    S_40_all=nc{'S_41'}(:);
    S_40_all(S_40_all>=1e35)=nan;
    STH_71_all=nc{'STH_71'}(:);
    STH_71_all(STH_71_all>=1e35)=nan;
    P_1_all=nc{'P_1'}(:);
    P_1_all(P_1_all>=1e35)=nan;
    site=nc.MOORING(:);
    inst_depth=nc.inst_height(:);
    inst_depth=nc.inst_depth(:);                   %instrument depth, m
    theFillValue=nc{'T_28'}.FillValue_(:);         %netcdf fill value
    delta_t = gmean(diff( jd_all));                  %data time step, s
    nc=close(nc);

    % 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;
    if f<=6
    ydata=[STH_71_all_int S_40_all_int C_51_all_int T_28_all_int hsig_north_comb];
    elseif f>6 && f<=10
    ydata=[STH_71_all_int S_40_all_int C_51_all_int T_28_all_int hsig_south_comb];
    else
    ydata=[STH_71_all_int S_40_all_int C_51_all_int T_28_all_int hsig_north_comb];
    end

    panelOrd=[1:min(size(ydata))];
    figure
    han=timeplt2(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');
        if h~=length(han)
            linehan=findobj(gca,'Type','line','color','b');
        end
    end

    % re-do xlabels
    axes(han(1))
    set(gca,'XTickLabel',{''})

    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),''])
    if inst_depth <= 1.0
        print(fullfile(plotDirectory,['ct_site',site,'_surface']),'-dpdf');
    else
        print(fullfile(plotDirectory,['ct_site',site,'_bottom']),'-dpdf');
        print(fullfile(plotDirectory,['ct_site',site,'_bottom']),'-dpng');
    end

    % get data ready for combined plot
    if f<=6
        STH_71_north_int_comb(:,f)=STH_71_all_int;
        S_40_north_int_comb(:,f)=S_40_all_int;
        C_51_north_int_comb(:,f)=C_51_all_int;
        T_28_north_int_comb(:,f)=T_28_all_int;
    elseif f>6 && f<=10
        STH_71_south_int_comb(:,f-6)=STH_71_all_int;
        S_40_south_int_comb(:,f-6)=S_40_all_int;
        C_51_south_int_comb(:,f-6)=C_51_all_int;
        T_28_south_int_comb(:,f-6)=T_28_all_int;
    else
        STH_71_tip_int_comb(:,f-10)=STH_71_all_int;
        S_40_tip_int_comb(:,f-10)=S_40_all_int;
        C_51_tip_int_comb(:,f-10)=C_51_all_int;
        T_28_tip_int_comb(:,f-10)=T_28_all_int;
    end

end
close all
%% now call timeplt.m to create a combined plot for the north site. 
xdata=jd_int;
ydata=[STH_71_north_int_comb S_40_north_int_comb,...
       C_51_north_int_comb T_28_north_int_comb hsig_north_comb];
panelOrd=[repmat(1,1,5),...
          repmat(2,1,5),...
          repmat(3,1,5),...
          repmat(4,1,5),...
          repmat(5,1,1)];
figure

han=timeplt(xdata,ydata,panelOrd);
clear global timeplt*
set(gcf,'ResizeFcn',[])
set(gcf,'Tag','')
dolandscape
fixpaper
LEG=legend(han(1),'Buoy microcat 1 mbs','Buoy seacat 4 mbs','Buoy PT 4 mbs',...
    'FloBee microcat 13 mbs',...
    'FloBee seacat 13 mbs','Minipod seacat 14 mbs');
% LEG=legend(han(1),'8531mc-a.nc 1 mbs','8532sc-a.nc 4 mbs',...
%     '8541mc-a.nc 13 mbs',...
%     '8543sc-a.nc 13 mbs','8554sc-a.nc 14 mbs');
 set(LEG,'Location','SouthOutside')   
 set(LEG,'Orientation','Horizontal')
figCol=get(gcf,'color');
set(LEG,'color',figCol)
set(LEG,'edgecolor',figCol)
pos=get(LEG,'position');
set(LEG,'position',[pos(1) pos(2)-0.07 pos(3) pos(4)-0.03])

% 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})
end

% re-do xlabels
axes(han(1))
oldhan=findobj(gca,'Tag','xtick labels');
delete(oldhan);
set(gca,'XTickLabel',{''})

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'],...
          ['North Site']})
print(fullfile(plotDirectory,['ct_north_sites']),'-dpdf');          

%% now call timeplt.m to create a combined plot for the south site. 
xdata=jd_int;
ydata=[STH_71_south_int_comb S_40_south_int_comb,...
       C_51_south_int_comb T_28_south_int_comb hsig_south_comb];
panelOrd=[repmat(1,1,3),...
          repmat(2,1,3),...
          repmat(3,1,3),...
          repmat(4,1,3),...
          repmat(5,1,1)];
figure

han=timeplt(xdata,ydata,panelOrd);
clear global timeplt*
set(gcf,'ResizeFcn',[])
set(gcf,'Tag','')
dolandscape
fixpaper
 LEG=legend(han(1),'Buoy seacat 4 mbs','Buoy PT 4 mbs',...
    'FloBee microcat 14 mbs','FloBee seacat 13 mbs');% set(LEG,'Location','EastOutside')   
%  LEG=legend(han(1),'8562sc-a.nc 4 mbs',...
%     '8571mc-a.nc 14 mbs','8573sc-a.nc 13 mbs');% set(LEG,'Location','EastOutside')   

 set(LEG,'Location','SouthOutside')   
 set(LEG,'Orientation','Horizontal')
figCol=get(gcf,'color');
set(LEG,'color',figCol)
set(LEG,'edgecolor',figCol)
pos=get(LEG,'position');
set(LEG,'position',[pos(1) pos(2)-0.07 pos(3) pos(4)-0.03])

% 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})
end

% re-do xlabels
axes(han(1))
oldhan=findobj(gca,'Tag','xtick labels');
delete(oldhan);
set(gca,'XTickLabel',{''})

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'],...
          ['South Site']})
print(fullfile(plotDirectory,['ct_south_sites']),'-dpdf');          

%% now call timeplt.m to create a combined plot for the tip site. 
xdata=jd_int;hsig_tip_comb=zeros(length(STH_71_tip_int_comb),1);
hsig_tip_comb(:)=nan;
ydata=[STH_71_tip_int_comb S_40_tip_int_comb,...
       C_51_tip_int_comb T_28_tip_int_comb hsig_tip_comb];
panelOrd=[repmat(1,1,1),...
          repmat(2,1,1),...
          repmat(3,1,1),...
          repmat(4,1,1),...
          repmat(5,1,1)];
figure

han=timeplt2(xdata,ydata,panelOrd);
clear global timeplt*
set(gcf,'ResizeFcn',[])
set(gcf,'Tag','')
dolandscape
fixpaper
 LEG=legend(han(1),'Micropod seacat 23 mbs');% set(LEG,'Location','EastOutside')   
%  LEG=legend(han(1),'8582sc-a.nc 23 mbs');% set(LEG,'Location','EastOutside')   
 set(LEG,'Location','SouthOutside')   
 set(LEG,'Orientation','Horizontal')
figCol=get(gcf,'color');
set(LEG,'color',figCol)
set(LEG,'edgecolor',figCol)
pos=get(LEG,'position');
set(LEG,'position',[pos(1) pos(2)-0.07 pos(3) pos(4)-0.03])

% 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})
end

% re-do xlabels
axes(han(1))
oldhan=findobj(gca,'Tag','xtick labels');
delete(oldhan);
set(gca,'XTickLabel',{''})

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'],...
          ['Tip Site']})
print(fullfile(plotDirectory,['ct_tip_sites']),'-dpdf');          
