% 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.,  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."
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 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_2014\candidate_final_files';
plotDirectory = 'I:\FireIsland\FI_14_WEBOFR\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 = {'9885advs-cal-trm.nc',...
    '9882pcs-cal-trm.nc','9916advs-cal-trmclean.nc','9917advs-cal-trm.nc'};
% waveFiles={'9892awWvs-p.nc',...
%     '9921whp-cal.nc'};
wavenameFiles={'Site 2','Site 2','Site 3','Site 3'};
nFiles = length(dataFiles);

% Define limits and x-ticks for the plots
startTime=julian(2014,2,5,00);
stopTime=julian(2014,05,7,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=2014:2014
    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(3,1,[.04 .3],[.1 .1],[.1 .1]) ;

H(1)=subplot(3,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:2 %the obs
    for s = 1
        serial = obs(p).serial{s}
        if isempty(serial); serial = 'No OBS'; end
        switch serial
            case 'OBS12069'
                AX(1) = subplot(3,1,2);
            case 'OBS2277'
                AX(2) = subplot(3,1,3);
                case 'No OBS'
                pause(1);
           otherwise
                break
        end
        H(1+p)=subplot(3,1,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 12069','Voltage (v)'});
        case 2
            ylhan=ylabel({'OBS 2277','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'});dolandscape;
print(fullfile(plotDirectory,['ofr2014-xxxx-FI14_site2_obs']),'-dpdf');
close
%% Plot the data at site 1
figure('pos',[100 100 1000 800])
orient landscape
H(1)=subplot(3,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 = 3:4 %the obs
    for s = 1
        serial = obs(p).serial{s}
        if isempty(serial); serial = 'No OBS'; end
        switch serial
            case 'OBS14204'
                AX(1) = subplot(3,1,2);
            case 'OBS14205'
                AX(2) = subplot(3,1,3);
                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 14204','Voltage (v)'});
        case 2
            ylhan=ylabel({'OBS 14205','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 3'})
print(fullfile(plotDirectory,['ofr2014-xxxx-FI14_site3_obs']),'-dpdf');