%parseAuxiliaryDataMarsRivers
%Edwin Kite

%Purpose:

if exist('WhichTransect','var') == 0; WhichTransect = 3; end

sbinsize = 10;% Stratigraphic bin size in m.
srange = [-300 300]; %Stratigraphic range to evaluate (relative to F1-F2 contact).
showregressionfigure = 1;
isschoice = 2; %For display purposes.
%***********************  0. Load auxiliary data:  ***********************
if WhichTransect == 1
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');
    %Not currently used
    %cl = shaperead('Transect_1_Channel_Segments_for_Channel_Alignment_Points.shp');
elseif WhichTransect == 2
    c = shaperead('channel_polygons_PSP_007474_1745_second_pass_for_channel_deposit_proportions_points.shp');
   mg = shaperead('Transect_2_Evidence_for_Channel_Migration_During_Deposition_Points.shp');
elseif WhichTransect == 3
    c = shaperead('channel_polygons_area_of_PSP_002002_1735_for_channel_deposit_proportions_points.shp');
   mg = []; %No evidence for channel migration during deposition for this transect.
end


%Load topography and interpolated stratigraphic surfaces (for density calculations):
if WhichTransect == 1 
    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);

elseif WhichTransect == 2
    t = imread('Transect_2_DTM_100m_Resolution_for_Density_Calculations.tif');
elseif WhichTransect == 3
    t = double(imread('Transect_3_DTM_100m_Resolution_for_Density_Calculations.tif'));
    t(t<-3000) = NaN; %Specific to Transect 3
end
%End of load topography and interpolated stratigraphic surfaces (for density calculations):

if (WhichTransect == 1 || WhichTransect == 2)
    %interpolated stratigraphic surfaces
    s(1).t = imread(strcat('Transect_',num2str(WhichTransect),'_interpolation_planar_100m_Resolution_for_Density_Calculations.tif'));
    s(2).t = imread(strcat('Transect_',num2str(WhichTransect),'_interpolation_quadratic_100m_Resolution_for_Density_Calculations.tif'));
    s(3).t = imread(strcat('Transect_',num2str(WhichTransect),'_interpolation_IDW_100m_Resolution_for_Density_Calculations.tif'));
    s(4).t = imread(strcat('Transect_',num2str(WhichTransect),'_interpolation_ordinary_kriging_spherical_semivariogram_100m_Resolution_for_Density_Calculations.tif'));
elseif (WhichTransect == 3)
    disp('Temporary kludge for stratigraphic surface elevation for Transect 3!')
    s(1).t = double(t.*0 + -2000);
    s(2).t = s(1).t; s(3).t = s(1).t; s(4).t = s(1).t;
end

%Align topography by defining an appropriate meshgrid
copyfile('Transect_1_DTM_100m_Resolution_for_Density_Calculations_B20G02_stereo_nodata_removed.tfw', ...
         'Transect_1_DTM_100m_Resolution_for_Density_Calculations.tfw')
te = load(strcat('Transect_',num2str(WhichTransect),'_DTM_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(t,2)*te(1)], ...
                   [te(6)+(te(4)/2):te(4):te(6)-(te(4)/2)+size(t,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.
       if (WhichTransect == 1 || WhichTransect == 2 || WhichTransect == 3)
       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
       end
       
       if WhichTransect == 1   
           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)));
       else %Only one pz option
              pz(l+1,1:sum(io)) = interp2(tx,ty,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

if WhichTransect == 1 
    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);
end

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

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

%***********************  1b. Find absolute frequency of evidence-for-channel-migration vs. strat. elevation. ***********************
%The same, but now for evidence-of-channel-migration polygons.
if isempty(mg) == 0
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.
   
       if (WhichTransect == 1 || WhichTransect == 2 || WhichTransect == 3)
       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
       end
       
       if WhichTransect == 1   
           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)));
       else %Only one pz option
              mgz(l+1,1:sum(io)) = interp2(tx,ty,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
   
if WhichTransect == 1
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);
end

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


if (WhichTransect == 1 || WhichTransect == 2 || WhichTransect == 3)
clear isinstratzbin mgwhichstratzbin
    %Aggregate/bin areas of channels.
%BIN 1: From srange(1) to [srange(1)+sbinsize].
for iss = 1:4; %interpolated stratigraphic surface
    zsindex = 0;
    inbin = zeros(size(mgx));
for zsbottom = [srange(1):sbinsize:srange(2)]; %bottom of statigraphic-elevation bin
    zsindex = zsindex + 1;
        isinstratzbin = find((mgz(:,:) - squeeze(mgziss(:,:,iss)) >  zsbottom & ...
                         (mgz(:,:) - squeeze(mgziss(:,:,iss)) < (zsbottom+sbinsize))));
        mgseh(zsindex,iss) = nansum(mga( [isinstratzbin] )); %channel-deposit stratigraphic elevation histogram
        if isempty(isinstratzbin) == 0;
            inbin([isinstratzbin]) = zsindex;
        end
    end
            mgwhichstratzbin(iss,:,:) = inbin;
end
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
if WhichTransect == 1
for l = 1:length(tHiRISE)
    mask(tHiRISE(l).t < 0 & tHiRISE(l).t > -9999) = 1;
end
else
    mask(t>-9999 & t<0) = 1;
end

%DEBUG CHECK
mystratbin = squeeze(whichstratzbin(isschoice,:,:));
figure;scatter(px(:),py(:),10,mystratbin(:))
hold on
scatter(px(mystratbin==27),py(mystratbin==27),10, ...
    mystratbin(mystratbin==27),'k')
hold on
contour(tx,ty,(mask.*((t - s(isschoice).t))>-40).*(((t - s(isschoice).t))<-30),[0.5 0.5])
contour(tx,ty,isnan(mask),[0.5 0.5])
    title('Channel deposit debug figure')
    
if exist('mgwhichstratzbin') ~=0
mystratbin = squeeze(mgwhichstratzbin(isschoice,:,:));
figure;scatter(mgx(:),mgy(:),10,mystratbin(:))
hold on
scatter(mgx(mystratbin==27),mgy(mystratbin==27),10, ...
    mystratbin(mystratbin==27),'k')
hold on
contour(tx,ty,(mask.*((t - s(isschoice).t))>-40).*(((t - s(isschoice).t))<-30),[0.5 0.5])
contour(tx,ty,isnan(mask),[0.5 0.5])
         title('Migration debug figure')
end
%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(isschoice).t))>-50).*(((t - s(isschoice).t))<-40),[0.5 0.5]);colorbar;set(gca,'ydir','reverse')

%BIN 1: From srange(1) to [srange(1)+sbinsize].
for iss = 1:4
eseh(:,iss) = histc(mask(:).*(t(:) - s(iss).t(:)),[srange(1):sbinsize:srange(2)]); %exposure stratigraphic-elevation histogram
end
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(:,isschoice))./eseh(:,isschoice))
hold on
plot([(srange(1)+sbinsize/2):sbinsize:(srange(2)+sbinsize/2)],eseh(:,isschoice)./max(eseh(:,isschoice)),'g')
plot([(srange(1)+sbinsize/2):sbinsize:(srange(2)++sbinsize/2)],squeeze(cdseh(:,isschoice))./max(squeeze(cdseh(:,isschoice))),'r')
if exist('mgseh') ~=0
    plot([(srange(1)+sbinsize/2):sbinsize:(srange(2)+sbinsize/2)],squeeze(mgseh(:,isschoice))./eseh(:,isschoice),'LineWidth',2)
    plot([(srange(1)+sbinsize/2):sbinsize:(srange(2)+sbinsize/2)],squeeze(mgseh(:,isschoice))./max(squeeze(mgseh(:,isschoice))),'g','LineWidth',2)
end


%***********************  3. Find channel-segment orientation versus stratigraphic elevation. ***********************
%TO BE IMPLEMENTED LATER

%***********************  4. Export summary data and display debug plots. ***********************
        
AuxiliaryDataSaveName = strcat('outputAuxiliaryDataMarsRivers_Transect_',num2str(WhichTransect))
if exist('mgseh') ~=0
    save(AuxiliaryDataSaveName,'mask','s','t','eseh','cdseh','sbinsize','srange','mgseh')
else
    save(AuxiliaryDataSaveName,'mask','s','t','eseh','cdseh','sbinsize','srange')
end