function monthly_stats_rivers(dataDirectory, plotDirectory)
%
%   monthly_stats_rivers(dataDirectory, plotDirectory)
%
%   Code to create the plots/tables of monthly statistics for the
%   river streamflow observations.

%Soupy Alexander, 4/7/04

%Do the scalar variables
theFiles = {'charles_modified.xls'; 'ipswich_modified.xls'; 'merrimack_modified.xls'; ...
        'parker_modified.xls'};
varNames = {'Charles Discharge'; 'Ipswich Discharge'; 'Merrimack Discharge'; ...
        'Parker Discharge'};
fileNames = {'charles'; 'ipswich'; 'merrimack'; 'parker'};
theUnits = {'ft^3/s'; 'ft^3/s'; 'ft^3/s'; 'ft^3/s'};
theLimits = {[0 3000]; [0 3000]; [0 55000]; [0 500];};
monthList = {'Jan'; 'Feb'; 'Mar'; 'Apr'; 'May'; 'Jun'; 'Jul'; 'Aug'; 'Sep'; 'Oct'; ...
        'Nov'; 'Dec'};

for indexFile = 1:length(theFiles)
    [theNums, theText] = xlsread(fullfile(dataDirectory, theFiles{indexFile}));
    theNums(1,:) = [];
    theTime = datenum2julian(theNums(:,1) + datenum('30-Dec-1899'));
    theData = theNums(:,2);
    
    gregTime = gregorian(theTime);
    post1990 = find(gregTime(:,1) >=1990);
    theTime = theTime(post1990);
    gregTime = gregTime(post1990, :);
    theData = theData(post1990);
    good = find(~isnan(theData));
    
    yearList = min(gregTime(:,1)):max(gregTime(:,1));
    
    %Figure out the monthly means
    for indexMonth = 1:12
        inMonth = find(gregTime(:,2) == indexMonth);
       
        if isempty(inMonth)
            monthPointsAll(indexMonth) = 0;
            continue
        end
        
        monthMeanAll(indexMonth) = nanmean(theData(inMonth));
        monthAnomAll(indexMonth) = 0;
        monthSTDAll(indexMonth) = nanstd(theData(inMonth));
        monthMinAll(indexMonth) = nanmin(theData(inMonth));
        monthMaxAll(indexMonth) = nanmax(theData(inMonth));
        monthPointsAll(indexMonth) = length(find(~isnan(theData(inMonth))));
    end
    
    yearMeanAll = nanmean(theData);
    
    monthMeanAll(13) = nanmean(theData);
    monthAnomAll(13) = 0;
    monthSTDAll(13) = nanstd(theData);
    monthMinAll(13) = nanmin(theData);
    monthMaxAll(13) = nanmax(theData);
    monthPointsAll(13) = length(find(~isnan(theData)));
    
    monthMean = repmat(NaN, length(yearList), 13);
    monthAnom = repmat(NaN, length(yearList), 13);
    monthSTD = repmat(NaN, length(yearList), 13);
    monthMin = repmat(NaN, length(yearList), 13);
    monthMax = repmat(NaN, length(yearList), 13);
    monthPoints = repmat(NaN, length(yearList), 13);
    
    
    %Do the monthly averages
    for indexYear = 1:length(yearList)
        targetYear = yearList(indexYear);
        
        inYear = find(gregTime(:,1) == targetYear);
        
        if isempty(inYear)
            monthPoints(indexYear, 13) = 0;
            continue
        end
        
        
        timeInYear = theTime(inYear);
        gregTimeInYear = gregorian(timeInYear);
        dataInYear = theData(inYear);
        
        monthMean(indexYear, 13) = nanmean(dataInYear);
        monthAnom(indexYear, 13) = nanmean(dataInYear) - yearMeanAll;
        monthSTD(indexYear, 13) = nanstd(dataInYear);
        monthMin(indexYear, 13) = nanmin(dataInYear);
        monthMax(indexYear, 13) = nanmax(dataInYear);
        monthPoints(indexYear, 13) = length(find(~isnan(dataInYear)));
        
        for indexMonth = 1:12
            inMonth = find(gregTimeInYear(:,2) == indexMonth);
            if isempty(inMonth)
                monthPoints(indexYear, indexMonth) = 0;
                continue
            end
            
            monthMean(indexYear,indexMonth) = nanmean(dataInYear(inMonth));
            monthAnom(indexYear, indexMonth) = nanmean(dataInYear(inMonth)) - ...
                monthMeanAll(indexMonth);
            monthSTD(indexYear, indexMonth) = nanstd(dataInYear(inMonth));
            monthMin(indexYear, indexMonth) = nanmin(dataInYear(inMonth));
            monthMax(indexYear, indexMonth) = nanmax(dataInYear(inMonth));
            monthPoints(indexYear, indexMonth) = length(find(~isnan(dataInYear(inMonth))));
        end
    end
    
    %Do the actual table creation/plotting
    fID = fopen(fullfile(plotDirectory, [fileNames{indexFile} '.txt']), 'w');
    fprintf(fID, [varNames{indexFile} '\n\n']);
    spaces = '\t\t\t\t\t\t';
    titles = '# Points\tMean\tAnomaly\tSTD\tMin\tMax\t';
    fprintf(fID, ['\t' repmat('Jan\t', 1, 7) repmat('Feb\t', 1, 7) ...
            repmat('Mar\t', 1, 7) repmat('Apr\t', 1, 7) repmat('May\t', 1, 7) ...
            repmat('Jun\t', 1, 7) repmat('Jul\t', 1, 7) repmat('Aug\t', 1, 7) ...
            repmat('Sep\t', 1, 7) repmat('Oct\t', 1, 7) repmat('Nov\t', 1, 7) ...
            repmat('Dec\t', 1, 7) repmat('Year\t', 1, 7) '\n']);
    fprintf(fID, ['\t' repmat(titles, 1, 13) '\n']);
    
    figure, orient landscape
    colList = {'b'; 'g'; 'y'; 'm'; 'k'; 'c'};
    colList = repmat(colList(:), 10, 1);
    symList = {'.'; 'o'; 'x'; '+'; 's'; 'd'; 'v'; 'p'};
    symList = repmat(symList(:), 10, 1);
    A = subplot(2,1,1); hold on
    B = subplot(2,1,2); hold on
    
    for indexYear = 1:length(yearList)
        
        if monthPoints(indexYear, 13) == 0
            continue
        end
        
        fprintf(fID, [num2str(yearList(indexYear), '%.0f') '\t']);
        
        for indexMonth = 1:13
            fprintf(fID, [num2str(monthPoints(indexYear, indexMonth), '%.0f') '\t']);
            if monthPoints(indexYear, indexMonth) ~= 0
                fprintf(fID, [num2str(monthMean(indexYear, indexMonth), '%.2f') '\t' ...
                        num2str(monthAnom(indexYear, indexMonth), '%.2f') '\t' ...
                        num2str(monthSTD(indexYear, indexMonth), '%.2f') '\t' ...
                        num2str(monthMin(indexYear, indexMonth), '%.2f') '\t' ...
                        num2str(monthMax(indexYear, indexMonth), '%.2f') '\t']);
            else
                fprintf(fID, '---\t---\t---\t---\t---\t');
            end
            
            if indexMonth ~= 13    
                subplot(2,1,1)
                plot(indexMonth, monthMean(indexYear, indexMonth), ...
                    [colList{indexYear} symList{indexYear}]);
                subplot(2,1,2)
                plot(indexMonth, monthAnom(indexYear, indexMonth), ...
                    [colList{indexYear} symList{indexYear}]);
            end
        end
        fprintf(fID, '\n');
    end
    
    
    fprintf(fID, 'All Years\t');
    
    for indexMonth = 1:13
        fprintf(fID, [num2str(monthPointsAll(indexMonth), '%.0f') '\t']);
        if monthPointsAll(indexMonth) ~= 0
            fprintf(fID, [num2str(monthMeanAll(indexMonth), '%.2f') '\t' ...
                    num2str(monthAnomAll(indexMonth), '%.2f') '\t' ...
                    num2str(monthSTDAll(indexMonth), '%.2f') '\t' ...
                    num2str(monthMinAll(indexMonth), '%.2f') '\t' ...
                    num2str(monthMaxAll(indexMonth), '%.2f') '\t']);
        else
            fprintf(fID, '---\t---\t---\t---\t---\t');
        end
        
        if indexMonth ~= 13    
            subplot(A)
            plot(indexMonth, monthMeanAll(indexMonth), 'r*');
        end
    end
    
    fclose(fID);
    
    subplot(A)
    xlim([0 13])
    ylim(theLimits{indexFile})
    set(gca, 'xtick', [1:12], 'xticklabel', monthList)
    ylabel(theUnits{indexFile})
    theSize = get(gca, 'position');
    theSize(3) = theSize(3)*.89;
    set(gca, 'position', theSize);
    title('Monthly Mean')
    box('on')
    
    subplot(B)
    xlim([0 13])
    ylim([-1*max(abs(ylim)) max(abs(ylim))])
    set(gca, 'xtick', [1:12], 'xticklabel', monthList)
    ylabel(theUnits{indexFile})
    theSize1 = get(gca, 'position');
    theSize1(3) = theSize1(3)*.89;
    set(gca, 'position', theSize1);
    title('Monthly Anomaly')
    box('on')
    
    legSize(1) = theSize1(1) + theSize1(3)*1.01;
    legSize(2) = theSize1(2);
    legSize(3) = theSize1(3)*0.09;
    legSize(4) = (theSize(2)+theSize(4)) - theSize1(2);
    theLeg = subplot('position', legSize);
    hold on
    
    for index = 1:length(yearList)
        plot(1, -1* yearList(index), [colList{index} symList{index}]);
        t = text(3, -1* yearList(index), num2str(yearList(index), '%.0f'));
        set(t, 'fontsize', 8);
    end
    plot(1, -1*(max(yearList)+1), 'r*')
    t = text(3, -1*(max(yearList)+1), 'All Years');
    set(t, 'fontsize', 8);
    axis('off')
    box('off')
    xlim([0 5])
    ylim([-1*(max(yearList)+2) -1*(min(yearList)-1)])
    
    suptitle(varNames{indexFile})
    
    print('-dpdf', fullfile(plotDirectory, fileNames{indexFile}))
    close
end