% met_ts.m  An mfile to plot the timeseries of South Carolina
%           air temperature, sea temperature, barometric pressure,
%           significant wave height, and average wave period from NDBC
%           moored buoys, NDBC C-MAN stations, and NOS water level stations.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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 04/27/06, show time series of temperature, pressure, and waves 
%                     
% Dependencies:
%       the netcdf toolbox
%       julian.m
%       timeplt.m
%       dolandscape.m
%       fixpaper.m
%       suptitle.m
%       smart_interp.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\MET';

% define data file names
dataFiles = { '41013h1-stdMet.nc',...    % NDBC moored buoy
              'FPSN7h1-stdMet.nc',...    % NDBC C-MAN station
              '8659897h1-stdMet.nc',...  % NOS water level station
              '8661070h1-stdMet.nc' };   % NOS waterlevel station
nFiles = length(dataFiles);

% define parameters for plots
ncvarnames = {'AT_21','T_25','BP_915','wh_4061','dwp_4063'};
names = {'air_temp','sea_temp','bp','Hs','T'};
startTime = julian(2003,10,25,00);        %plot start stime
stopTime = julian(2004,04,25,00);         %plot stop time
jd_int = startTime:1/24:stopTime;
ndt = length(jd_int);
col = get(0,'defaultaxescolororder');     %plot color order
ylab = { {'Peak Wave','Period (sec)'},...   %plot ylabels
         {'Significant','Wave Height (m)'},...
         {'Barometric','Pressure (mbar)'},...
         {'Sea Surface','Temperature','(degrees C)'},...
         {'Air','Temperature','(degrees C)'} };
timediff = stopTime-startTime;                                %plot xlimits extend 5 percent
xlims = [(startTime-timediff*0.05) (stopTime+timediff*0.05)]; %past start time and stop time
ylims = [ 2 20;...   %plot ylimits
          0 5;...
       1000 1040;...
          0 30;...
          -10 30];
xt = [];
for yy = 2003:2004
    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

% initialize variables
air_temp_all = nan(ndt,nFiles);
sea_temp_all = nan(ndt,nFiles);
bp_all = nan(ndt,nFiles);
Hs_all = nan(ndt,nFiles);
T_all = nan(ndt,nFiles);

for f=1:nFiles
    
        nc = netcdf(fullfile(dataDirectory,dataFiles{f}));
        sta_type = nc.STATION_TYPE(:);
        sta_id = nc.STATION_ID(:);
        station = [sta_type,' ',sta_id]; %station ID
        
        % get time
        jd = nc{'time'}(:) + (nc{'time2'}(:)/3600/1000/24);
        metidx(1) = find(jd(:) >= startTime,1,'first');
        metidx(2) = find(jd(:) <= stopTime,1,'last');
        jd = jd(metidx(1):metidx(2));
        
        % get data, nan the bad stuff, and interpolate onto 
        % standard time array
        for n = 1:length(ncvarnames);
            ncobj = nc{ncvarnames{n}};
            if isempty(ncobj)
                data_int = nan(size(jd_int'));
            else
                theFillValue = ncobj.FillValue_(:);
                data = ncobj(metidx(1):metidx(2));
                bads = find(data >= theFillValue);
                data(bads) = nan;
                data_int = smart_interp(jd, data, jd_int, 3)';
            end
            eval([names{n},'= data_int;'])
        end
        nc = close(nc);

        % call timeplt.m to create individual plots
        % (clearing global timeplt*, the ResizeFcn,
        % and the Tag are necessary to release
        % timeplt.m's control over your figure!)
        xdata = jd_int;
        ydata = [T Hs bp sea_temp air_temp];
        panelOrd = [1 2 3 4 5]; 
        figure
        han = timeplt(xdata,ydata,panelOrd);
        clear global timeplt*
        set(gcf,'ResizeFcn',[])
        set(gcf,'Tag','')
        dolandscape
        fixpaper
 
        % reformat xlims, ylims, ylabels, and line color
        for h = 1:length(han)
            axes(han(h))
            grid on
            set(gca,'xlim',xlims,'xtick',xt)
            set(gca,'ylim',ylims(h,:))
            ylabhan = get(gca,'Ylabel');
            set(ylabhan,'String',ylab{h})
            pos=get(ylabhan,'position');
            set(get(gca,'Ylabel'),'position',[julian(2003,10,7) pos(2) pos(3)]) %align ylabels nicely
            linehan1 = findobj(gca,'Type','line','color','b');
            set(linehan1,'color',col(f,:))
            linehan2 = findobj(gca,'Type','line','color','k');
            if ~isempty(linehan2)
                delete(linehan2)
            end
            if h == 4
                if f == 2 || f == 3 
                    text(julian(2003,11,5),15,'No Data')
                end
            elseif h == 1
                if f == 3 || f == 4
                    text(julian(2003,11,5),12,'No Data')
                end
            elseif h == 2
                if f == 3 || f == 4
                    text(julian(2003,11,5),3,'No Data')
                end
            end                
        end

        % reformat 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)}; %2 line label
            text(xt(tic),1,xtl,'horizontalalignment','center',...
                                 'verticalalignment','top')
        end
        
        % give it a title and export to .pdf
        set(gcf,'PaperPositionMode','Auto')
        suptitle({station,' Temperature, Pressure, and Wave Observations'});
        print(fullfile(plotDirectory,['temp_press_wvs_',sta_id]),'-dpdf')

        % get data ready for combined plot
        air_temp_all(:,f) = air_temp;
        sea_temp_all(:,f) = sea_temp;
        bp_all(:,f) = bp;
        Hs_all(:,f) = Hs;
        T_all(:,f) = T;

end            

% now call timeplt.m to create a combined plot.
ydata =  [T_all Hs_all bp_all sea_temp_all air_temp_all];
panelOrd = [repmat(1,1,nFiles),...
            repmat(2,1,nFiles),...
            repmat(3,1,nFiles),...
            repmat(4,1,nFiles),...
            repmat(5,1,nFiles)];
figure
han = timeplt(xdata, ydata, panelOrd);
clear global timeplt*
set(gcf,'ResizeFcn',[])
set(gcf,'Tag','')
LEG=legend(han(5),strvcat('NDBC Moored','Buoy 41013'),...
                  strvcat('NDBC C-Man','Station FPSN7'),...
                  strvcat('NOS Water Level','Station 8659897'),...
                  strvcat('NOS Water Level','Station 8661070'));
set(LEG,'Location','NorthOutside')   
set(LEG,'Orientation','Horizontal')
figCol = get(gcf,'color');
set(LEG,'color',figCol)
set(LEG,'edgecolor',figCol)
dolandscape
fixpaper
pos=get(LEG,'position');
set(LEG,'position',[pos(1) pos(2)+.1 pos(3) pos(4)])

% reformat xlims, ylims, ylabels
for h = 1:length(han)
    axes(han(h))
    grid on
    set(gca,'xlim',xlims,'xtick',xt)
    set(gca,'ylim',ylims(h,:))
    ylabhan=get(gca,'Ylabel');
    set(ylabhan,'String',ylab{h})
    pos=get(ylabhan,'position');
    set(get(gca,'Ylabel'),'position',[julian([2003,10,7,0,0,0]) pos(2) pos(3)]) %align ylabels nicely
    linehan2 = findobj(gca,'Type','line','color','k');
    if ~isempty(linehan2)
        delete(linehan2)
    end
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)}; %2 line label
    text(xt(tic),1,xtl,'horizontalalignment','center',...
        'verticalalignment','top')
end

% give it a title and export to .pdf
set(gcf,'PaperPositionMode','Auto')
suptitle([{'All Stations','Temperature, Pressure, and Wave Observations'}])
print(fullfile(plotDirectory,['temp_press_wvs_all']),'-dpdf')