% river_ts.m  An mfile to plot the timeseries of river discharge data from
%             Van Swamp and Tranters Creek for the North Carolina
%             Open-File report.
%
% 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., 2012, Carolinas Coastal Change Processes Project Data 
%  Report for Observations near Diamond Shoals, NC from February 2010: 
%  U.S. Geological Survey Open-File Report 2012-xxxx.
%
% Program written in Matlab 2012b v8.0.0 
% Program ran on PC with Windows 7 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:
%       julian.m

clear
more off

% define data and output directories
dataDirectory = 'D:\BNA\data\Carolinas\webpage\dataDVD\Hatteras2010DVD\files'; 
%D:\BNA\data\Carolinas\webpage\dataDVD\Hatteras2010DVD\files';
plotDirectory = 'D:\BNA\data\Carolinas\webpage\dataDVD\Hatteras2010DVD\PLOTS';  
%D:\BNA\data\Carolinas\webpage\dataDVD\Hatteras2010DVD\PLOTS';

% define data file names
dataFiles={ 'VanSwamp_02_2012.nc',...
            'TrantersCreek_02_2012.nc' };
nFiles=length(dataFiles);
names={'Van_Swamp','Tranters_Creek'};
% define parameters for plots
startTime=julian(2010,2,1,00);
stopTime=julian(2010,2,28,00);
col=get(0,'defaultaxescolororder');     %plot color order
ylab={ {'Van Swamp','Discharge (ft^3/s)'},...   %plot ylabels
       {'Tranters Creek','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 10;...   %plot ylimits
        0 10 ];
xt=[];
for yy=2010:2010
    for mm=2:2
        for dd=1:5:28
            tic=julian(yy,mm,dd,00);
            if tic>=xlims(1) & tic<=xlims(2)
                xt=[xt; tic];
            end
        end
    end
end
% initialize variables
jd=cell(nFiles,1);
rivDis=cell(nFiles,1);

for f=1:nFiles
    
        river=(fullfile(dataDirectory,dataFiles{f}));
        jd{f}=julian([ncread(river,'Year'),ncread(river,'Month'),...
            ncread(river,'Day'),ncread(river,'Hour'),ncread(river,'Minute'),zeros(length(ncread(river,'Minute')))]);
        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
        
            ncobj = ncread(river,'Gage Height');
            theFillValue = 999;
            eval([names{f},'= ncobj(rividx(1):rividx(2));'])
            eval(['bads = find(',names{f},'(:) >= theFillValue);'])
            eval([names{f},'(bads) = nan;'])
        
        
        % create plot
        xdata = jd{f};
        eval(['ydata = ',names{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),-1,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
