% 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, Brandy N., Warner, John C., List, Jeff H., Martini, 
% Marinna A., Montgomery, Ellyn T., Traykovski, P., Traykovski, Peter, Voulgaris, George 
% 2014, Coastal Change Processes Project data report for oceanographic 
% observations near Fire Island, New York, February through May  
% 2014 : U.S. Geological Survey Open-File Report 2014-xxxx, available 
% online at http://pubs.usgs.gov/of/2014/xxxx/.
%
% Program written in Matlab v8.0.0.755
% Program ran on PC with Windows 7 64bit

%
% "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:
%       julian.m
%       gregorian.m
%       smart_interp.m
%       PL66TN.m
%       suptitle.m
%       tight_subplot
%       draw.m
%       gmean.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
filt=false;
%define data and output directories
dataDirectory = 'c:\home\data\validation\FI14\candidate_final_files';
plotDirectory = 'C:\Users\emontgomery\doc\OFR_FI14\FI_14_WEBOFR\plotsem';

% define data file names
dataFiles={'9852mc-a.cdf','9871mc-a.cdf','9901mc-a.cdf',...
    '9932mc-a.cdf','9951mc-a.cdf','9971mc-a.cdf',...
    '9991mc-a.cdf','10012mc-a.cdf',...
    '9862mc-a.cdf','9881sc-a.cdf','9912sc-a.cdf',...
    '9982sc-a.cdf','10002sc-a.cdf','10023mc-a.cdf','10042sc-a.cdf'};

names={'Site 1s','Site 2s','Site 3s','Site 4s','Site 5s','Siet 6s','Site 7s','Site 8s',...
    'Site 1m','Site 2f','Site 3f','Site 6m','Site 7m','Site 8m','Site 9'};

waveFiles={'9861whp-calBS-trm.nc','9891whp-cal404-trm.nc','9921p-cal404-trm.nc',...
    '9981whp-cal404-trm.nc','10022sgWvs-p_trim-d.nc'};
wavenameFiles={'Site 1','Site 2','Site 3','Site 6','Site 8'};


%%
% waveFiles={'9861whp-cal.nc','9892awWvs-p.nc',...
%     '9921whp-cal.nc','9981whp-cal.nc','10001whp-cal.nc','10021whp-cal.nc','10041whp-cal.nc'};
% wavenameFiles={'Site 1','Site 2','Site 3','Site 6','Site 7','Site 8','Site 9'};

nFiles=length(dataFiles);
ncVars={'T_28','C_51','S_41'};
vars={'C','Sm','psu'};
% define parameters for plots

startTime=julian(2014,2,8,00);                   
stopTime=julian(2014,05,5,00);         
col=colormap;     %plot color order
ylab={ {'Significant Wave','Height (m)'}...   %plot ylabels
    {'Temperature','(\circC)'},{'Conductivity','(S/m)'},{'Salinity','(psu)'}};
timediff=stopTime-startTime;                                %plot xlimits extend 5 percent
xlims=[startTime-timediff*0.025 stopTime+timediff*0.025]; %past start time and stop time
ylims=[ 0   5;...
   0 15;...   %plot ylimits
    2.6 4;...
    30 35];
xt=[];
for yy=2014:2014
     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
delta_t = 900/(3600*24); %for 15 min. data
jd_int = [startTime:delta_t:stopTime]';  % this is the datenum equiv of jd_press,
maxGap = 1;
%
for f=1:length(waveFiles)
    ncwaves=(fullfile(dataDirectory,waveFiles{f}));
    time=double(ncread(ncwaves,'time'));
    time2=double(ncread(ncwaves,'time2'));
    jd_waves = time + (time2/3600/1000/24);
    waveDelta_t=gmean(diff(jd_waves));
    hsig=squeeze(ncread(ncwaves,'wh_4061'));
    try
        theFillValue=ncreadatt(ncwaves,'wh_4061','_FillValue');
    catch
        theFillValue=1.0E35;
    end
    hsig(hsig >= theFillValue)=nan;
    hsig(hsig <= 0)=nan;
    hsig_int(:,f) = smart_interp([jd_waves], [hsig],...
        jd_int,3);
end

% initialize variables
jd=zeros(length(jd_int),nFiles);
T_28=zeros(length(jd_int),nFiles);
C_51=zeros(length(jd_int),nFiles);
S_41=zeros(length(jd_int),nFiles);

%%
for f=1:length(dataFiles)
    clear ncobj;
    % Load the data
    nc=(fullfile(dataDirectory,dataFiles{f}));
    time=double(ncread(nc,'time'));
    time2=double(ncread(nc,'time2'));
    jd_press=time + (time2/3600/1000/24);
    delta_t = mean(diff(jd_press));
    for v=1:length(ncVars)
        try
        ncobj=squeeze(ncread(nc,ncVars{v}));
        data=ncobj(:);
        try
            theFillValue=ncreadatt(nc,ncVars{v},'_FillValue');
        catch
            theFillValue= 1.0E35;
        end
            data=ncobj;
            bads=find(data>=theFillValue);
            data(bads)=nan;
            dmean_data = data - gmean(data);
            data_int = smart_interp(jd_press,data,jd_int,3);       
            eval([ncVars{v},'_int(:,f) = data_int;'])
        catch
            display ([ncVars{v} ' is not in file ' nc]);
             eval([ncVars{v},'_int(:,f) = nan(size(jd_int));'])
        end
    end   
end
clear data bads dmean_data 
%prune the colormap
div1=floor(length(col)/nFiles);
col1=col([1:4:32],:);
col2=col([39:3:60],:);
cols=[col1; col2];
cols=flipud(cols)

%% now call timeplt.m to create a combined plot for the north site. 
    clf
set(gcf,'position',[100 25 850 1100])%8.5x11
H = tight_subplot(4,1,[.04 .3],[.1 .1],[.1 .2]) ;

H(1)=subplot(4,1,1,'align');
    set(H(1),'nextplot','add')
    set(H(1),'box','on')
    H(2)=subplot(4,1,2,'align');
    set(H(2),'nextplot','add')
    set(H(2),'box','on')
    H(3)=subplot(4,1,3,'align');
    set(H(3),'nextplot','add')
    set(H(3),'box','on')
    H(4)=subplot(4,1,4,'align');
    set(H(4),'nextplot','add')
    set(H(4),'box','on')
    
%%
plot(H(1), jd_int, hsig_int);
[LEGH,OBJH,OUTH,OUTM]=legend(H(1),wavenameFiles,'Location','NorthEast');%,'Orientation','Horizontal');
        set(LEGH,'Fontsize',8);
   title(H(1),'Significant Wave Height Fire Island 2014');
%%
 set(H(2),'ColorOrder',cols);  set(H(3),'ColorOrder',cols);  set(H(4),'ColorOrder',cols);

plot(H(2), jd_int, T_28_int);%legend(H(2),names,'Location','NorthEast','Orientation','Horizontal');
    title(H(2),'Temperature Degrees Celsius');
   plot(H(3), jd_int, C_51_int);%legend(H(3),names,'Location','NorthEast','Orientation','Horizontal');
    title(H(3),'Conductivity');
    plot(H(4), jd_int, S_41_int);
    title(H(4),'Salinity (PSU)');

   % extend xlim to allow for the legend
   xlims(2)=xlims(2)+17; 
for h=1:length(H)
          set(H(h),'xlim',xlims,'xtick',xt)
           set(H(h),'xticklabel',{[datestr(datenum(gregorian(xt)),'mmm dd yyyy')]})
         set(H(h),'ylim',ylims(h,:))
        ylabhan=get(H(h),'Ylabel');
        set(ylabhan,'String',ylab{h})
        pos=get(ylabhan,'position');
        
end

    [LEG,OBJH,OUTH,OUTM]=legend(H(3),names,'Location','NorthEast');%,'Orientation','Horizontal');
%         figCol=get(gcf,'color');
%         set(LEG,'color',figCol)
%         set(LEG,'edgecolor',figCol)
            set(LEG,'Fontsize',8);
         legpos=get(LEG,'position');
        axpos=get(gca,'position');
        set(LEG,'position',[axpos(1)+axpos(3)-.105 axpos(2)+.115 legpos(3) legpos(4)])
        set(LEG,'XAxisLocation','top')
        set(LEG,'ColorOrder',cols)

set(gcf,'PaperPositionMode','Auto','PaperOrientation','portrait')

% give it a title and export to .pdf
title(H(1),{['Salinity and Temperature Observations Fire island 2014']})
print(fullfile(plotDirectory,['ofr2014-xxxx-FI14_ct_sites']),'-dpdf');          




