%parseAuxiliaryDataMarsRivers
%Edwin Kite

sbinsize = 10;% Stratigraphic bin size in m.
srange = [-300 300]; %Stratigraphic range to evaluate (relative to F1-F2 contact).
showregressionfigure = 1;
%***********************  0. Load auxiliary data:  ***********************
c = shaperead('channel_polygons_Transect_1_for_channel_deposit_proportions_points.shp');
    for l = find([c(:).Aeolis1z]==0);
        c(l).Aeolis1z = NaN;
    end
    for l = find([c(:).B20G02z]==0);
        c(l).B20G02z = NaN;
    end
mg = shaperead('Transect_1_Evidence_for_Channel_Migration_During_Deposition_Points.shp');

%Load topography and interpolated stratigraphic surfaces (for density calculations):

tCTX = imread('Transect_1_DTM_100m_Resolution_for_Density_Calculations_B20G02_stereo_nodata_removed.tif');
tHiRISE(1).t = imread('Transect_1_DTM_100m_Resolution_for_Density_Calculations_PSP_006683_1740_PSP_010322_1740.tif');
tHiRISE(2).t = imread('Transect_1_DTM_100m_Resolution_for_Density_Calculations_ESP_017548_1740_ESP_019104_1740.tif');
tHiRISE(3).t = imread('Transect_1_DTM_100m_Resolution_for_Density_Calculations_ESP_019038_1740_ESP_019882_1740.tif');


%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.
    tHiRISE(1).t(tHiRISE(1).t<-9999) = NaN;       
    overlap = find(tCTX~=0 & isnan(tCTX)==0 & tHiRISE(1).t~=0 & isnan(tHiRISE(1).t)==0);
    ctxonhirise = polyfit(tHiRISE(1).t(overlap),...
                            tCTX(overlap),1)
    ctxonhirise =  fliplr(ctxonhirise);

        if showregressionfigure == 1; 
            figure
            scatter(tHiRISE(1).t(overlap),tCTX(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
        end
    t = tCTX;

             t(isnan(t)==1 & tHiRISE(1).t < 0 & tHiRISE(1).t > -9999) ...
= tHiRISE(1).t(isnan(t)==1 & tHiRISE(1).t < 0 & tHiRISE(1).t > -9999)*ctxonhirise(2) + ctxonhirise(1);

%interpolated stratigraphic surfaces
s(1).t = imread('Transect_1_interpolation_planar_100m_Resolution_for_Density_Calculations.tif');
s(2).t = imread('Transect_1_interpolation_quadratic_100m_Resolution_for_Density_Calculations.tif');
s(3).t = imread('Transect_1_interpolation_IDW_100m_Resolution_for_Density_Calculations.tif');
s(4).t = imread('Transect_1_interpolation_ordinary_kriging_spherical_semivariogram_100m_Resolution_for_Density_Calculations.tif');

%Align topography by defining an appropriate meshgrid
te = load('Transect_1_interpolation_IDW_100m_Resolution_for_Density_Calculations.tfw');%topography extent
%Long axis in Transect 1 is y pos in m

[tx,ty] = meshgrid([te(5)+(te(1)/2):te(1):te(5)-(te(1)/2)+size(tCTX,2)*te(1)], ...
                   [te(6)+(te(4)/2):te(4):te(6)-(te(4)/2)+size(tCTX,1)*te(4)]); %Note 50m offsets

disp('Loading of auxiliary data complete')

%***********************  1a. Find absolute frequency of channel polygons vs. strat. elevation. ***********************

%'tesselate'
%Split each channel polygon into triangles, find the area of those triangles, and
%assign to each the mean elevation of the vertex points.

warning off
for l = 0:max([c(:).ORIG_FID])%loop over individual polygons. Note zero index
    cc = c([c(:).ORIG_FID] == l); %extract points from one polygon.
    checksum(l+1) = sum([cc(:).X]);
    if rem(l,100) == 0; disp(l); end
if length(cc)>=3 %At least one triangle?
    dtpolygon = delaunayTriangulation([[cc(:).X]',[cc(:).Y]'], ...
                [1:length([cc(:).X])'; 2:length([cc(:).X]') 1]'); %All edges are 'constrained' edges for DT 
    io = isInterior(dtpolygon);
    dt(l+1).dt = dtpolygon;
    %Important to prevent spurious values of elevation:
    
    hold on
    for mindex = 1:sum(io)
        okio = find(io==1);
        m = okio(mindex);
        pa(l+1,mindex)   = polyarea(dt(l+1).dt.Points(dt(l+1).dt.ConnectivityList(m,:),1), ...
                  dt(l+1).dt.Points(dt(l+1).dt.ConnectivityList(m,:),2));
        px(l+1,mindex)       = nanmean([cc(dt(l+1).dt.ConnectivityList(m,:)).X]);
        py(l+1,mindex)       = nanmean([cc(dt(l+1).dt.ConnectivityList(m,:)).Y]);
        line([cc(dt(l+1).dt.ConnectivityList(m,:)).X cc(dt(l+1).dt.ConnectivityList(m,1)).X], ...
             [cc(dt(l+1).dt.ConnectivityList(m,:)).Y cc(dt(l+1).dt.ConnectivityList(m,1)).Y],'Color',[.7 .7 .7])
        hold on
        
    end
       %For consistency, the pz is assigned to the elevation on the 100m
       %smoothed DTM.
       %Failure to do this means that channel areas that form ridges are
       %assigned to the wrong zstrat bin.
       pzctx(l+1,1:sum(io)) = interp2(tx,ty,tCTX,px(l+1,1:sum(io)),py(l+1,1:sum(io)));
    pzhirise(l+1,1:sum(io)) = interp2(tx,ty,tHiRISE(1).t,px(l+1,1:sum(io)),py(l+1,1:sum(io)));

       for iss = 1:4 %interpolated stratigraphic surface
           pziss(l+1,1:sum(io),iss) = interp2(tx,ty,s(iss).t,px(l+1,1:sum(io)),py(l+1,1:sum(io)));
       end
       
    
        line([cc(:).X],[cc(:).Y],'Color','k')
end
end
warning on
    title('Channel deposit debug figure')

disp('Debug Kludge!')
pz = pzctx;
      pz(isnan(pz)==1 & pzhirise<0 & pzhirise>-9999) = ...
pzhirise(isnan(pz)==1 & pzhirise<0 & pzhirise>-9999)*ctxonhirise(2) + ctxonhirise(1);

pa(pa<10) = NaN; %Remove polygon areas that are unimportantly small.
px(px==0) = NaN;
py(py==0) = NaN;
pz(pz==0) = NaN;

zsindex = 0;
%Aggregate/bin areas of channels.
%BIN 1: From srange(1) to [srange(1)+sbinsize].
whichstratzbin = zeros(size(px));
for zsbottom = [srange(1):sbinsize:srange(2)]; %bottom of statigraphic-elevation bin
    zsindex = zsindex + 1;
    for iss = 1:4; %interpolated stratigraphic surface
        isinstratzbin = ((pz(:,:) - squeeze(pziss(:,:,iss)) >  zsbottom & ...
                         (pz(:,:) - squeeze(pziss(:,:,iss)) < (zsbottom+sbinsize))));
        cdseh(zsindex,iss) = nansum(pa( [isinstratzbin] )); %channel-deposit stratigraphic elevation histogram
        whichstratzbin([isinstratzbin]) = zsindex;
    end
end

%***********************  1b. Find absolute frequency of evidence-for-channel-migration vs. strat. elevation. ***********************
%The same, but now for evidence-of-channel-migration polygons.
clear dt
figure
warning off
for l = 0:max([mg(:).ORIG_FID])%loop over individual polygons. Note zero index
    mgc = mg([mg(:).ORIG_FID] == l); %extract points from one polygon.
    if rem(l,100) == 0; disp(l); end
if length(mgc)>=3 %At least one triangle?
    dtpolygon = delaunayTriangulation([[mgc(:).X]',[mgc(:).Y]'], ...
                [1:length([mgc(:).X])'; 2:length([mgc(:).X]') 1]'); %All edges are 'constrained' edges for DT 
    io = isInterior(dtpolygon);
    dt(l+1).dt = dtpolygon;
    %Important to prevent spurious values of elevation:
    
    hold on
    for mindex = 1:sum(io)
        okio = find(io==1);
        m = okio(mindex);
        mga(l+1,mindex)   = polyarea(dt(l+1).dt.Points(dt(l+1).dt.ConnectivityList(m,:),1), ...
                   dt(l+1).dt.Points(dt(l+1).dt.ConnectivityList(m,:),2));
        mgx(l+1,mindex)       = nanmean([mgc(dt(l+1).dt.ConnectivityList(m,:)).X]);
        mgy(l+1,mindex)       = nanmean([mgc(dt(l+1).dt.ConnectivityList(m,:)).Y]);
        line([mgc(dt(l+1).dt.ConnectivityList(m,:)).X mgc(dt(l+1).dt.ConnectivityList(m,1)).X], ...
             [mgc(dt(l+1).dt.ConnectivityList(m,:)).Y mgc(dt(l+1).dt.ConnectivityList(m,1)).Y],'Color',[.7 .7 .7])
        hold on
        
    end
       %For consistency, the pz is assigned to the elevation on the 100m
       %smoothed DTM.
       %Failure to do this means that channel areas that form ridges are
       %assigned to the wrong zstrat bin.
       mgzctx(l+1,1:sum(io)) = interp2(tx,ty,tCTX,mgx(l+1,1:sum(io)),mgy(l+1,1:sum(io)));
    mgzhirise(l+1,1:sum(io)) = interp2(tx,ty,tHiRISE(1).t,mgx(l+1,1:sum(io)),mgy(l+1,1:sum(io)));

       for iss = 1:4 %interpolated stratigraphic surface
           mgziss(l+1,1:sum(io),iss) = interp2(tx,ty,s(iss).t,mgx(l+1,1:sum(io)),mgy(l+1,1:sum(io)));
       end
       
    
        line([mgc(:).X],[mgc(:).Y],'Color','k')
end
end
warning on
        title('Migration debug figure')
    
disp('Debug Kludge!')
mgz = mgzctx;
      mgz(isnan(mgz)==1 & mgzhirise<0 & mgzhirise>-9999) = ...
mgzhirise(isnan(mgz)==1 & mgzhirise<0 & mgzhirise>-9999)*ctxonhirise(2) + ctxonhirise(1);

mga(pa<10) = NaN; %Remove polygon areas that are unimportantly small.
mgx(px==0) = NaN;
mgy(py==0) = NaN;
mgz(pz==0) = NaN;

zsindex = 0;
%Aggregate/bin areas of channels.
%BIN 1: From srange(1) to [srange(1)+sbinsize].
mgwhichstratzbin = zeros(size(mgx));
for zsbottom = [srange(1):sbinsize:srange(2)]; %bottom of statigraphic-elevation bin
    zsindex = zsindex + 1;
    for iss = 1:4; %interpolated stratigraphic surface
        isinstratzbin = ((mgz(:,:) - squeeze(mgziss(:,:,iss)) >  zsbottom & ...
                         (mgz(:,:) - squeeze(mgziss(:,:,iss)) < (zsbottom+sbinsize))));
        mgseh(zsindex,iss) = nansum(mga( [isinstratzbin] )); %migration-evidence stratigraphic elevation histogram
        mgwhichstratzbin([isinstratzbin]) = zsindex;
    end
end


%***********************  2. Find normalized density of channel polygons (dividing by exposed area). ***********************
mask = NaN.*zeros(size(t)); %1 for area searched (i.e. area of HiRISE DTMs) and zero for areas not searched
for l = 1:length(tHiRISE)
    mask(tHiRISE(l).t < 0 & tHiRISE(l).t > -9999) = 1;
end


%DEBUG CHECK
figure;scatter(px(:),py(:),10,whichstratzbin(:))
hold on
scatter(px(whichstratzbin==26),py(whichstratzbin==26),10,whichstratzbin(whichstratzbin==26),'r')
hold on
contour(tx,ty,(mask.*((t - s(2).t))>-50).*(((t - s(2).t))<-40),[0.5 0.5])
contour(tx,ty,isnan(mask),[0.5 0.5])
%END OF DEBUG CHECK

figure;contourf(mask.*(t - s(2).t),[srange(1):sbinsize:srange(2)]);colorbar;set(gca,'ydir','reverse')
figure;contourf((mask.*((t - s(2).t))>-50).*(((t - s(2).t))<-40),[0.5 0.5]);colorbar;set(gca,'ydir','reverse')

%BIN 1: From srange(1) to [srange(1)+sbinsize].
eseh = histc(mask(:).*(t(:) - s(2).t(:)),[srange(1):sbinsize:srange(2)]); %exposure stratigraphic-elevation histogram
eseh = 100.*100.*eseh; %pixel size correction (want 100m x 100m bins)
figure;plot([(srange(1)+sbinsize/2):sbinsize:(srange(2)+sbinsize/2)],squeeze(cdseh(:,2))./eseh)
hold on
plot([(srange(1)+sbinsize/2):sbinsize:(srange(2)+sbinsize/2)],eseh./max(eseh),'g')
plot([(srange(1)+sbinsize/2):sbinsize:(srange(2)++sbinsize/2)],squeeze(cdseh(:,2))./max(squeeze(cdseh(:,2))),'r')
plot([(srange(1)+sbinsize/2):sbinsize:(srange(2)+sbinsize/2)],squeeze(mgseh(:,2))./eseh,'LineWidth',2)
plot([(srange(1)+sbinsize/2):sbinsize:(srange(2)+sbinsize/2)],squeeze(mgseh(:,2))./max(squeeze(mgseh(:,2))),'g','LineWidth',2)


(t - s(2).t); 


%***********************  3. Find channel-segment orientation versus stratigraphic elevation. ***********************

%***********************  4. Export summary data and display debug plots. ***********************