% river_ts.m  An mfile to plot the timeseries of river discharge data from
%             the Cape Fear and Pee Dee Rivers for the South Carolina
%             Open-File report.
%
% Written by Charlene Sullivan
% USGS Woods Hole Science Center
% Woods Hole, MA 02540
% csullivan@usgs.gov
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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."
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
% Dependencies:
%       the netcdf toolbox
%       julian.m
%       timeplt.m
%       dolandscape.m
%       fixpaper.m
%       suptitle.m

% C. Sullivan   11/17/05,   version 1.0 
% Discharge data downloaded from http://waterdata.usgs.gov/nwis/rt and
% converted to netcdf by rivDis2nc.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\RIVER';

% define data file names
dataFiles={ '2105769d1-river.nc',...
            '2131010d1-river.nc' };
nFiles=length(dataFiles);

% define parameters for plots
ncvarnames = {'flow'};
names = {'rivDis'};
startTime=julian(2003,10,25,00);        %plot start time
stopTime=julian(2004,04,25,00);         %plot stop time  
%startTime=julian(2003,10,25,00);        %plot start time - fall
%stopTime=julian(2004,01,25,00);         %plot stop time - fall
%startTime=julian(2004,01,25,00);        %plot start time - spring
%stopTime=julian(2004,04,25,00);         %plot stop time - spring
col=get(0,'defaultaxescolororder');     %plot color order
ylab={ {'Cape Fear River','Discharge (m^3/s)'},...   %plot ylabels
       {'Pee Dee River','Discharge (m^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 500;...   %plot ylimits
        0 500 ];
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
jd=cell(nFiles,1);
rivDis=cell(nFiles,1);

for f=1:nFiles
    
        nc=netcdf(fullfile(dataDirectory,dataFiles{f}), 'nowrite');
        
        jd{f}=nc{'time'}(:) + (nc{'time2'}(:)/3600/1000/24);
        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 = nc{ncvarnames{n}};
            theFillValue = ncobj.FillValue_(:);
            eval([names{n},'{f} = ncobj(rividx(1):rividx(2));'])
            eval(['bads = find(',names{n},'{f}(:) >= theFillValue);'])
            eval([names{n},'{f}(bads) = nan;'])
        end

        nc=close(nc);
        
        % 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,['pd_cf_flow']),'-dpdf')
        end
end