%aggregateTransects.m
%Edwin Kite

%Purpose: 
%(1) Aggregate data gathered from multiple Transects onto
%a single stratigraphic column. 
%(2) Align multiple data types (e.g., meander wavelengths and channel
%widths).
%(3) (TO BE IMPLEMENTED LATER) Launch scripts for producing
%publication-quality figures.

%Requirements: If summary data files do not already exist,
%aggregateTransects will cause them to be generated, launching the appropriate
%m-scripts as needed.

%******************** * CONTROL PARAMETERS ********************************
iss_for_display = 2;%Interpolation of stratigraphic surface. %1 = Planar; 2 = Quadratic; 3 = IDW; 4 = Kriged.
if iss_for_display == 4; error('Error bars for kriging not yet incorporated'); end

%**** CHECK FOR EXISTENCE OF OUTPUT DATA FILES, PROCESSING IF NEEDED ******
for WhichTransect = 1:3
if exist(strcat('outputMeanderWavelengthMeasurementsMarsRivers_Transect_',num2str(WhichTransect),'.mat'),'file') ~= 2
    parseMeanderWavelengthMeasurementsMarsRivers
    save WhichTransect WhichTransect iss_for_display
    clear
    load WhichTransect
end
end

for WhichTransect = 1:3
if exist(strcat('outputWidthMeasurementsMarsRivers_Transect_',num2str(WhichTransect),'.mat'),'file') ~= 2
    parseWidthMeasurementsMarsRivers
    save WhichTransect WhichTransect iss_for_display
    clear
    load WhichTransect
end
end

for WhichTransect = 1:3
if exist(strcat('outputAuxiliaryDataMarsRivers_Transect_',num2str(WhichTransect),'.mat'),'file') ~= 2
    parseAuxiliaryDataMarsRivers
    save WhichTransect WhichTransect iss_for_display
    clear
    load WhichTransect
end
end

%************ END OF CHECK FOR EXISTENCE OF OUTPUT DATA FILES *************

%************ ADDITIONAL CONTROL PARAMETERS *******************************

bigstratbinwidth = 40; %Width of 'big' stratigraphic-elevation bins (m).
bigbinlist = [-300:bigstratbinwidth:300];

%******* AGGREGATE STRATIGRAPHIC INFORMATION ON WIDTHS ********************

loadPostprocessedChannelWidthMeasurementsMarsRivers
    
    awi = 0; %aggregate-width index
    clear aw awerr awlog awn
    
%This might be changed so as to report log errors.
for bigbin = bigbinlist
    awi = awi + 1;
    wtbin = find( stratz > bigbin & stratz < (bigbin + bigstratbinwidth));%widths, this stratigraphic bin

        awn(awi) = length(wtbin);
         aw(awi) = nanmean([wall(wtbin).cw]); %aggregate lambda
      awerr(awi) = sqrt( ( sum(( [wall(wtbin).cwdev    ].^2+[wall(wtbin).cw].^2 )) ...
                    / length(wtbin))...
                    - aw(awi).^2);  
      awlog(awi) = nanmean([wall(wtbin).cwlog]);
   awlogerr(awi) = sqrt( ( sum(( [wall(wtbin).cwdevlog    ].^2+[wall(wtbin).cwlog].^2 )) ...
                    / length(wtbin))...
                    - awlog(awi).^2);  
end


%***************** DISPLAY SUMMARY WIDTH FIGURE ***************************
figure

save iss_for_display iss_for_display bigstratbinwidth bigbinlist aw awerr awlog awn awlogerr awlog

% load outputBootstrapRiverDepositMeasurementsMarsRivers_6_JUNE_2014_TEMPORARY
% 
% plot(bigbinlist, w3l2sig_high,'Color',[.7 .7 .7])
% hold on
% plot(bigbinlist, w3l2sig_low,'Color',[.7 .7 .7])
% plot(bigbinlist, w3lmedian,'b','LineWidth',2)
% 
load outputBootstrapRiverDepositMeasurementsMarsRivers_log_errors_production

plot(bigbinlist, exp(w3l2sig_high),'Color',[.7 .7 .7])
hold on
plot(bigbinlist, exp(w3l2sig_low),'Color',[.7 .7 .7])
plot(bigbinlist, exp(w3lmedian),'b','LineWidth',2)


load iss_for_display


for trs = 1:3
    if trs == 1; trsstring = 'x'; elseif trs == 2; trsstring = 'o'; elseif trs == 3; trsstring = 'd';end

       for icw = 1:length(w(trs).t)
            w(trs).tzadjusted(icw) = [w(trs).t(icw).czs(iss_for_display)];
            w(trs).tzdevadjusted(icw) = [w(trs).t(icw).czsdev(iss_for_display)];
       end
    %h = errorbar(w(trs).tzadjusted,[w(trs).t(:).cw],[w(trs).t(:).cwdev],strcat(trsstring,'k'),'LineWidth',1,'MarkerSize',0.001)
    %h = errorbar(w(trs).tzadjusted,exp([w(trs).t(:).cwlog]),[w(trs).t(:).cwdev],strcat(trsstring,'k'),'LineWidth',1,'MarkerSize',0.001)
    %Added 20 June 2014
    h= errorbar(w(trs).tzadjusted,exp([w(trs).t(:).cwlog]), ...
         exp([w(trs).t(:).cwlog]) - exp([w(trs).t(:).cwlog]-[w(trs).t(:).cwdevlog]), ...
         exp([w(trs).t(:).cwlog] + [w(trs).t(:).cwdevlog]) - exp([w(trs).t(:).cwlog]),strcat('k',trsstring),'MarkerSize',0.001)

    
    
    hold on;grid on;errorbar_tick(h,0);
    set(gca,'FontSize',18)
    title({'Channel widths (full error accounting)'; ...
        'crosses = Transect 1, circles = Transect 2, diamonds = Transect 3'},'Interpreter','Latex','FontSize',18)
    xlabel('Stratigraphic elevation (m)','Interpreter','Latex'); ylabel('Channel width (m)','Interpreter','Latex')
    
    %scatter(w(trs).tzadjusted,[w(trs).t(:).cw],100,strcat(trsstring,'k'),'LineWidth',1)
    %Added 20 June 2014
    scatter(w(trs).tzadjusted,exp([w(trs).t(:).cwlog]),100,strcat(trsstring,'k'),'LineWidth',1)

end
    set(gca,'yscale','log')
    ylim([4 120])
    view([90 90]); set(gca,'xdir','reverse') 
    %errorbar(bigbinlist+(bigstratbinwidth/2),aw,awerr,'xr','LineWidth',3)
    %Added 20 June 2014
    errorbar(bigbinlist+(bigstratbinwidth/2),exp(awlog), ...
       exp(awlog) - exp(awlog-awlogerr), exp(awlog+awlogerr) - exp(awlog) ,'xr','LineWidth',3)

    saveas(gca,strcat('figureAggregateTransects_1.eps'),'epsc');

    save iss_for_display iss_for_display bigstratbinwidth bigbinlist aw awerr awlog awn awlogerr
   
    


%*********** END OF AGGREGATE STRATIGRAPHIC INFORMATION ON WIDTHS *********
clear
load iss_for_display 

%*************** AGGREGATE STRATIGRAPHIC INFORMATION ON WAVELENGTHS *******


loadPostprocessedMeanderWavelengthMeasurementsMarsRivers

ali = 0; %aggregate-lambda index
    clear al alerr allog aln
     
for bigbin = bigbinlist
    ali = ali + 1;
    rhmtbin = find( stratz > bigbin & stratz < (bigbin + bigstratbinwidth));%replicable half-meanders, this stratigraphic bin

        aln(ali) = length(rhmtbin);
         al(ali) = nanmean([lallr(rhmtbin).lambda]); %aggregate lambda
      alerr(ali) = sqrt( ( sum(( [lallr(rhmtbin).err    ].^2+[lallr(rhmtbin).lambda].^2 )) ...
                    / length(rhmtbin))...
                    - al(ali).^2);
                
      allog(ali) = nanmean([lallr(rhmtbin).lambdalog]);
   allogerr(ali) = sqrt( ( sum(( [lallr(rhmtbin).errlog    ].^2+[lallr(rhmtbin).lambdalog].^2 )) ...
                    / length(rhmtbin))...
                    - allog(ali).^2);
end

%***************** DISPLAY SUMMARY WAVELENGTHS FIGURE *********************
figure

% load outputBootstrapRiverDepositMeasurementsMarsRivers_6_JUNE_2014_TEMPORARY
% plot(bigbinlist, 2.*l32sig_high,'Color',[.7 .7 .7])
% hold on
% plot(bigbinlist, 2.*l32sig_low,'Color',[.7 .7 .7])
% plot(bigbinlist, 2.*l3median,'b','LineWidth',2)

load outputBootstrapRiverDepositMeasurementsMarsRivers_log_errors_production
plot(bigbinlist, 2.*exp(l32sig_high),'Color',[.7 .7 .7])
hold on
plot(bigbinlist, 2.*exp(l32sig_low),'Color',[.7 .7 .7])
plot(bigbinlist, 2.*exp(l3median),'b','LineWidth',2)

load iss_for_display %resets 'bigbinlist'


for trs = 1:3
        if trs == 1; trsstring = 'x'; elseif trs == 2; trsstring = 'o'; elseif trs == 3; trsstring = 'd';end
    
        for irhm = 1:length(l(trs).chm)
            l(trs).czadjusted(irhm) = [l(trs).chm(irhm).z] - [l(trs).chm(irhm).cz(iss_for_display)];
        end
        %errorbar([l(trs).czadjusted],2.*[l(trs).chm(:).lambda],2.*[l(trs).chm(:).err],strcat('k',trsstring),'MarkerSize',0.001)
        %Added 20 June 2014
       % errorbar([l(trs).czadjusted],2.*exp([l(trs).chm(:).lambdalog]),2.*[l(trs).chm(:).err],strcat('k',trsstring),'MarkerSize',0.001)
%Added 20 June 2014
        errorbar([l(trs).czadjusted],2.*exp([l(trs).chm(:).lambdalog]), ...
         2.*exp([l(trs).chm(:).lambdalog]) - 2.*exp([l(trs).chm(:).lambdalog]-[l(trs).chm(:).errlog]), ...
         2.*exp([l(trs).chm(:).lambdalog] + [l(trs).chm(:).errlog]) - 2.*exp([l(trs).chm(:).lambdalog]),strcat('k',trsstring),'MarkerSize',0.001)
 
        
        hold on
       
        %scatter([l(trs).czadjusted],2.*[l(trs).chm(:).lambda],100,'k',trsstring)        
        %Added 20 June 2014
        scatter([l(trs).czadjusted],2.*exp([l(trs).chm(:).lambdalog]),100,'k',trsstring)
        
        box on;grid on; set(gca,'yscale','log');ylim([100 2000])
        title({'Meander wavelength versus corrected elevation (full error accounting)'; ...
            'crosses = Transect 1, circles = Transect 2, diamonds = Transect 3'},'Interpreter','Latex')
        set(gca,'FontSize',18);xlabel('Stratigraphic elevation (m)','Interpreter','Latex');
        ylabel('Meander wavelength (m)','Interpreter','Latex')
end
        view([90 90]); set(gca,'xdir','reverse')
%                errorbar(bigbinlist+(bigstratbinwidth/2),2.*al,2.*alerr,'xr','LineWidth',3)        
%Added 20 June 2014
    errorbar(bigbinlist+(bigstratbinwidth/2),2.*exp(allog), ...
       2.*exp(allog) - 2.*exp(allog-allogerr), 2.*exp(allog+allogerr) - 2.*exp(allog) ,'xr','LineWidth',3)
                
                saveas(gca,strcat('figureAggregateTransects_2b.eps'),'epsc');
        



        

%*****************  DISPLAY SUMMARY 'COMPARISON' FIGURE: ******************

figure %Following http://www.mathworks.com/help/matlab/creating_plots/graph-with-multiple-x-axes-and-y-axes.html
                errorbar(bigbinlist+(bigstratbinwidth/2),2.*exp(allog), ...
                    2.*exp(allog) - 2.*exp(allog-allogerr), 2.*exp(allog+allogerr) - 2.*exp(allog),'xr','LineWidth',3)
                ylabel('Meander wavelength (m)', 'FontSize',18,'Interpreter','Latex')
                xlabel('Stratigraphic elevation (m)', 'FontSize',18,'Interpreter','Latex')
                haxes1 = gca; set(gca,'FontSize',16);%set(gca,'xlim',[-200 350]);
                set(haxes1,'XColor','r','YColor','r')
                  for l = 1:length(awn)
                    if aw(l) > 0  
                       text(bigbinlist(l)+(bigstratbinwidth/2)-15,((aw(l)*10)+(2.*al(l)))/2,num2str(aln(l)),'Color','r','FontSize',18,'HorizontalAlignment','Center')
                    else
                       text(bigbinlist(l)+(bigstratbinwidth/2)-15,2.*al(l),num2str(aln(l)),'Color','r','FontSize',18,'HorizontalAlignment','Center')                      
                    end
                  end
                hold on
                haxes1_pos = get(haxes1,'Position')
                    errorbar(bigbinlist+(bigstratbinwidth/2)+5,exp(awlog)*10, ...
                         10*(exp(awlog) - exp(awlog-awlogerr)), 10*(exp(awlog+awlogerr) - exp(awlog)),'xk','LineWidth',3,'Parent',haxes1)

                    for l = 1:length(awn)
                        if al(l)>0
                            text(bigbinlist(l)+(bigstratbinwidth/2)+20,((aw(l)*10)+(2.*al(l)))/2,num2str(awn(l)),'Color','k','FontSize',18,'Parent',haxes1,'HorizontalAlignment','Center')
                        else
                            text(bigbinlist(l)+(bigstratbinwidth/2)+20,aw(l)*10,num2str(awn(l)),'Color','k','FontSize',18,'Parent',haxes1,'HorizontalAlignment','Center')
                        end
                    end
                    
                                  haxes2 = axes('Position',haxes1_pos,...
              'XAxisLocation','top',...
              'YAxisLocation','right',...
              'Color','none');
                    ylabel('Channel width (m)', 'FontSize',18,'Interpreter','Latex','Parent',haxes2)
                    xlabel('Stratigraphic elevation (m)', 'FontSize',18,'Interpreter','Latex','Parent',haxes2)  
                  set(haxes1,'xdir','reverse');set(haxes2,'xdir','reverse')
                  set(haxes2,'yticklabel',1./10.*(str2num(get(haxes1,'yticklabel'))));
                                    set(haxes1,'view',[90 90]); set(haxes2,'view',[90 90]);
                  set(haxes2,'Position',get(haxes1,'Position'))
                  set(haxes2,'xlim',get(haxes1,'xlim')); 
                  set(haxes2,'ylim',get(haxes1,'ylim'));
                  set(gca,'FontSize',16)    
                  
        saveas(gca,strcat('figureAggregateTransects_3.eps'),'epsc');
        
%DISPLAY SUMMARY SINOUSITY FIGURE
%figure
%scatter(stratzc,log10([lallc(:).sinu]-1))

%************************** AGGREGATE AUXILIARY DATA **********************

clear
load iss_for_display 


%AGGREGATE AUXILIARY INFORMATION

for trs = 1:3 %ts = transect
load(strcat('outputAuxiliaryDataMarsRivers_Transect_',num2str(trs)))
    aux(trs).cdseh = cdseh; %l for lambda, i.e. wavelength
    aux(trs).eseh = eseh;
    if exist('mgseh')~=0; aux(trs).mgseh = mgseh;end
    aux(trs).sbinsize = sbinsize;
    aux(trs).srange = srange;
    clear cdseh eseh mgseh 
end

seseh = sum([[aux(1).eseh(:,iss_for_display)]';[aux(2).eseh(:,iss_for_display)]';aux(3).eseh(:,iss_for_display)']); %sum stratigraphic elevation histogram
scdseh = sum([[aux(1).cdseh(:,iss_for_display)]';[aux(2).cdseh(:,iss_for_display)]';aux(3).cdseh(:,iss_for_display)']); %sum channel deposit histogram
smgseh = sum([[aux(1).mgseh(:,iss_for_display)]';[aux(2).mgseh(:,iss_for_display)]']); %sum migration-evidence histogram


%************* DISPLAY SUMMARY AUXILIARY DATA FIGURE **********************
figure;plot([srange(1):sbinsize:srange(2)],(aux(1).eseh(:,iss_for_display)))
hold on
plot([srange(1):sbinsize:srange(2)],aux(2).eseh(:,iss_for_display),'r')
plot([srange(1):sbinsize:srange(2)],aux(3).eseh(:,iss_for_display),'m')
plot([srange(1):sbinsize:srange(2)],seseh,'k')
plot([srange(1):sbinsize:srange(2)],aux(1).cdseh(:,iss_for_display),'b','LineWidth',2)
plot([srange(1):sbinsize:srange(2)],aux(2).cdseh(:,iss_for_display),'r','LineWidth',2)
plot([srange(1):sbinsize:srange(2)],scdseh        ,'g','LineWidth',2)
plot([srange(1):sbinsize:srange(2)],smgseh,'g','LineWidth',1)
title('Area of inspected (HiRISE DTM) outcrop versus stratigraphic elevation','FontSize',18,'Interpreter','Latex')
view([90 90]); set(gca,'xdir','reverse'); ylabel('Area in m^2','FontSize',18);xlabel('Stratigraphic elevation (m)','FontSize',18);
set(gca,'FontSize',18)
    saveas(gca,strcat('figureAggregateTransects_4.eps'),'epsc');

figure
plot([srange(1):sbinsize:srange(2)],scdseh./seseh,'g','LineWidth',2)
hold on
plot([srange(1):sbinsize:srange(2)],smgseh./seseh,'g','LineWidth',1)
view([90 90]); set(gca,'xdir','reverse')
grid on; ylabel({'Areal Fraction of Channel Deposits (thick green line)';'or Areal Fraction of Lateral-Migration Evidence (thin green line)'},'FontSize',18,'Interpreter','Latex'); 
xlabel('Stratigraphic elevation (m)','FontSize',18,'Interpreter','Latex');
set(gca,'FontSize',18)
    saveas(gca,strcat('figureAggregateTransects_5.eps'),'epsc');
