% obs_ts.m  An mfile to plot the timeseries of OBS data from North
%           Carolina.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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."
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% B. Armstrong 6/20/14, Raw Hydra datalogger data is processed and saved to netCDF files.
%                       It is assumed that the netCDF files are stored in a 
%                       directory defined by the variable 'dataDirectory'.
% Dependencies:
%       smart_interp.m
%       suptitle.m
%       gmean.m
%       replace.m
%       gregorian.m
%       julian.m
%       draw.m

clear
more off

% Define data and output directories
dataDirectory = '\\IGSAGIEGGS-FILE\Stg-proc\data\FireIsland_2012\candidate_final_files';
plotDirectory = 'I:\FireIsland\FI_12_DataDVD\PLOTS';

% Define OBS data filenames (OBS data stored w/
% ADV and PCADP data).  Second OBS (OBS 1555) at
% site 6 had no good data.
dataFiles = {'9203pcs-cal.nc','9205advs-cal.nc','9206advs-cal.nc',...
    '9233pcs-cal.nc','9235advs-cal.nc','9236advs-cal.nc'};
waveFiles={'9241whp-cal-trm.nc','9261awWvs-p-ll-trm.nc',...
    '9281whp-cal-trm.nc','9303sgWvs-p-trm.nc','9321awWvs-p-ll-trm.nc'};
wavenameFiles={'Site 2','Site 3','Site 4','Site 5','site 6'};
nFiles = length(dataFiles);

% Define limits and x-ticks for the plots
startTime=julian(2012,1,27,00);
stopTime=julian(2012,04,19,00);
timediff=stopTime-startTime;                                %plot xlimits extend 5 percent
xlims=[(startTime-timediff*0.05) (stopTime+timediff*0.05)]; %beyond start time and stop time
xt=[];
for yy=2012:2012
    for mm=1:12
        tic=julian(yy,mm,01,00);
        if tic>=xlims(1) & tic<=xlims(2)
            xt=[xt; tic];
        end
    end
end
ylims=[0 1]; % OBS plot ylimits

delta_t = 1/24;
jd_int = (startTime:delta_t:stopTime)';
nTime = length(jd_int);
maxGap = 3;

% Initialize structure array for OBS data
obs=struct('time',[],'data',{},'serial',{},'fill',{},'site',[]);
%% 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
%% Load the OBS data
for f = 1:length(dataFiles)
    nc = (fullfile(dataDirectory,dataFiles{f}));
    obs(f).time = double(ncread(nc,'time')) + (double(ncread(nc,'time2'))/3600/1000/24);
    % 1st OBS
    try
    ncobj = ncread(nc,'NEP2_56');
    obs(f).data{1} = ncobj(:);
        obs(f).serial{1} = ncreadatt(nc,'NEP2_56','serial');
        try
            obs(f).fill{1}=ncreadatt(nc,'NEP2_56','_FillValue');
        catch
            obs(f).fill{1}=1.0E35
        end
        bads = find(obs(f).data{1} >= obs(f).fill{1});
        obs(f).data{1}(bads) = nan; % set fill values to NaN
    catch
        display(['no OBS']);
        obs(f).data{1} = nan(size(obs(f).time));
        obs(f).serial{1} = '';    
    end
    % Set the site
    try
        obs(f).site = ncreadatt(nc,'NEP2_56','initial_sensor_height');
    catch
        obs(f).site = ncreadatt(nc,'NEP1_56','initial_sensor_height');
    end
end
%% Plot the data at site 1
figure('pos',[100 100 1000 800])
H = tight_subplot(4,1,[.04 .3],[.1 .1],[.1 .1]) ;

H(1)=subplot(4,1,1); % the waves
plot(jd_int,hsig_int)
ylhan=ylabel({'Significant Wave','Height (m)','All Sites'});
legend(wavenameFiles,'FontSize',10,'Location','NorthWest','Orientation','Horizontal');
pos=get(ylhan,'position');ylim([0 5]);
set(gca,'xlim',xlims)
set(gca,'xtick',xt)
set(gca,'xticklabel',[])
hold on;
for p = 1:3 %the obs
    for s = 1
        serial = obs(p).serial{s}
        if isempty(serial); serial = 'No OBS'; end
        switch serial
            case 'OBS2181'
                AX(1) = subplot(4,1,2);
            case 'OBS1498'
                AX(2) = subplot(4,1,3);
            case 'OBS1503'
                AX(3) = subplot(4,1,4);
                case 'No OBS'
                pause(1);
           otherwise
                break
        end
        H(1+p)=plot(obs(p).time, obs(p).data{s});
        hold on
    end
end
for a = 1:length(AX)
    axes(AX(a));
    % set xlimits and labels
    set(gca,'xlim',xlims)
    set(gca,'xtick',xt)
    set(gca,'xticklabel',[])
    switch a
        case 1
            ylhan=ylabel({'OBS 2181','Voltage (v)'});
        case 2
            ylhan=ylabel({'OBS 1498','Voltage (v)'});
        case 3
            ylhan=ylabel({'OBS 1503','Voltage (v)'});
    end
    pos=get(ylhan,'position');
    if a == length(AX)
        for t=1:length(xt)
            xtg=gregorian(xt(t));
            xtd=datenum(xtg);
            xtl={[datestr(xtd,3),' 1'];datestr(xtd,10)}; %2 line label
            xtlpos=-0.002;
            text(xt(t),xtlpos,xtl,'horizontalalignment','center',...
                'verticalalignment','top')
        end
    end
    set(gca,'box','on')
end
title(H(1),{'OBS Observations Site 1'});dolandscape;
print(fullfile(plotDirectory,['site1_obs']),'-dpdf');
close
%% Plot the data at site 1
figure('pos',[100 100 1000 800])
orient landscape
H(1)=subplot(4,1,1); % the waves
plot(jd_int,hsig_int)
ylhan=ylabel({'Significant Wave','Height (m)','All Sites'});
legend(wavenameFiles,'FontSize',10,'Location','NorthWest','Orientation','Horizontal');
pos=get(ylhan,'position');ylim([0 5]);
set(gca,'xlim',xlims)
set(gca,'xtick',xt)
set(gca,'xticklabel',[])
hold on;
for p = 4:6 %the obs
    for s = 1
        serial = obs(p).serial{s}
        if isempty(serial); serial = 'No OBS'; end
        switch serial
            case 'OBS2277'
                AX(1) = subplot(4,1,2);
            case 'OBS694'
                AX(2) = subplot(4,1,3);
            case 'OBS1428'
                AX(3) = subplot(4,1,4);
                case 'No OBS'
                pause(1);
           otherwise
                break
        end
        plot(obs(p).time, obs(p).data{s});
        hold on
    end
end
for a = 1:length(AX)
    axes(AX(a));
    % set xlimits and labels
    set(gca,'xlim',xlims)
    set(gca,'xtick',xt)
    set(gca,'xticklabel',[])
    switch a
        case 1
            ylhan=ylabel({'OBS 2277','Voltage (v)'});
        case 2
            ylhan=ylabel({'OBS 694','Voltage (v)'});
        case 3
            ylhan=ylabel({'OBS 1428','Voltage (v)'});
    end
    pos=get(ylhan,'position');
    if a == length(AX)
        for t=1:length(xt)
            xtg=gregorian(xt(t));
            xtd=datenum(xtg);
            xtl={[datestr(xtd,3),' 1'];datestr(xtd,10)}; %2 line label
            xtlpos=-0.002;
            text(xt(t),xtlpos,xtl,'horizontalalignment','center',...
                'verticalalignment','top')
        end
    end
    set(gca,'box','on')
end
title(H(1),{'OBS Observations Site 2'})
print(fullfile(plotDirectory,['site2_obs']),'-dpdf');