% river_ts.m  An mfile to plot the timeseries of river discharge data from
%             the Cape Fear and Pee Dee Rivers for the North Carolina
%             Open-File report.
%
% Adapted by Brandy Armstrong
% USGS Woods Hole Science Center
% Woods Hole, MA 02540
% barmstrong@usgs.gov
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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
%       julian.m
%       timeplt.m
%       dolandscape.m
%       fixpaper.m
%       suptitle.m

clear
more off

% define data and output directories
dataDirectory = 'D:\BNA\data\Carolinas\webpage\dataDVD\Hatteras2009DVD\supporting data';
plotDirectory = 'D:\BNA\data\Carolinas\webpage\dataDVD\Hatteras2009DVD\PLOTS\MET';

% define data file names
dataFiles={ '02084557_data',...
            '02092500_data' };
nFiles=length(dataFiles);

% define parameters for plots
ncvarnames = {'flow'};
names = {'rivDis'};
startTime=julian(2009,1,8,00);
stopTime=julian(2009,05,15,00);
col=get(0,'defaultaxescolororder');     %plot color order
ylab={ {'Van Swamp','Discharge (ft^3/s)'},...   %plot ylabels
       {'Trent River','Discharge (ft^3/s)'} };
timediff=stopTime-startTime;                                %plot xlimits extend 5 percent
xlims=[(startTime-timediff*0.05) (stopTime+timediff*0.05)]; %past start time and stop time
ylims=[ 0 40;...   %plot ylimits
        0 500 ];
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];
        end
    end
end
% initialize variables
jd=cell(nFiles,1);
rivDis=cell(nFiles,1);

for f=1:nFiles
    
        river=dlmread(fullfile(dataDirectory,dataFiles{f}));
        jd{f}=julian([river(:,2),river(:,3),river(:,4),river(:,5),river(:,6),river(:,7)]);
        rividx(1)=find(jd{f}(:) >= startTime,1,'first');
        rividx(2)=find(jd{f}(:) <= stopTime,1,'last');
        
        %trim time
        jd{f}=jd{f}(rividx(1):rividx(2));
        delta_t=gmean(diff(jd{f}));     %time step of the data
        
        %get data and nan the bad stuff
        for n = 1:length(ncvarnames);
            ncobj = river(:,11);
            theFillValue = 999;
            eval([names{n},'{f} = ncobj(rividx(1):rividx(2));'])
            eval(['bads = find(',names{n},'{f}(:) >= theFillValue);'])
            eval([names{n},'{f}(bads) = nan;'])
        end
        
        % create plot
        xdata = jd{f};
        ydata = rivDis{f};
        subplot(2,1,f)
        plot(xdata, ydata)
        ylabel(ylab{f});
        set(gca,'ylim',ylims(f,:))
        set(gca,'xlim',xlims,'xtick',xt,'xticklabel',[])
        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),-20,xtl,'horizontalalignment','center',...
                                 'verticalalignment','top')
        end
        grid on
        orient landscape
        
        % give it a title and export to .pdf
        if f == nFiles
            suptitle(['Average Daily Streamflow'])
            print(fullfile(plotDirectory,['sg_flow']),'-dpdf')
        end
end
