%bracketing_crosscheck.m
%Edwin Kite

%Purpose: Compare meander widths and wavelengths when these have been
%measured for the same channel. Based on results reported for Earth, the
%width/wavelength ratio shouldn't change dramatically with elevation and
%the ratio should be about 15 (depends on slope, grainsize, e.t.c.)

%Desired future changes: (1) Make everything part of a 'bc' (bracketing crosscheck) structure. 
%(2) Incorporate all Transects into the comparison.
%(3) Show All Data (not just bracketed data) histograms on the common axes.
%See 4 June 2014 staple-in inside 'Rivers 2' notebook.

%Queries on hand-inspection of output (most of these have since been accounted for):
% Ccid 88 didn't register, why
% ccid 87,88 didn't *show*, why
% ccid 85, why the big blue circle, why offset
% ccid 84, bp id 185, seems like a spurious bank pair
% ccid 66, bp id 193, big blue circle is at the center of a bnka not the center of a bank pair
% ccid 49, appears bracketed, why no show?
% ccid 48, is the reason for the big blue circle misplaced meanx meany? should be median
% ccid 42, why no show?
% ccid 35, appears bracketed, why no show?
% ccid 31, appears bracketed, why no show?
% ccid 30, why no show?
% ccid 27, bankpair ID 133, twisted, manually scrub this one!
% ccid 25 and 24, why no show?
% ccid 11, appears bracketed, why no show?
% ccid 10, bp id 47, blue circle appears misplaced
% ccid 8 & 7, appears bracketed, why no show?
% ccid 1, why no show?

savedebugfigures = 1;
showcheckplot = 0;
iss_for_display = 2;

WhichTransect = 3;
loadWidthMeasurementsMarsRivers

%Identify the wavelength associated with width measurements that "bracketed" a meander
load(strcat('outputWidthMeasurementsMarsRivers_Transect_',num2str(WhichTransect)));
load(strcat('outputMeanderWavelengthMeasurementsMarsRivers_Transect_',num2str(WhichTransect)));

assert(sum([ismember({'rhm','d','chm'},who)]) == 3,'Meander wavelength data file incomplete')

%
mw = [rhm(:).lambda]; %meander wavelength
ccid = [rhm(:).m];

with_unit_assignments = shaperead(strcat('Transect_',num2str(WhichTransect),'_Pick_1_Fourth_Pass_Centerlines_for_Meander_Wavelength_Measurements.shp'));
    
    assert((length(d(1).ctl) == length(with_unit_assignments)));
    ccu = [with_unit_assignments.UnitAssign]; %channel-centerline unit assignment
    mun = ccu([rhm(:).m]); %meander unit
clear with_unit_assignments
%


%Align origins
assert(origin(1) == minx); assert(origin(2) == miny);

%Align coordinate systems
%assert(sum(d(1).origin < origin) == 2)
px = px - (d(1).origin(1) - origin(1)); %Points along banks (for width measurements)
py = py - (d(1).origin(2) - origin(2));

figure
scatter(px,py,10,'k')
hold on
scatter([d(1).vx],[d(1).vy],10,'r')


%For each channel width, determine if there is a meander-wavelength point
%inside the convex hull. Only bother with first pick:

    d(1).bracketed = zeros(size(d(1).vx));
    widthlambda    = zeros(size(d(1).vx));
    widthlambdaz   = zeros(size(d(1).vx));
    lambda(1:(max(bid/2))) = 0;
 lambdaerr(1:(max(bid/2))) = 0;
   lambdac(1:(max(bid/2))) = 0;
lambdacerr(1:(max(bid/2))) = 0;
   zlambda(1:(max(bid/2))) = 0;
      acid(1:(max(bid/2))) = 0;
      unitfromlambda(1:(max(bid/2))) = 0;

for nn = 1:(max(bid)/2) %Channel consisting of two banks
    if s(11).cqual(nn) <=4 %Discard 'candidate' (=low) quality data

    thesebanks = find(bid==((nn*2) ) | bid==((nn*2)-1));
        w(nn) = s(11).cw(nn); %Using the "uniform dotting" width measurement strategy
        
        
       cx(nn) = nanmedian(px(bid==(nn*2) | bid==((nn*2)-1)));
       cy(nn) = nanmedian(py(bid==(nn*2) | bid==((nn*2)-1)));
       
    if isempty(thesebanks)==0
        for m = 1:max(size(d(1).p)) %For each channel

        this_convex_hull = convhull(px(thesebanks),py(thesebanks));
        IN = inpolygon(d(1).p(m).x,d(1).p(m).y, ...
             px(thesebanks(this_convex_hull)), ...
             py(thesebanks(this_convex_hull)));
               d(1).p(m).bracketed = IN;
        if sum(IN)>0
            if median(d(1).p(m).arhm(IN)) > -1
                  ahm(nn) = median(d(1).p(m).arhm(IN)); %associated half-meander
                 acid(nn) = m; %associated channel ID
               lambda(nn) = rhm(   median(d(1).p(m).arhm(IN))).lambda; %Closest (measured along channel centerline) half-meander
            lambdaerr(nn) = rhm(   median(d(1).p(m).arhm(IN))).err;    %Closest (measured along channel centerline) half-meander            
              zlambda(nn) = rhm(   median(d(1).p(m).arhm(IN))).z;
              %lzstrat(nn) = rhm(   median(d(1).p(m).arhm(IN))).cz(iss_for_display);
       unitfromlambda(nn) = mun(ahm(nn)); %unit assignment for this channel
               widthlambda(thesebanks) = rhm(median(d(1).p(m).arhm(IN))).lambda;
              widthlambdaz(thesebanks) = rhm(median(d(1).p(m).arhm(IN))).z;
              
              
              lambdac(nn) = chm(m).lambda;                             %Channel average
           lambdacerr(nn) = chm(m).err;                                %Channel average
            end
        end
        
        end
    end
    end
end



w(133) = NaN; lambda(133) = NaN; lambdac(133) = NaN; 
disp('Manually excluding bank-pair ID 133, from channel-ID 27! See notes from 21FEB2014, Rivers notebook')

%Check plot
if showcheckplot == 1
hold on
for m = 1:max(size(d(1).p))
figure
    scatter(d(1).p(m).x,d(1).p(m).y,5,max(1,d(1).p(m).arhm - min(d(1).p(m).arhm)),'o')
    axis equal; set(gca,'FontSize',18); hold on
text(mean(d(1).p(m).x),mean(d(1).p(m).y),num2str(m),'FontWeight','bold','FontSize',18)
xlimits = get(gca,'xlim');
ylimits = get(gca,'ylim');
scatter(px,py,'k.');
set(gca,'xlim',xlimits)
set(gca,'ylim',ylimits)
box on
scatter(px(widthlambdaz~=0),py(widthlambdaz~=0),30,'k','Filled'),%widthlambda(widthlambdaz~=0)
scatter(cx(lambda>0),cy(lambda>0),100,lambda(lambda>0).*2./w(lambda>0),'Filled')
    colorbar
    for n = 1:length(cx)
        if acid(n) == m
     text(cx(n),cy(n),{strcat('bankpair ID:',num2str(n)); ...
                       strcat('width(m):',num2str(w(n))); ...
                       strcat('widtherr(m):',num2str(s(11).cwdev(n))); ...
                       strcat('lambda(m):',num2str(lambda(n))); ...
                       strcat('lambdaerr(m):',  num2str(lambdaerr(n)))}, ...
                       'FontAngle','oblique','FontSize',12,'Color','r');
        end
    end
    xlabel('distance x (m)'); ylabel('distance y (m)');
    title(strcat('Bracketing Crosscheck: Channel Centerline ID:',num2str(m)))
    if savedebugfigures == 1 
            %save to eps color file
            saveas(gca, ...
            strcat('DebugParseTransect1BracketingCrosscheck',num2str(m),'.eps'), ...
            'epsc');
    end
end
end



[b bint] = regress(lambda(lambda>0)',w(lambda>0)')

figure
hold on
for k = 1:length(w)
    if lambda(k) > 0 & lambdac(k) ~= lambda(k)
        line([w(k) w(k)],[2.*lambda(k) 2.*lambdac(k)],'Color','k')
    end
end
scatter(w(lambda>0)',2.*lambdac(lambda>0)',100,zlambda(lambda>0)','MarkerFaceColor','w')
scatter(w(lambda>0)',2.*lambda(lambda>0)',100,zlambda(lambda>0)','Filled')
line([10 60],[110 660],'Color','k')

set(gca,'FontSize',18)
grid on; colorbar
text(2.5+w(lambdac>0)',2.*lambdac(lambdac>0)',num2str(acid(lambdac>0)'))
xlabel('width (m)');ylabel('meander wavelength (m)');box on
title({'Closed circles: individual replicable half meanders.';
       'Open circles: Channel averages.'})
xlim([0 65])

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

figure
hold on
for k = 1:length(w)
    if lambda(k) > 0 & lambdac(k) ~= lambda(k)
        line([w(k) w(k)],[2.*lambda(k) 2.*lambdac(k)],'Color','k')
    end
end
scatter(w(lambda>0)',2.*lambdac(lambda>0)',100,unitfromlambda(lambda>0)','MarkerFaceColor','w')
scatter(w(lambda>0)',2.*lambda(lambda>0)',100,unitfromlambda(lambda>0)','Filled')
line([10 120],[10 120].*11,'Color','g')
line([10 120],[10 120].*8,'Color','g')
line([10 120],[10 120].*15,'Color','g')


set(gca,'FontSize',18)
grid on; colorbar
text(2.5+w(lambdac>0)',2.*lambdac(lambdac>0)',num2str(acid(lambdac>0)'))
xlabel('width (m)');ylabel('meander wavelength (m)');box on
title({'Closed circles: individual replicable half meanders.';
       'Open circles: Channel averages.'})
xlim([0 110])

saveas(gca, ...
            strcat('DebugParseTransect1BracketingCrosscheckAWithUnits.eps'), ...
            'epsc');       
        
        
        
figure;
plot(sort(2.*lambdac(lambdac>0)./w(lambdac>0)),'r')
hold on;grid on
plot(sort(2.*lambda(lambdac>0)./w(lambdac>0)),'b')
[Y,I] = sort(2.*lambdac(lambdac>0)./w(lambdac>0));
ok_lambda = lambda(lambdac>0);
ok_w = w(lambdac>0);
ok_cpstyle = [s(11).cpstyle(lambdac>0)];
scatter([1:length(ok_w)],sort(2.*lambdac(lambdac>0)./w(lambdac>0)),200,log(ok_cpstyle(I))'./log(2),'Filled')
title('R values versus preservation style'); xlabel('half-meander #');ylabel('R')     
        
figure;
errorbar(zlambda(lambda>0)',2.*lambda(lambda>0)'./w(lambda>0)',2.*lambdaerr(lambda>0)'./w(lambda>0)','.')
hold on; colorbar; set(gca,'FontSize',18)
scatter(zlambda(lambda>0)',2.*lambda(lambda>0)'./w(lambda>0)',100,w(lambda>0)','Filled')
text(10+zlambda(lambda>0)',2.*lambda(lambda>0)'./w(lambda>0)',num2str(acid(lambda>0)'))
set(gca,'FontSize',18)
xlim([-2400 -2050])
grid on
xlabel('elevation (m)');ylabel('meander wavelength / channel width '); title('color is width (m)'); box on
            saveas(gca, ...
            strcat('DebugParseTransect1BracketingCrosscheckB.eps'), ...
            'epsc');

        %R versus elevation style, with units assigned
        figure;
plot(sort(2.*lambdac(lambdac>0)./w(lambdac>0)),'k')
%plot(sort(2.*lambda(lambdac>0)./w(lambdac>0)),'k')
hold on;grid off
line([0 43],[16 16],'Color','g')
line([0 43],[8 8],'Color','g')
set(gca,'FontSize',34);set(gca,'xminortick','on');set(gca,'yminortick','on')
[Y,I] = sort(2.*lambdac(lambdac>0)./w(lambdac>0));
ok_lambda = lambda(lambdac>0);
ok_w = w(lambdac>0);
ok_cpstyle = s(11).cpstyle(lambdac>0);
ok_unitfromlambda = unitfromlambda(lambdac>0);
scatter([1:length(ok_w)],sort(2.*lambdac(lambdac>0)./w(lambdac>0)),100,ok_unitfromlambda(I)','Filled')
scatter(3,30,100,2,'Filled'); text(5,30,'F2','Interpreter','latex','FontSize',34)
scatter(3,35,100,1,'Filled'); text(5,35,'F1','Interpreter','latex','FontSize',34)
text(2,14,{'criterion for paleodischarge'},'Interpreter','latex','FontSize',24)
text(2,10,{'(Burr et al. 2010)'},'Interpreter','latex','FontSize',24)

axis square

caxis([1 4]);xlim([0 43])
title('$\lambda$ / w values versus unit','Interpreter','latex'); xlabel('half-meander \#','Interpreter','latex');ylabel('$\lambda$ / w','Interpreter','latex')     
            saveas(gca, ...
            strcat('DebugParseTransect1BracketingCrosscheck_R_versus_units_assigned.eps'), ...
            'epsc');
  

figure;
errorbar(zlambda(lambda>0)',2.*lambda(lambda>0)'./w(lambda>0)',2.*lambdaerr(lambda>0)'./w(lambda>0)','.')
hold on; colorbar; set(gca,'FontSize',18)
scatter(zlambda(lambda>0)',2.*lambda(lambda>0)'./w(lambda>0)',100,w(lambda>0)','Filled')
text(10+zlambda(lambda>0)',2.*lambda(lambda>0)'./w(lambda>0)',num2str(acid(lambda>0)'))
set(gca,'FontSize',18)
xlim([-2400 -2050])
grid on
xlabel('elevation (m)');ylabel('meander wavelength (m)'); title('color is width (m)'); box on

        
        
%Raw elevation:
figure;scatter(zlambda(lambda>0)',2.*lambda(lambda>0)'./w(lambda>0)',100,log(s(11).cpstyle(lambda>0)')./log(2),'Filled')
grid on; colorbar; set(gca,'FontSize',18)
text(10+zlambda(lambda>0)',2.*lambda(lambda>0)'./w(lambda>0)',num2str(acid(lambda>0)'))
xlim([-2400 -2050])
xlabel('raw elevation (m)');ylabel('meander wavelength / channel width'); title('color is log2(preservation style)');box on
            saveas(gca, ...
            strcat('DebugParseTransect1BracketingCrosscheckC.eps'), ...
            'epsc');

figure;scatter(zlambda(lambda>0)',2.*lambda(lambda>0)'./w(lambda>0)',100,s(11).cqual(lambda>0)','Filled')
grid on; colorbar; set(gca,'FontSize',18)
text(10+zlambda(lambda>0)',2.*lambda(lambda>0)'./w(lambda>0)',num2str(acid(lambda>0)'))
%xlim([-2400 -2050])
xlabel('raw elevation (m)');ylabel('meander wavelength / channel width'); title('color is cpqual');box on
            saveas(gca, ...
            strcat('DebugParseTransect1BracketingCrosscheckD.eps'), ...
            'epsc');
        
%Corrected elevation:
figure;scatter(lzstrat(lambda>0)',2.*lambda(lambda>0)'./w(lambda>0)',100,log(s(11).cpstyle(lambda>0)')./log(2),'Filled')
grid on; colorbar; set(gca,'FontSize',18)
text(10+lzstrat(lambda>0)',2.*lambda(lambda>0)'./w(lambda>0)',num2str(acid(lambda>0)'))
%xlim([-2400 -2050])
xlabel('elevation relative to F1/F2 contact (m)');ylabel('meander wavelength / channel width'); title('color is log2(preservation style)');box on
            saveas(gca, ...
            strcat('DebugParseTransect1BracketingCrosscheckC_CorrectedElevation.eps'), ...
            'epsc');

figure;scatter(lzstrat(lambda>0)',2.*lambda(lambda>0)'./w(lambda>0)',100,s(11).cqual(lambda>0)','Filled')
grid on; colorbar; set(gca,'FontSize',18)
text(10+lzstrat(lambda>0)',2.*lambda(lambda>0)'./w(lambda>0)',num2str(acid(lambda>0)'))
%xlim([-2400 -2050])
xlabel('elevation relative to F1/F2 contact (m)');ylabel('meander wavelength (m) / width'); title('color is cpqual');box on
            saveas(gca, ...
            strcat('DebugParseTransect1BracketingCrosscheckD_CorrectedElevation.eps'), ...
            'epsc');
        
        

figure;scatter(s(11).cqual(lambda>0)',2.*lambda(lambda>0)'./w(lambda>0)',100,w(lambda>0)','Filled')
grid on; colorbar; set(gca,'FontSize',18)
text(0.2+s(11).cqual(lambda>0)',2.*lambda(lambda>0)'./w(lambda>0)',num2str(acid(lambda>0)'))
xlim([1 9])
xlabel('preservation quality');ylabel('meander wavelength / channel width'); title('color is width (m)');box on
            saveas(gca, ...
            strcat('DebugParseTransect1BracketingCrosscheckE.eps'), ...
            'epsc');

figure;scatter(s(11).cqual(lambda>0)',2.*lambda(lambda>0)'./w(lambda>0)',100,log(s(11).cpstyle(lambda>0))'./log(2),'Filled')
grid on; colorbar; set(gca,'FontSize',18)
text(0.2+s(11).cqual(lambda>0)',2.*lambda(lambda>0)'./w(lambda>0)',num2str(acid(lambda>0)'))
xlabel('preservation quality');ylabel('meander wavelength / channel width'); title('color is log2(preservation style)');box on
xlim([1 9])
saveas(gca, ...
            strcat('DebugParseTransect1BracketingCrosscheckF.eps'), ...
            'epsc');
