% obs_ts.m  An mfile to plot the (optionally filtered)
%           timeseries of current velocity data from South
%           Carolina ADV data.    

%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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 03/18/05, Raw ADV 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'.
%                       To plot low-pass filtered OBS data, set the
%                       switch 'filt' to 'yes', otherwise (filt='no')
%                       unfiltered OBS data is plotted.
%
% 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\DATAFILES';
plotDirectory = 'E:\CSULLIVAN\SOUTH_CAROLINA\PLOTS\OBS';

% Define OBS data filenames (OBS data stored w/
% ADV and PCADP data)
dataFiles = {'adv7282vp-cal2.nc','adv7482vp-cal2.nc',...
             'adv7261vp-cal2.nc','adv7461vp-cal2.nc',...
             '728pcvp-cal2.nc','748pcvp-cal2.nc',...
             'adv7262vp-cal2.nc','adv7462vp-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 2]; % OBS plot ylimits

% Define whether to filter and/or interpolate
% data 
filtr = 'no';
intrp = 'no';
switch intrp
    case {'yes','y','Yes','Y'}
        % 'jd_int' is the time array onto which all data is
        % interpolated.
        int_delta_t = 1/24; % 1 hour time spacing
        obs_time_int = [startTime:int_delta_t:stopTime]';
        % data gaps larger than 'maxGap' hours are not
        % interpolated over
        maxGap = 3;
end

% Significant wave height from site 2 is
% shown in the top panel of all plots
wvs1 = netcdf('7221adw-a.nc');
wvs2 = netcdf('7421adw-a.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;

% Initialize variables for OBS data
obs_time = cell(nFiles,1);
obs1_data = cell(nFiles,1);
obs2_data = cell(nFiles,1);
obs_delta_t = cell(nFiles,1);
obs1_serial = cell(nFiles,1);
obs2_serial = cell(nFiles,1);

for f = 1:length(dataFiles)
        % Load the OBS data
        if ~isempty(dataFiles{f})
            obs = netcdf(fullfile(dataDirectory,dataFiles{f}), 'nowrite');
            obs_time{f} = obs{'time'}(:) + obs{'time2'}(:)/3600/1000/24;
            % 1st OBS
            obs1_data{f} = obs{'NEP1_56'}(:);
            if isempty(obs1_data{f})
                obs1_data{f} = nan.*ones(size(obs_time{f}));
            else
                obs1_serial{f} = obs{'NEP1_56'}.serial(:);
                theFillValue = obs{'NEP1_56'}.FillValue_(:);
                badInd = find(obs1_data{f} >= theFillValue);
                obs1_data{f}(badInd) = nan;
            end
            % 2nd OBS
            obs2_data{f} = obs{'NEP2_56'}(:);
            if isempty(obs2_data{f})
                obs2_data{f} = nan.*ones(size(obs_time{f}));
            else
                obs2_serial{f} = obs{'NEP2_56'}.serial(:);
                theFillValue = obs{'NEP2_56'}.FillValue_(:);
                badInd = find(obs2_data{f} >= theFillValue);
                obs2_data{f}(badInd) = nan;
            end
            obs_delta_t{f} = gmean(diff(obs_time{f})); %time step, julian days
            site=obs.DESCRIPT(6); %site #
            obs=close(obs);
        end
            
        if mod(f,2)==0
            
            % Set the site
            if f==2 || f==6
                site = '7A';
            else
                site = '6';
            end
            
            % Combine data from deployment 1 and
            % deployment 2
            obs_time_all = [obs_time{f-1}; nan; obs_time{f}];
            obs1_data_all = [obs1_data{f-1}; nan; obs1_data{f}];
            obs2_data_all = [obs2_data{f-1}; nan; obs2_data{f}];
    
            % Get serial numbers for OBS's
            if strcmp(obs1_serial{f-1},obs1_serial{f})
                serial1 = obs1_serial{f};
            elseif isempty(obs1_serial{f-1}) && ~isempty(obs1_serial{f})
                serial1 = obs1_serial{f};
            elseif ~isempty(obs1_serial{f-1}) && isempty(obs1_serial{f})
                serial1 = obs1_serial{f-1};
            else
                serial1 = 'No 1st OBS';
            end
            if strcmp(obs2_serial{f-1},obs2_serial{f})
                serial2 = obs2_serial{f};
            elseif isempty(obs2_serial{f-1}) && ~isempty(obs2_serial{f})
                serial2 = obs2_serial{f};
            elseif ~isempty(obs2_serial{f-1}) && isempty(obs2_serial{f})
                serial2 = obs2_serial{f-1};
            else
                serial2 = 'No 2nd OBS';
            end
                    
            % (Optionally) Interpolate OBS data
            switch intrp
                case {'yes','y','Yes','Y'}
                    disp('Interpolating wave data')
                    hsig_int = smart_interp(wave_time, hsig, obs_time_int, maxGap);
                    hsig = hsig_int;
                    wave_time = obs_time_int;
                    disp('Interpolating OBS data')
                    obs1_data_all_int = smart_interp(obs_time_all, obs1_data_all, obs_time_int, maxGap);
                    obs1_data_all = obs1_data_all_int;
                    obs2_data_all_int = smart_interp(obs_time_all, obs2_data_all, obs_time_int, maxGap);
                    obs2_data_all = obs2_data_all_int;
                    obs_time_all = obs_time_int;
                    delta_t = int_delta_t;
                otherwise
                    disp('Not interpolating data')
%                     delta_t = obs_delta_t{f}; %TODO:  different delta-t's?
            end
            
%             % (Optionally) Filter OBS data
%             switch filt
%                 case {'yes','y','Yes','Y'}
%                     disp('Low-pass filtering OBS data');
%                     obs1_data_all=PL66TN(obs_data_all, delta_t*24, 33);
%                     obs2_data_all=PL66TN(obs_data_all, delta_t*24, 33);
%                 otherwise
%                     disp('Not filtering data');
%             end
            % Plot
            figure('pos',[100 100 1000 800])
            orient landscape
            AX(1) = subplot(3,1,1);
            plot(wave_time,hsig)
            title('Significant Wave Height from Site 2')
            AX(2) = subplot(3,1,2);
            plot(obs_time_all, obs1_data_all)
            title(serial1)
            AX(3) = subplot(3,1,3);
            plot(obs_time_all, obs2_data_all)
            title(serial2)
            for a = 1:length(AX)
                ax = AX(a);
                % set ylimits and labels
                if a>=2
                    set(ax,'ylim',ylims)
                    set(get(ax,'Ylabel'),'String','OBS voltage (v)')
                else
                    set(ax,'ylim',[0 2])
                    set(get(ax,'Ylabel'),'String',{'Signifiant wave','height (m)'})
                end
                % set xlimits and labels
                set(ax,'xlim',xlims)
                set(ax,'xtick',xticks)
                set(ax,'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.5;
                        text(xticks(t),xtlpos,xtl,'horizontalalignment','center',...
                            'verticalalignment','top')
                    end
                   
                end
                set(ax,'box','on')
                suptitle({['Site ',site,' OBS Observations'],...
                          [dataFiles{f-1},', ',dataFiles{f}]})
            end

%             % combine data from all sites in an array for plotting
%             ind=f/2;
%             for var=1:nVars
%                 eval([varList{var},'_all_comb(:,ind)=',...
%                       varList{var},'_all_int;'])
%             end

        end
        
        
end

% %%%%%%%%%%%%%%%%%%%%%% PLOT %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % decimate the wind
% uv_wind_all_int_dec = nan*ones(size(uv_wind_all_int));
% uv_wind_all_int_dec(1:5:end) = uv_wind_all_int(1:5:end);
% % call timeplt.m 
% xdata=jd_int;
% ydata=[vstd_all_comb,...
%        ustd_all_comb,...
%        v_all_comb,...
%        u_all_comb,...
%        brange_all_comb,...
%        pressure_all_comb]; %uv_wind_all_int_dec];
% panelOrd=[repmat(1,1,nFiles/2),...
%           repmat(2,1,nFiles/2),...
%           repmat(3,1,nFiles/2),...
%           repmat(4,1,nFiles/2),...
%           repmat(5,1,nFiles/2),...
%           repmat(6,1,nFiles/2)];%repmat(6,1,1)];
% figure
% han=timeplt(xdata,ydata,panelOrd,ylims);
% clear global timeplt*
% set(gcf,'ResizeFcn',[])
% set(gcf,'Tag','')
% LEG=legend(han(5),'Site 7 ADV1','Site 7 ADV2','Site 6 ADV1','Site 6 ADV2');
% set(LEG,'Location','NorthOutside')   
% set(LEG,'Orientation','Horizontal')
% figCol=get(gcf,'color');
% set(LEG,'color',figCol)
% set(LEG,'edgecolor',figCol)
% dolandscape
% fixpaper
% legpos=get(LEG,'position');
% pos=get(han(6),'position');
% set(LEG,'position',[pos(1) pos(2)+pos(4) pos(3) legpos(4)])
% 
% % reformat xlims, ylabels
% for h=1:length(han)
%     axes(han(h))
%     set(gca,'xlim',xlims,'xtick',xt)
%     ylabhan=get(gca,'Ylabel');
%     set(ylabhan,'String',ylab{h})
%     pos=get(ylabhan,'position');
%     set(get(gca,'Ylabel'),'position',[julian([2003,10,6,0,0,0]) pos(2) pos(3)])
% end
% 
% % re-do xlabels
% axes(han(1))
% oldhan=findobj(gca,'Tag','xtick labels');
% delete(oldhan);
% for tic=1:length(xt)
%     xtg=gregorian(xt(tic));
%     xtd=datenum(xtg);
%     xtl={[datestr(xtd,3),' 1'];datestr(xtd,10)};
%     xtlpos=ylims(1,1)-1.5;
%     text(xt(tic),xtlpos,xtl,'horizontalalignment','center',...
%         'verticalalignment','top')
% end
% 
% 
% % give it a title and export to .pdf
% switch filt
%     case 'yes'
%         set(gcf,'PaperPositionMode','Auto')
%         suptitle([{'Low-Pass Filtered Wind and Bottom Current Observations',...
%                    'All Sites'}]);
%         print(fullfile(plotDirectory,['adv_all_lp']),'-dpdf');
%     otherwise
%         set(gcf,'PaperPositionMode','Auto')
%         suptitle([{'Wind and Bottom Current Observations',...
%                    'All Sites'}]);
%         print(fullfile(plotDirectory,['adv_all']),'-dpdf');
% end