clear all
close all
clc

grid_res = 2.5;
load('MAT\rtk_data');
load('MAT\sb25m_data');
load('D:\werf_je\Werk\1.201.819-Proefsuppletie_Galgeplaat_2e_evaluatie\Data_analyse\Cirkel_bestek\cb_data');
titlenum = {'a','b','c','d','a','b','c','d','a','b','c','d','a','b','c','d'};

%Gecombineerde set maken
len1 = length(rtk_data.datum);
len2 = length(sb25m_data.datum);
%
rtk_sb25m_data.datum(1:len1) = rtk_data.datum;
rtk_sb25m_data.datum(len1+1:len1+len2) = sb25m_data.datum;
%
rtk_sb25m_data.naam(1:len1) = rtk_data.naam;
rtk_sb25m_data.naam(len1+1:len1+len2) = sb25m_data.naam;
%
rtk_sb25m_data.x(1:len1) = rtk_data.x;
rtk_sb25m_data.x(len1+1:len1+len2) = sb25m_data.x;
%
rtk_sb25m_data.y(1:len1) = rtk_data.y;
rtk_sb25m_data.y(len1+1:len1+len2) = sb25m_data.y;
%
rtk_sb25m_data.z(1:len1) = rtk_data.z;
rtk_sb25m_data.z(len1+1:len1+len2) = sb25m_data.z;
%
rtk_sb25m_data.type(1:len1) = {'RTK'};
rtk_sb25m_data.type(len1+1) = {'singlebeam 25 m'};
rtk_sb25m_data.type(len1+2:len1+4) = {'gecombineerde meting'};
rtk_sb25m_data.type(len1+5:len1+len2) = {'singlebeam 25 m'};

%Even sorten op datum
[hm ind] = sort([rtk_sb25m_data.datum{:}]);
rtk_sb25m_data.datum = rtk_sb25m_data.datum(ind);
rtk_sb25m_data.naam = rtk_sb25m_data.naam(ind);
rtk_sb25m_data.x = rtk_sb25m_data.x(ind);
rtk_sb25m_data.y = rtk_sb25m_data.y(ind);
rtk_sb25m_data.z = rtk_sb25m_data.z(ind);
rtk_sb25m_data.type= rtk_sb25m_data.type(ind);

for i = 1:length(rtk_sb25m_data.naam)
    clear X Y IN
    figure
    [C,h] = contourf(rtk_sb25m_data.x{i},rtk_sb25m_data.y{i},rtk_sb25m_data.z{i});
    hold on
    axis equal
    set(h,'LevelList',[-1.5:0.25:0.75])
    axis([0.548 0.564 3.972 3.988]*1e5)
    plot(cb_data.x,cb_data.y,'k-','LineWidth',1.0)
    xlabel('xRD (m)')
    ylabel('yRD (m)') 
    %dum = rem(i,4+max(1-rem(i,4),0));
    title(['(',titlenum{i},') bodemligging [m NAP], ',datestr(rtk_sb25m_data.datum{i}),', ',rtk_sb25m_data.type{i}])
    caxis([-1.5 0.75])
    hcb = colorbar;
    set(hcb,'ytick',[-1.5:0.25:0.75])
    saveas(gcf,['PNG\rtk_sb25m_',num2str(i,'%02.0f'),'.png'])
    saveas(gcf,['EPS\rtk_sb25m_',num2str(i,'%02.0f'),'.eps'],'psc2')
    %Kubering
    [X,Y] = meshgrid(rtk_sb25m_data.x{i},rtk_sb25m_data.y{i});
    [IN] = inpolygon(X,Y,cb_data.x,cb_data.y);
    vol(i) = nansum(nansum(IN.*rtk_sb25m_data.z{i}))*grid_res^2; 
end
combi_vol(:,1) = [rtk_sb25m_data.datum{1:4}];
combi_vol(:,2) = vol(1:4)-vol(1);
save MAT\combi_vol combi_vol
%%%
ind1 = [1,5:7,9,11:12];
sb25m_vol(:,1) = [rtk_sb25m_data.datum{ind1}];
sb25m_vol(:,2) = vol(ind1)-vol(1);
save MAT\sb25m_vol sb25m_vol
%%%
ind2 = [1,8,10,13];
rtk_vol(:,1) = [rtk_sb25m_data.datum{ind2}];
rtk_vol(:,2) = vol(ind2)-vol(1);
save MAT\rtk_vol rtk_vol
%%%
