%loadWidthMeasurementsMarsRivers
%Edwin Kite
%Purpose: Imports ArcGIS picks into MATLAB for extraction of river-deposit
%dimensions

showregressionfigure = 1;
%******************************* LOAD WIDTH DATA **************************
if WhichTransect == 1
        wua = shaperead('Transect_1_Channel_Width_Measurements_Fourth_Pass_Points.shp');

        bankunit = [wua.UnitAssign];
        pstyle     = [wua.PStyle]';
        pqual      = [wua.PQuality]';  
        bid        = [wua.ORIG_FID]'+1;
        unitassign = [wua.UnitAssign]';
        tienumber  = [wua.TieNumber]';
        
        assert(min(diff(bid))==0,'check bank IDs')
        assert(max(diff(bid))==1,'check bank IDs')
        
        minx     = min([wua.X]); miny = min([wua.Y]);
        px       = ([wua.X]') - minx;
        py       = ([wua.Y]') - miny;
        origin   = [minx,miny];


        hiz =  [wua.Aeolis1z]'; %Note that the x038882z and x104548z elevations have note been checked and so shouldn't be used.
        ctxz = [wua.B20G02z]';
 

%Find z by (1) regressing HIZ on CTXZ for points 
%for which those overlap (2) applying the resulting correction to
%HIZ points where there is no CTXZ point.
    overlap = find(hiz~=0 & ctxz~=0);
    ctxonhirise = robustfit(hiz(overlap),ctxz(overlap)) 
        if showregressionfigure == 1; 
            figure
            scatter(hiz(overlap),ctxz(overlap));grid minor
            xlabel('hizstrat');ylabel('ctxzstrat');
            hold on
            line([-3000 0],[-3000*ctxonhirise(2)+ctxonhirise(1) ...
                                0*ctxonhirise(2)+ctxonhirise(1)]) %ctx value predicted using hirise value
            line([-3000 0],[-3000 0])
        end
            z = ctxz;
            z(z==0) = hiz(ctxz==0).*ctxonhirise(2) + ctxonhirise(1);
    
    residuals = ctxz(ctxz~=0 & hiz~=0) - ...
               (hiz(ctxz~=0  & hiz~=0)*ctxonhirise(2) + ctxonhirise(1));
    disp(strcat('std dev of CTX-HiRISE residuals (m): ',num2str(std(residuals))));
    
            z(0==(ctxz) & 0==(hiz)) = NaN;
        %%%%
        
        zc(:,1) = [wua.F1F2pl]; %contact elevation
        zc(:,2) = [wua.F1F2qd]; %contact elevation
        zc(:,3) = [wua.F1F2idw]; %contact elevation
        zc(:,4) = [wua.F1F2kr]; %contact elevation
        zc(:,5) = [wua.F1F2krv]; %contact elevation
        
        
        
elseif (WhichTransect == 2 || WhichTransect == 3)
       %wua = with unit assignments
       if WhichTransect == 2
        wua = shaperead('Transect_2_Channel_Width_Measurements_Fourth_Pass_Points.shp');
        bankunit   = [wua.UnitAssign]';
       elseif WhichTransect == 3
        wua = shaperead('Transect_3_Channel_Width_Measurements_Fourth_Pass_Points.shp');
        bankunit = [wua.UnitAssign]'.*0 + 3; %Force all banks to Unit F3.
       end
       
        pstyle     = [wua.PStyle]';
        pqual      = [wua.PQuality]';  
        bid        = [wua.ORIG_FID]'+1;
        unitassign = [wua.UnitAssign]';
        tienumber  = [wua.TieNumber]';
        
        assert(min(diff(bid))==0,'check bank IDs')
        assert(max(diff(bid))==1,'check bank IDs')
        
        minx     = min([wua.X]); miny = min([wua.Y]);
        px       = ([wua.X]') - minx;
        py       = ([wua.Y]') - miny;
        origin   = [minx,miny];
        
        whichdtm = 1+zeros(size(px)); %Transect 2 consists of a single HiRISE DTM
    if WhichTransect == 2
        hiz      = [wua.x474497z]';
    elseif WhichTransect == 3
        hiz      = [wua.x002279z]';
    end
    
        z        = hiz; %CTXz information is not available for these transects
        
    if WhichTransect == 2        
        %Assign stratigraphic elevations    
        zc(:,1) = [wua.F1F2pl]; %contact elevation
        zc(:,2) = [wua.F1F2qd]; %contact elevation
        zc(:,3) = [wua.F1F2idw]; %contact elevation
        zc(:,4) = [wua.F1F2kr]; %contact elevation
        zc(:,5) = [wua.F1F2krv]; %contact elevation
    elseif WhichTransect == 3
        disp('Implementing temporary kludge for Transect 3 stratigraphic elevations! (fixed -2km F1F2 elevation)')
        zc(:,1) = -2000 + zeros(size(wua));
        zc(:,2) = zc(:,1); zc(:,3) = zc(:,1); zc(:,4) = zc(:,1);
        zc(:,5) = zeros(size(wua)); %No stratigraphic interpolation
    end        
end

clear a overlap

disp('Exiting load_channel_widths')