% obs_ts.m  An mfile to plot the timeseries of OBS data from North
%           Carolina.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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
%       smart_interp.m
%       PL66TN.m
%       timeplt.m
%       dolandscape.m
%       suptitle.m

clear
more off

% Define data and output directories
dataDirectory = 'X:\data\hatteras09\BBV_reviewed';
plotDirectory = 'D:\BNA\data\Carolinas\webpage\dataDVD\Hatteras2009DVD\mfiles';

% Define OBS data filenames (OBS data stored w/
% ADV and PCADP data).  Second OBS (OBS 1555) at
% site 6 had no good data.
dataFiles = {'8544pcs-cal.nc','8545advs-cal.nc',...
    '8575advs-cal.nc','8576advs-cal.nc'};
waveFiles={'8542p-cal.nc','8542p-cal.nc',...
    '8572p-cal.nc','8572p-cal.nc'};

nFiles = length(dataFiles);

% Define limits and x-ticks for the plots
startTime=julian(2009,1,8,00);
stopTime=julian(2009,05,15,00);
timediff = stopTime-startTime; % OBS plot xlimits extend 5 percent beyond 'startTime'
xlims = [(startTime-timediff*0.05) (stopTime+timediff*0.05)]; % and 'stopTime'
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
ylims=[0 1]; % OBS plot ylimits


% Initialize structure array for OBS data
obs=struct('time',[],'data',{},'serial',{},'fill',{},'site',[]);

% Load the OBS data
for f = 1:length(dataFiles)

    ncw=netcdf(['X:\data\hatteras09\candidate_final_files\',waveFiles{f},], 'nowrite');
    jd_waves{f} = ncw{'time'}(:) + (ncw{'time2'}(:)/3600/1000/24);
    waveDelta_t=gmean(diff(jd_waves{f}));
    hsig{f}=ncw{'wh_4061'}(:);
    theFillValue=ncw{'wh_4061'}.FillValue_(:);
    badInd=find(hsig{f} >= theFillValue);
    hsig{f}(badInd)=nan;
    badInd=find(hsig{f} <= 0);
    hsig{f}(badInd)=nan;
    ncw=close(ncw);

    nc = netcdf(fullfile(dataDirectory,dataFiles{f}), 'nowrite');
    obs(f).time = nc{'time'}(:) + nc{'time2'}(:)/3600/1000/24;
    % 1st OBS
    ncobj = nc{'NEP1_56'};
    if isempty(ncobj)
        obs(f).data{1} = nan(size(obs(f).time));
        obs(f).serial{1} = '';
    else
        obs(f).data{1} = ncobj(:);
        obs(f).serial{1} = ncobj.serial(:);
        obs(f).fill{1} = ncobj.FillValue_(:);
        bads = find(obs(f).data{1} >= obs(f).fill{1});
        obs(f).data{1}(bads) = nan; % set fill values to NaN
    end
    % 2nd OBS
    ncobj = nc{'NEP2_56'};
    if isempty(ncobj)
        obs(f).data{2} = nan(size(obs(f).time));
        obs(f).serial{2} = '';
    else
        obs(f).data{2} = ncobj(:);
        obs(f).serial{2} = ncobj.serial(:);
        obs(f).fill{2} = ncobj.FillValue_(:);
        bads = find(obs(f).data{2} >= obs(f).fill{2});
        obs(f).data{2}(bads) = nan; % set fill values to NaN
    end
    % Set the site
    obs(f).site = dataFiles{f}(1:4);
end
ncclose
%% Plot the data at North site
figure('pos',[100 100 1000 800])
orient landscape
subplot(3,1,1); % the waves
plot(jd_waves{1},hsig{1})
ylhan=ylabel({'North Site','Significant Wave','Height (m)'});
pos=get(ylhan,'position');
set(gca,'xlim',xlims)
set(gca,'xtick',xt)
set(gca,'xticklabel',[])
hold on;
for p = 1:2 %the obs
    for s = 1:2
        serial = obs(p).serial{s}
        if isempty(serial); serial = 'No OBS'; end
        switch serial
            case 'OBS2182'
                AX(1) = subplot(3,1,2);
            case 'OBS7409'
                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 2182','Voltage (v)'});
        case 2
            ylhan=ylabel({'OBS 7409','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
suptitle({'OBS Observations','North Site'})
print(fullfile(plotDirectory,['Northsite_obs']),'-dpdf');

%% Plot the data at South site
figure('pos',[100 100 1000 800])
orient landscape
subplot(3,1,1); % the waves
plot(jd_waves{3},hsig{3})
ylhan=ylabel({'South Site','Significant Wave','Height (m)'});
pos=get(ylhan,'position');
set(gca,'xlim',xlims)
set(gca,'xtick',xt)
set(gca,'xticklabel',[])
hold on;
for p = 3:4 %the obs
    for s = 1:2
        serial = obs(p).serial{s}
        if isempty(serial); serial = 'No OBS'; end
        switch serial
            case 'OBS2181'
                AX(1) = subplot(3,1,2);
            case 'OBS2183'
                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 2181','Voltage (v)'});
        case 2
            ylhan=ylabel({'OBS 2183','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
suptitle({'OBS Observations','South Site'})
print(fullfile(plotDirectory,['Southsite_obs']),'-dpdf');
