% Basins Test Code
% ----------------
% This worked well enough for what I needed to do.

% Uncomment to read the data (and create the S1-12 variables in the workspace)
%{
S1 = shaperead("hybas_eu_lev01_v1c.shp");
S2 = shaperead("hybas_eu_lev02_v1c.shp");
S3 = shaperead("hybas_eu_lev03_v1c.shp");
S4 = shaperead("hybas_eu_lev04_v1c.shp");
S5 = shaperead("hybas_eu_lev05_v1c.shp");
S6 = shaperead("hybas_eu_lev06_v1c.shp");
S7 = shaperead("hybas_eu_lev07_v1c.shp");
S8 = shaperead("hybas_eu_lev08_v1c.shp");
S9 = shaperead("hybas_eu_lev09_v1c.shp");
S10 = shaperead("hybas_eu_lev10_v1c.shp");
S11 = shaperead("hybas_eu_lev11_v1c.shp");
S12 = shaperead("hybas_eu_lev12_v1c.shp");
%}
%%return

% First figure just for basic visialisation
%figure(1)
%clf
% Examle basins
%mapshow(S2(3));
%mapshow(S5(326));
% Results S5 - [327,329,330,331,332]
%return
%geoshow(ax, S)
% X: [-8 0]
% Y: [54 60]
%xlim([-8 0])
%ylim([54 60])

% Vounding box of scotland for empirical results
scotland_bb=[-8, 54; 0, 60];

function bbox = convertBoundingBox(bb_in)
    % Converts the verted bounding box to the one required by bboxOverlapRatio
    bbox=[bb_in(1,1),bb_in(1,2),diff(bb_in(:,1)), diff(bb_in(:,2))];
end

% placeholder for the list of basin indexies
results=[];

% Adjsut this to look at different basins.
data=S12;

% convreting my crude coords to the correct format
% Scotland Bounding Box
sbb=convertBoundingBox(scotland_bb);

% looping over all the basins in that S level
for j = 1:size(data,1)

    % adjusting the bounding box format
    basin_bb=data(j).BoundingBox;
    bbb=convertBoundingBox(basin_bb);
    
    % calculating the overlap. Threshold should be meaningless, as the basins 
    % are << than scotland so should be "in" or "out" (so this will be ~1 or 0)
    % the "min" is essential.
    if (bboxOverlapRatio(bbb, sbb, "Min")>0.5)
        %{
        % for plotting and saving every approprate basin for manual overlap
        % Bad idea.
        clf
        mapshow(data(j));
        xlim([-8 0])
        ylim([54 60])
        txt = sprintf('%06d.png',j);
        saveas(gcf,txt)
        %}
        % saving the indecies of correct basins
        results=[results j];
        
        % Endo check (all were 0 for scotland!)
        if data(j).ENDO==1
            disp(j)
            disp("IS ENDO")
        end
    end
end

% Saving the big figure
% started at S9 and worked back to S5
% manually uncomment and change the data value above.
% not worht the effort of automating, only needed it once.
%
figure (2)
%clf
% Clf for only first one (need to overwrite after that)
hold on
for j = results
    %S5
    %mapshow(data(j), "FaceAlpha",0, "LineWidth",1, "EdgeColor","k") %
    %S6
    %mapshow(data(j), "FaceAlpha",0, "LineWidth",0.2, "EdgeColor","k") % 
    %S7
    %mapshow(data(j), "FaceAlpha",0, "LineWidth",0.2, "EdgeColor"," b") %
    %S8
    %mapshow(data(j), "FaceAlpha",0, "LineWidth",0.2, "EdgeColor","g" , "LineStyle",":")
    %S9
    mapshow(data(j), "FaceAlpha",0, "LineWidth",0.2, "EdgeColor","y" ,"LineStyle",":")
end

xlim([-8 0])
ylim([54 60])
pbaspect([1,1,1])
saveas(gcf,"combined9.png")
exportgraphics(gcf,'combinedHR9.png',"Resolution",1200)

% Notes of results.
% S4 is outline of scotland, so not really needed.
% Results S3 - [27]
% Results S4 - [115]
% Results S5 - [327,329,330,331,332]
% S6 - 1009
% S7 - 51 basins

%
