%parse_Transect_1_Channel_Width_Measurements_Second_Pass_Points_for_LPSC_2014
%This is preliminary: the right way to do this is to dot points at equal
%distances along both banks (rather than just measuring the width at
%vertices).


csvread Transect_1_Channel_Width_Measurements_Second_Pass_Points_csv.csv

a = ans;

ctxz = a(:,18);
hiz = max(a(:,17),a(:,16));

z = min(ctxz,hiz);
z(z<-99999) = NaN;

dz = z - a(:,15); %minus kriged contact

dz(z<-99999) = NaN;

%check image:
scatter(a(:,12),a(:,13),10,dz);

%take the square root of the semivariance:
rootsemi = sqrt(a(:,14));

%check image:
figure
scatter(a(:,12),a(:,13),10,rootsemi);

%find channel widths (pointwise)

bid = a(:,11); %bank ID

whichdtm = (a(:,16)>-99999)+1;%which DTM

dcls_temp = NaN;
for l = 2:(size(a,1)-1) %omit endpoints of banks
    l = l
    if bid(l+1) == bid(l) && bid(l-1) == bid(l)
    %find distance and coordinates of closest line segment
        P = [a(l,12),a(l,13),0];
    for m = 2:(size(a,1)-1) 
        if bid(m+1) == bid(m) && l~=m && bid(m)~=bid(l) %DO consider the first line segment (for which bid(m-1)~=bid(m))
    %row vectors!
            Q1 = [a(m,12),a(m,13),0];
            Q2 = [a(m+1,12),a(m+1,13),0];
            dcls_temp_v(m,:) = norm(cross(Q2-Q1,P-Q1))/norm(Q2-Q1);
            dcls_temp_line(m) = max(dcls_temp_v(m,:));
    %segment is finite
            dQ1(m) = sqrt((Q1(1) - P(1))^2 + (Q1(2) - P(2))^2);
            dQ2(m) = sqrt((Q2(1) - P(1))^2 + (Q2(2) - P(2))^2);
            dcls_temp(m) = median([dQ1(m),dQ2(m),dcls_temp_line(m)]);
        else
            dcls_temp(m) = NaN;%skip
        end
    end
    
    [C,I] = nanmin(dcls_temp);
    dcls(l) = C;
    cls(l) = I;
    clsbank(l) = bid(I);
    
     %We want to suppress the case where the closest point is a
            %vertex at the END of a line segment
            if dQ1(I) == dcls_temp(I) && bid(I-1)~=bid(I) 
                dcls(l) = NaN; cls(l) = NaN; clsbank(l) = NaN;
            elseif dQ2(I) == dcls_temp(I) && bid(I+2)~=bid(I)
                dcls(l) = NaN; cls(l) = NaN; clsbank(l) = NaN;
            end
    
    else
        cls(l) = NaN;
        dcls(l) = NaN;
        clsbank(l) = NaN;
        %skip
    end
    
end

dcls(size(dcls,2):size(a,1)) = NaN;

%now find median width, dz, and rootsemi for each channel
for nn = 1:(max(bid)/2) %each channel has two banks; they are assumed to be in order of this pairing!
    cw(nn)   = nanmedian(dcls(bid==(nn*2)-2 | bid==((nn*2)-1)));
    css(nn)  = sum(bid==((nn*2) - 2) | bid==((nn*2)-1));%channel sample size (for calculating pooled variance, later)
    cwdev(nn) = nanstd(dcls(bid==(nn*2)-2 | bid==(nn*2)-1));
    cdz(nn) = nanmedian(dz((bid==(nn*2)-2 | bid==(nn*2)-1) & isnan(dcls')==0));
    cqual(nn) = nanmedian(a((bid==(nn*2)-2 | bid==(nn*2)-1),6)); %note that FID == a(:,2)
    cz(nn) = nanmedian(z((bid==(nn*2)-2 | bid==(nn*2)-1) & isnan(dcls')==0));
    crootsemi(nn) = nanmedian(rootsemi((bid==(nn*2)-2 | bid==(nn*2)-1) & isnan(dcls')==0));
    cwhichdtm(nn) = nanmean(whichdtm(bid==(nn*2)-2 | bid==(nn*2)-1));
    cpstyle(nn) = nanmean(a((bid==(nn*2)-2 | bid==(nn*2)-1),5));
end

%Kludge fix for mixed preservation styles
cpstyle(cpstyle==9) = 8;
cpstyle(cpstyle==10) = 8;


figure

    errorbar(cdz,cw,cwdev,'xk','LineWidth',2)
    hold on
    scatter(cdz,cw,50,'co','fill')
    scatter(cdz(cdz<0),cw(cdz<0),50,'bo','fill')
    scatter(cdz(cdz>0 & cdz<50),cw(cdz>0 & cdz<50),50,'wo','fill')
    scatter(cdz(cdz>0 & cdz<50),cw(cdz>0 & cdz<50),50,'ko')

    %errorbar(cdz(cwhichdtm==2),cw(cwhichdtm==2),cwdev(cwhichdtm==2),'xr')

    %herrorbar(cdz,cw,crootsemi,'x')
xlabel('Relative Time (Height Above F1/F2 Contact (m))','FontSize',20)
ylabel('River Channel Width (m)','FontSize',20)
set(gca,'FontSize',20)
set(gca,'xminortick','on')
set(gca,'yminortick','on')
ylim([0 70])
xlim([-100 150])
%grid on
box on
line([-100,0],[29.3,29.3],'Color','b','LineWidth',2)
line([50,150],[15.8616,15.8616],'Color','c','LineWidth',2)
line([0,50],[ 26.8001, 26.8001],'Color','k','LineWidth',2)


%Quick summation of probability distribution (assuming uncorrelated, Guassian error
%bars!) From Wikipedia, "Multivariate normal distribution"
p = zeros(601,401);
 for nn = 1:(max(bid)/2) %add the probability kernel of each observation
oldp = p;
nn
            sigmax = crootsemi(nn);
            sigmay = cwdev(nn);
            normfactor = 1/(2*pi*sigmax*sigmay);
            mux = cdz(nn);
            muy = cw(nn);
            if sum(isnan([sigmax,sigmay,mux,muy]))==0
                for x = -300:1:300 %height above F1/F2 contact (m)
                    for y = -100:1:300 %river channel width (m)

                            p(x+301,y+101) = p(x+301,y+101) +  normfactor*exp(-1/2*(((x-mux)^2/(sigmax^2))+((y-muy)^2/(sigmay^2))));
                    end
                end
            end
                %For debug only
                %check that the additional probability sums to 1 (note that
                %the discretisation and the cutoffs will introduce errors!)
                %if abs((sum(p(:)) - sum(oldp(:)))-1) > 0.01
                %    error('MVN Dist Checksum Error!')
                %end
 end

 [X,Y] = meshgrid(-300:1:300,-100:1:300);
%contour(X,Y,p');grid on %show the sum of all kernels

%show 1 sigma and 2 sigma error windows for channel width as a function of
%stratigraphic position

pcs = cumsum(p,2);
ps = sum(p');
for l = 1:601 %cumsum normalized
pcsn(l,:) = pcs(l,:)./ps(l);
end

hold on
contour(X,Y,pcsn',[0.05,(0.31731/2),1-((0.31731/2)),0.95]);

%find pooled variance for 50m stratigaphic-position bins: (Note to self:
%should preferably use *fractional* variance?)
stratbin = 1;
for stratz = -200:50:200
    to_use(stratbin).v = find(cdz>stratz & cdz <stratz+50);
    ninbin(stratbin) = length(to_use(stratbin).v);
    pooledvar(stratbin) = sqrt(sum(css(to_use(stratbin).v).*(cwdev(to_use(stratbin).v).^2))/sum(css(to_use(stratbin).v)));
    stratbin = stratbin + 1;
end

%Overall pooled variance:
opooledvar = sqrt(sum(css.*(cwdev.^2))/sum(css))


