% obs_ts.m  An mfile to plot the timeseries of OBS data from South
%           Carolina.    
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Use of this program is described in:
%
% Sullivan, C.M., Warner, J.C., Martini, M.A., Voulgaris, G., 
% Work, P.A., Haas, K.A., and Hanes, D.H. (2006) 
% South Carolina Coastal Erosion Study Data Report for Observations
% October 2003 - April 2004., USGS Open-File Report 2005-1429.
%
% 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."
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

% C. Sullivan 06/23/06, Reformat plot title for consistency w/ the other
%                       plots in the data dvd. Update paths and data file
%                       names.
% C. Sullivan 09/13/05, OBS data is processed and saved to netCDF files by
%                       the USGS Hydratools Toolbox for Matlab.  It is 
%                       assumed that the netCDF files are stored in a 
%                       directory defined by the variable 'dataDirectory'.
%             05/24/06, Update paths to files and file names.  Update to
%                       show Hs from site 2 in the top panel of the plot. 
%                      
% Dependencies:  
%       the netcdf toolbox
%       smart_interp.m
%       PL66TN.m
%       timeplt.m
%       dolandscape.m
%       suptitle.m

clear
more off

% Define data and output directories
dataDirectory = 'E:\CSULLIVAN\SOUTH_CAROLINA\DATA_DVD\DATAFILES';
plotDirectory = 'E:\CSULLIVAN\SOUTH_CAROLINA\DATA_DVD\PLOTS\OBS';

% Define OBS data filenames (OBS data stored w/
% ADV and PCADP data).  Second OBS (OBS 1555) at
% site 6 had no good data.
dataFiles = {'adv7282vp-cal2.nc','adv7482vp-cal2.nc',...
             '728pcvp-cal5.nc','748pcvp-cal6.nc',...
             'adv7261vp-cal2.nc','adv7461vp-cal2.nc'};
nFiles = length(dataFiles);

% Define limits and x-ticks for the plots 
startTime = julian(2003,10,25,00); % OBS plot start time
stopTime = julian(2004,04,25,00); % OBS plot stop time  
timediff = stopTime-startTime; % OBS plot xlimits extend 5 percent beyond 'startTime'
xlims = [(startTime-timediff*0.05) (stopTime+timediff*0.05)]; % and 'stopTime'
xticks = [];
for yy = 2003:2004
    for mm = 1:12
        tick = julian(yy,mm,01,00); 
        if tick >= xlims(1) & tick <= xlims(2)
            xticks = [xticks; tick];   % fir xticks on the first of every month
        end
    end
end
ylims=[0 1]; % OBS plot ylimits

% Significant wave height from site 2 is
% shown in the top panel of all plots
wvs1 = netcdf('7221adwp-cal.nc');
wvs2 = netcdf('7421adwp-cal.nc');
wave_time = [wvs1{'time'}(:) + wvs1{'time2'}(:)/3600/1000/24;...
             wvs2{'time'}(:) + wvs2{'time2'}(:)/3600/1000/24];
hsig = [wvs1{'wh_4061'}(:); wvs2{'wh_4061'}(:)];
theFillValue = wvs1{'wh_4061'}.FillValue_(:);
badInd = find(hsig >= theFillValue);
hsig(badInd) = nan;
close(wvs1)
close(wvs2)

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

% Load the OBS data
for f = 1:length(dataFiles)
        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
        if f <= 4
            obs(f).site = '7A';
        else
            obs(f).site = '6';
        end
end

% Plot the data at site 7A
figure('pos',[100 100 1000 800])
orient landscape
subplot(5,1,1); % the waves
plot(wave_time,hsig) 
ylhan=ylabel({'Site 2','Significant Wave','Height (m)'});
pos=get(ylhan,'position');
set(get(gca,'Ylabel'),'position',[julian(2003,10,8) pos(2) pos(3)]) %align ylabels nicely
set(gca,'xlim',xlims)
set(gca,'xtick',xticks)
set(gca,'xticklabel',[])
for p = 1:4 %the obs
    for s = 1:2
        serial = obs(p).serial{s};
        if isempty(serial); serial = 'No OBS'; end
        switch serial
            case 'OBS1244'
                AX(1) = subplot(5,1,2);
            case 'OBS2277'
                AX(2) = subplot(5,1,3);
            case 'OBS2181'
                AX(3) = subplot(5,1,4);
            case 'OBS2182'
                AX(4) = subplot(5,1,5);
            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',xticks)
    set(gca,'xticklabel',[])
    switch a
        case 1
           ylhan=ylabel({'OBS 1244','Voltage (v)'});
        case 2
           ylhan=ylabel({'OBS 2277','Voltage (v)'});
        case 3
           ylhan=ylabel({'OBS 2181','Voltage (v)'});
        case 4
            ylhan=ylabel({'OBS 2182','Voltage (v)'});
    end
    pos=get(ylhan,'position');
    set(get(gca,'Ylabel'),'position',[julian(2003,10,8) pos(2) pos(3)]) %align ylabels nicely
    if a == length(AX)
        for t=1:length(xticks)
            xtg=gregorian(xticks(t));
            xtd=datenum(xtg);
            xtl={[datestr(xtd,3),' 1'];datestr(xtd,10)}; %2 line label
            xtlpos=-0.002;
            text(xticks(t),xtlpos,xtl,'horizontalalignment','center',...
                'verticalalignment','top')
            text(julian(2004,3,1),0.025,'No Data','horizontalalignment',...
                 'center')
        end
    end
    set(gca,'box','on')
end
suptitle({'OBS Observations','Site 7A'})
print(fullfile(plotDirectory,['site7A_obs']),'-dpdf');

% Plot the data at site 6
figure('pos',[100 100 1000 800])
orient landscape
subplot(2,1,1); % the waves
plot(wave_time,hsig) 
ylabel({'Site 2','Significant Wave','Height (m)'})
set(gca,'xlim',xlims)
set(gca,'xtick',xticks)
set(gca,'xticklabel',[])
for p = 5:6 %the obs
    serial = obs(p).serial{1};
    subplot(2,1,2);
    plot(obs(p).time, obs(p).data{1});
    hold on
    ylhan=ylabel({'OBS 1551','Voltage (v)'});
end
pos=get(ylhan,'position');
set(get(gca,'Ylabel'),'position',[julian(2003,10,8) pos(2) pos(3)]) %align ylabels nicely
% set xlimits and labels
set(gca,'xlim',xlims)
set(gca,'xtick',xticks)
set(gca,'xticklabel',[])
if a == length(AX)
    for t=1:length(xticks)
        xtg=gregorian(xticks(t));
        xtd=datenum(xtg);
        xtl={[datestr(xtd,3),' 1'];datestr(xtd,10)}; %2 line label
        xtlpos=-0.05;
        text(xticks(t),xtlpos,xtl,'horizontalalignment','center',...
            'verticalalignment','top')
        text(julian(2003,12,1),1.5,'No Data','horizontalalignment',...
            'center')
    end
end
set(gca,'box','on')
suptitle({'OBS Observations','Site 6'})
print(fullfile(plotDirectory,['site6_obs']),'-dpdf');