% 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., Voulgaris, George 
% 2014, Coastal Change Processes Project data report for oceanographic 
% observations near Fire Island, New York, January to April  
% 2012 : 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 = '\\IGSAGIEGGS-FILE\Stg-proc\data\FireIsland_2012\candidate_final_files';
plotDirectory = 'I:\FireIsland\FI_12_DataDVD\PLOTS';

% define data file names
dataFiles={'9202sc-a.cdf','9221mc-a.cdf','9232sc-a.cdf',...
    '9262sc-a.cdf','9271mc-a.cdf','9282mc-a-clean.cdf',...
    '9291mc-a.cdf','9302mc-a.cdf','9311mc-a.cdf','9322sc-a.cdf',...
    '9332sc-a.cdf'};
        names={'Site 1f','Site 1s','Site 2f',...
            'Site 3','Site 4s','Site 4',...
            'Site 5s','Site 5','Site 6s','Site 6',...
            'Site 7'};
waveFiles={'9241whp-cal-trm.nc','9261awWvs-p-ll-trm.nc',...
            '9281whp-cal-trm.nc','9303sgdataWvs-p-trm.nc','9321awWvs-p-ll-trm.nc'};
wavenameFiles={'Site 2m','Site 3m','Site 4','Site 5','site 6'};

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

startTime=datenum(gregorian(julian(2012,1,27,00)));                   
stopTime=datenum(gregorian(julian(2012,05,9,00)));         
col=get(0,'defaultaxescolororder');     %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.05) (stopTime+timediff*0.05)]; %past start time and stop time
ylims=[ 0   4;...
    22 28;...   %plot ylimits
    4 14;...
    3 4];
xt=[];
for yy=2012:2012
    for mm=1:12
        tic=datenum(gregorian(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
startTime=julian(2012,1,27,00);                   
stopTime=julian(2012,05,9,00);         
jd_int = [startTime:1/24:stopTime]';
maxGap = 1;

% Load site wave data
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_40=zeros(length(jd_int),nFiles);
STH_71=zeros(length(jd_int),nFiles);

%%
for f=1:nFiles
    f
    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

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

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), datenum(gregorian(jd_int(5:end))), hsig_int(5:end,:));legend(H(1),wavenameFiles,'Location','NorthEast');%,'Orientation','Horizontal');
    %title(H(1),'Significant Wave Height Fire Island 2012');
    plot(H(2), datenum(gregorian(jd_int)), T_28_int);legend(H(2),names,'Location','NorthEast');%,'Orientation','Horizontal');
    %title(H(2),'Temperature Degrees Celsius');
    plot(H(3), datenum(gregorian(jd_int)), C_51_int);legend(H(3),names,'Location','NorthEast');%,'Orientation','Horizontal');
    %title(H(3),'Conductivity');
    plot(H(4), datenum(gregorian(jd_int)), S_41_int);legend(H(4),names,'Location','NorthEast');%,'Orientation','Horizontal');
    %title(H(4),'Salinity (PSU)');

   

for h=1:length(H)
%         axes(H(h))
        set(H(h),'xlim',xlims,'xtick',xt)
        set(H(h),'xticklabel',{[datestr(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
set(gcf,'PaperPositionMode','Auto','PaperOrientation','portrait')

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




