%% Simple example show the Volumetric changes in wet area
clear
%close all
clc

savefigs=1;

%% control zones

% Cross shore 
B_offs_all = 1200;  % full domain
B_ons_all = 0;


%%  load TUD survey data

Mpath='F:\check_outs\tudelft\vlugtenburg\';
%Mpath='C:\Madeschipper\check_out\OpenEarthRawData\tudelft\vlugtenburg\scripts\Matthieu'

load([Mpath 'AllVbeachsurveys.mat'])
load([Mpath 'datesandpairs.mat'])


%% Load Jarkus data
load([Mpath 'scripts\Matthieu\AnalyseJarkus\JarkusVB.mat'])
disp('data loaded')



%% Fix Jetski data
disp('--- precheck')
% find spots where there is no data at or befor the x=0 point
for is=1:length(Vsurvey);
    disp(['survey  ' num2str(is)])
    for it=1:22
        if isnan(Vsurvey(is,1).transects.elevation2(301:311,it))
            disp(['error.. no data near -9 m NAP ' num2str(it)])
        end
        if isnan(Vsurvey(is,1).transects.elevation2(51:61,it))
            disp(['error.. no data near dune fence ' num2str(it)])
        end
        
    end
end

% Fill data with surveys from before/after
%%
Vsurvey(8,1).transects.elevation2(61,21)=Vsurvey(7,1).transects.elevation2(61,21);
Vsurvey(8,1).transects.elevation2(61,22)=Vsurvey(7,1).transects.elevation2(61,22);
%%
Vsurvey(1,1).transects.elevation2(301:311,1:2)=Vsurvey(2,1).transects.elevation2(301:311,1:2);
Vsurvey(1,1).transects.elevation2(301:311,1:2)=Vsurvey(2,1).transects.elevation2(301:311,1:2);
Vsurvey(4,1).transects.elevation2(301:311,1:4)=Vsurvey(3,1).transects.elevation2(301:311,1:4);
Vsurvey(4,1).transects.elevation2(301:311,1:4)=Vsurvey(3,1).transects.elevation2(301:311,1:4);
Vsurvey(9,1).transects.elevation2(51:61,22)=Vsurvey(8,1).transects.elevation2(51:61,22);
Vsurvey(10,1).transects.elevation2(301:311,1)=Vsurvey(9,1).transects.elevation2(301:311,1);
Vsurvey(11,1).transects.elevation2(51:61,22)=Vsurvey(10,1).transects.elevation2(51:61,22);
Vsurvey(12,1).transects.elevation2(51:61,22)=Vsurvey(13,1).transects.elevation2(51:61,22);
Vsurvey(14,1).transects.elevation2(51:61,22)=Vsurvey(13,1).transects.elevation2(51:61,22);
Vsurvey(15,1).transects.elevation2(301:311,1:2)=Vsurvey(14,1).transects.elevation2(301:311,1:2);
Vsurvey(19,1).transects.elevation2(301:311,1:2)=Vsurvey(20,1).transects.elevation2(301:311,1:2);
Vsurvey(18,1).transects.elevation2(301:311,1:8)=Vsurvey(19,1).transects.elevation2(301:311,1:8);
Vsurvey(17,1).transects.elevation2(301:311,1:5)=Vsurvey(18,1).transects.elevation2(301:311,1:5);
Vsurvey(22,1).transects.elevation2(51:61,22)=Vsurvey(21,1).transects.elevation2(51:61,22);
Vsurvey(25,1).transects.elevation2(51:61,22)=Vsurvey(24,1).transects.elevation2(51:61,22);
Vsurvey(28,1).transects.elevation2(51:61,22)=Vsurvey(27,1).transects.elevation2(51:61,22);
Vsurvey(31,1).transects.elevation2(51:61,22)=Vsurvey(30,1).transects.elevation2(51:61,22);

Vsurvey(35,1).transects.elevation2(51:61,22)=Vsurvey(34,1).transects.elevation2(51:61,22);

for is=1:length(Vsurvey);
    disp(['survey  ' num2str(is)])
    for it=1:22
        if isnan(Vsurvey(is,1).transects.elevation2(301:311,it))
            disp(['error.. no data near -9 m NAP ' num2str(it)])
        end
        if isnan(Vsurvey(is,1).transects.elevation2(51:61,it))
            disp(['error.. no data near dune fence ' num2str(it)])
        end
    end
end

%%%  No error should be left here in the workspace

% fix SoS (adjust 0.5 % by hand)
if 0
    Vsurvey(7,1).transects.elevation2(Vsurvey(7,1).transects.elevation2 <0)=Vsurvey(7,1).transects.elevation2(Vsurvey(7,1).transects.elevation2 <0)*0.995;
    Vsurvey(10,1).transects.elevation2(Vsurvey(10,1).transects.elevation2 <0)=Vsurvey(10,1).transects.elevation2(Vsurvey(10,1).transects.elevation2 <0)*1.005;
    Vsurvey(12,1).transects.elevation2(Vsurvey(12,1).transects.elevation2 <0)=Vsurvey(12,1).transects.elevation2(Vsurvey(12,1).transects.elevation2 <0)*0.995;
    Vsurvey(13,1).transects.elevation2(Vsurvey(13,1).transects.elevation2 <0)=Vsurvey(13,1).transects.elevation2(Vsurvey(13,1).transects.elevation2 <0)*1.005;
    Vsurvey(14,1).transects.elevation2(Vsurvey(14,1).transects.elevation2 <0)=Vsurvey(14,1).transects.elevation2(Vsurvey(14,1).transects.elevation2 <0)*0.995;
    Vsurvey(16,1).transects.elevation2(Vsurvey(16,1).transects.elevation2 <0)=Vsurvey(16,1).transects.elevation2(Vsurvey(16,1).transects.elevation2 <0)*1.005;
    Vsurvey(17,1).transects.elevation2(Vsurvey(17,1).transects.elevation2 <0)=Vsurvey(17,1).transects.elevation2(Vsurvey(17,1).transects.elevation2 <0)*0.995;
    Vsurvey(18,1).transects.elevation2(Vsurvey(18,1).transects.elevation2 <0)=Vsurvey(18,1).transects.elevation2(Vsurvey(18,1).transects.elevation2 <0)*0.995;
    Vsurvey(19,1).transects.elevation2(Vsurvey(19,1).transects.elevation2 <0)=Vsurvey(19,1).transects.elevation2(Vsurvey(19,1).transects.elevation2 <0)*0.995;
    % Vsurvey(29,1).transects.elevation2(Vsurvey(29,1).transects.elevation2 <0)=Vsurvey(29,1).transects.elevation2(Vsurvey(29,1).transects.elevation2 <0)*0.995;
     Vsurvey(34,1).transects.elevation2(Vsurvey(34,1).transects.elevation2 <0)=Vsurvey(34,1).transects.elevation2(Vsurvey(34,1).transects.elevation2 <0)*1.005;
      Vsurvey(35,1).transects.elevation2(Vsurvey(35,1).transects.elevation2 <0)=Vsurvey(35,1).transects.elevation2(Vsurvey(35,1).transects.elevation2 <0)*0.995;
       Vsurvey(36,1).transects.elevation2(Vsurvey(36,1).transects.elevation2 <0)=Vsurvey(36,1).transects.elevation2(Vsurvey(36,1).transects.elevation2 <0)*0.995;
filename=[ 'TotalVolume_corr']
else
    filename=[ 'TotalVolume_uncorr']
end


%% Jarkus method
% V vertical (planview boxes)
b_b=-10;
b_up=10;

dy=-1*diff(Vsurvey(1,1).transects.y); % alongshore distance between transects
for i_date=1:length(Vsurvey)
    for i_transect=1:22;
        xi=Vsurvey(i_date).transects.cr_sh;
        zi=Vsurvey(i_date).transects.elevation2(:,i_transect);
       if i_transect>=10
                b_up=4;  % cut off parts in south with vakantiehuisjes restaurants etc
       else
                b_up=10;
       end
        
        Vsurvey(i_date).VolV(i_transect)=jarkus_getVolume(xi, zi, b_up, b_b ,B_ons_all,B_offs_all);
        
        % Alongshore width represented by the transect (outside transects taken for half)
        if i_transect==1
            dy_tr=0.5*dy(1);
        end
        if i_transect==22
            dy_tr=0.5*dy(21);
        end
        if i_transect>=2 && i_transect<=21
            dy_tr=0.5*dy(i_transect-1)+0.5*dy(i_transect);
        end
        
        Vsurvey(i_date).transects.dy_tr(i_transect)=dy_tr; % store in Vsurvey
        Vsurvey(i_date).VolV_tr(i_transect)=Vsurvey(i_date).VolV(i_transect)*dy_tr; % Volume per transect.
        
    end
    Vsurvey(i_date).Vtotal=sum(Vsurvey(i_date).VolV_tr);  % Volume in survey area (from B_ons_all to B_off_all )
    Vtotal(i_date)=sum(Vsurvey(i_date).VolV_tr);% idem, but stored in a vector to use locally
end

%% Calculatre volumes Jarkus
for i_yr=1:length(Jarkus)
for i_transect=1:22%22;
xi=Vsurvey(i_date).transects.cr_sh;
zi=Jarkus(i_yr).transects.elevation(i_transect,:);
Jarkus(i_yr).VolV(i_transect)=jarkus_getVolume(xi, zi, b_up, b_b ,B_ons_all,B_offs_all);

if i_transect==1
    dy_tr=0.5*dy(1);
end
if i_transect==22
    dy_tr=0.5*dy(21);
end
if i_transect>=2 && i_transect<=21
    dy_tr=0.5*dy(i_transect-1)+0.5*dy(i_transect);
end

Jarkus(i_yr).transects.dy_tr(i_transect)=dy_tr;

Jarkus(i_yr).VolV_tr(i_transect)=Jarkus(i_yr).VolV(i_transect)*dy_tr;

end
Jarkus(i_yr).Vtotal=sum(Jarkus(i_yr).VolV_tr);
VJtotal(i_yr)=sum(Jarkus(i_yr).VolV_tr);
end


%% Figure Volume changes wrt to Jarkus 2007 
% Jarkus method
 figure
    set(gcf,'position',[ 189         372        1094         321])
set(gcf,'PaperPosition',[0 0  35 10])

windowSize = 3;
Vtotal_f=filtfilt(ones(1,windowSize)/windowSize,1,Vtotal)
plot(datenum(dates(pairs(:,1),:)),(Vtotal_f-VJtotal(1))/1E6,'r', 'linewidth' ,2)
hold on;
plot(datenum(dates(pairs(:,1),:)),(Vtotal-VJtotal(1))/1E6,'.-b')

%%% Plot Jarkus surveys
% plot(Jarkus(1,1).Mdate,VJtotal(1)-VJtotal(1),'*k')
% plot(Jarkus(2,1).Mdate,VJtotal(2)-VJtotal(1),'*b') Bad data only till -6
% plot(Jarkus(3,1).Mdate,VJtotal(3)-VJtotal(1),'*m') % During construction
% plot(Jarkus(4,1).Mdate,VJtotal(4)-VJtotal(1),'vk')


holdax=axis

x = [datenum(['01/01/2010'],'mm/dd/yyyy'); datenum(['12/31/2010'],'mm/dd/yyyy');datenum(['12/31/2010'],'mm/dd/yyyy') ; datenum(['01/01/2010'],'mm/dd/yyyy')];
y = [ 600; 600;0; 0];
tcolor(1,1,1:3) = 0.7*[1 1 1];
patch(x,y,tcolor,'facealpha',0.3,'EdgeColor','none')

x = [datenum(['01/01/2012'],'mm/dd/yyyy'); datenum(['12/31/2012'],'mm/dd/yyyy');datenum(['12/31/2012'],'mm/dd/yyyy') ; datenum(['01/01/2012'],'mm/dd/yyyy')];
y = [ 600; 600;0; 0];
tcolor(1,1,1:3) = 0.7*[1 1 1];
patch(x,y,tcolor,'facealpha',0.3,'EdgeColor','none')

axis(holdax)
datetick('x','mmm') 
%title('Sediment volume at Vlugtenburg beach')
grid on
xlim([datenum(['06/01/2009'],'mm/dd/yyyy') datenum(['12/31/2012'],'mm/dd/yyyy')])

text(datenum(['06/01/2010'],'mm/dd/yyyy'),2.42,'2010','color','w','fontsize',18,'FontWeight','bold')
text(datenum(['04/10/2012'],'mm/dd/yyyy'),2.42,'2012','color','w','fontsize',18,'FontWeight','bold')
text(datenum(['09/01/2009'],'mm/dd/yyyy'),2.42,'2009','color',0.7*[1 1 1],'fontsize',18,'FontWeight','bold')
text(datenum(['06/01/2011'],'mm/dd/yyyy'),2.42,'2011','color',0.7*[1 1 1],'fontsize',18,'FontWeight','bold')


xlabel('Date')
ylabel({'Additional volume in survey area';'wrt June 2007 [10^6 m^3] '})
%legend(' Monthly surveys','3 months moving average','location','northeast')


hold off

if savefigs==1

print ([ Mpath 'scripts\Matthieu\VolumeAnalysis\images\' filename '.jpg'],'-djpeg')
end


break
% 
% 
% %% plot may 2010 vs Jarkus
%  
% figure
% subplot (131)
% pcolor (Vsurvey(2,1).transects.y,Vsurvey(2,1).transects.cr_sh,Jarkus(4).transects.elevation2(:,:)')
%  shading flat
%  ylim([0 1200])
%  caxis([-9 5])
% 
%  subplot(132)
%  pcolor(Vsurvey(2,1).transects.y,Vsurvey(2,1).transects.cr_sh,Vsurvey(12).transects.elevation2)
%  shading flat
%  ylim([0 1200])
%   caxis([-9 5])
%   colorbar
%   
%  subplot(133)
%  pcolor(Vsurvey(2,1).transects.y,Vsurvey(2,1).transects.cr_sh,Jarkus(4,1).transects.elevation2.'-Vsurvey(12).transects.elevation2)
%  shading flat
%   ylim([0 1200])
%  caxis([-1 1])
%  colormap(jet)
%  colorbar
% 
% 
% 
% %% Jarkus method
%  figure
% 
%   plot(datenum(dates(pairs(1,1),:)),(muV(1)-muV(1)),'or')
%   hold on;
%   bar(datenum(dates(pairs(:,1),:)),DV_total)
% % 
%  plot(datenum(dates(pairs(:,1),:)),DV1_total,'*--g')
%  plot(datenum(dates(pairs(:,1),:)),DV2_total,'o--m')
%  plot(datenum(dates(pairs(:,1),:)),DV3_total,'+--k')
%   plot(datenum(dates(pairs(:,1),:)),DV4_total,'v--b')
%   plot(datenum(dates(pairs(:,1),:)),DV5_total,'^--r')
%   ylim([-250 200])
%  xlim([ datenum(dates(pairs(1,1),:)) datenum(dates(pairs(end,1),:)) ])
%  datetick('x') 
%  title('Jarkus method, Boxes')
% 
% 
%  
%  grid on
%  ylabel('\Delta V per m alongshore wrt June 2009')
%  legend('First survey','Total area','Area I','Area II','Area III','Area IV','Area V','location','northwest')
% 
%  set(gcf,'PaperPosition',[0  0  25 20])
% hold off
% 
% if savefigs==1
% filename=[ 'Jarkus Volume control area y=  ' num2str(B_south) '_' num2str(B_north) ', x=  ' num2str(B_ons_all) '_' num2str(B_offs_all)]
% print ([ Mpath '\VolumeAnalysis\images\' filename '.jpg'],'-djpeg')
% end
% 
% 
% 
% 
% 
% break
% 
% %% Jarkus method
% % V vertical (planview boxes)
% b_b=0.9;
% b_up=1.1;
% 
% for i_date=1:length(Vsurvey)
% for i_transect=1:22;
% xi=Vsurvey(i_date).transects.cr_sh;
% zi=Vsurvey(i_date).transects.elevation2(:,i_transect);
% VolMCL(i_date,i_transect)=jarkus_getVolume(xi, zi, b_up, b_b ,0,1200);
% X1m(i_date,i_transect)=VolMCL(i_date,i_transect)/0.2;
% end
% end
% 
% %% fig control area
% for i_date=2:length(Vsurvey)
% figure(46)
% pcolor(Vsurvey(i_date).transects.y(ind_voly),Vsurvey(i_date).transects.cr_sh(ind_volx),Temp(i_date).SedEroC)
% 
% colorbar; caxis([-0.15 .15]); ca=colormap(jet) ; ca(30:34,:)=1 ; colormap(ca); colorbar
% xlabel('Alongshore distance [m]'); ylabel('Cross shore distance [m]')
% end
% 
% 
% 
% 
% 
% 
% %% Figure first and last topo and areas of interest
% figure
% 
% subplot(3,2,[1 3])
% pcolor(xi,yi,altitude_i(:,:,1));  % plot first survey
% shading interp
% hold on;
% [C,h] = contour(xi,yi,altitude_i(:,:,1),[-9:1:5],'--k');
% clabel(C,h,'LabelSpacing',72,'FontSize',8,'Color','k','Rotation',0) ;
% 
% rectangle('Position', [B_south B_ons_all (B_north-B_south) (B_offs_all-B_ons_all)],'LineWidth',2,'LineStyle','--')
% rectangle('Position', [B_south B_ons_1+2 (B_north-B_south) (B_offs_1-B_ons_1)-2],'EdgeColor','y')
% rectangle('Position', [B_south B_ons_2 (B_north-B_south) (B_offs_2-B_ons_2)-5],'EdgeColor','r')
% rectangle('Position', [B_south B_ons_3 (B_north-B_south) (B_offs_3-B_ons_3)-5],'EdgeColor','m')
% rectangle('Position', [B_south B_ons_4 (B_north-B_south) (B_offs_4-B_ons_4)-5],'EdgeColor','b')
% rectangle('Position', [B_south B_ons_5 (B_north-B_south) (B_offs_5-B_ons_5)-5],'EdgeColor','g')
% text(B_south+10, B_offs_1-30,'Area I','FontAngle','italic','color','w')
% text(B_south+10, B_offs_2-30,'Area II','FontAngle','italic','color','w')
% text(B_south+10, B_offs_3-30,'Area III','FontAngle','italic','color','k')
% text(B_south+10, B_offs_4-30,'Area IV','FontAngle','italic','color','k')
% text(B_south+10, B_offs_5-30,'Area V','FontAngle','italic','color','w')
% axis equal
% ylim([-50 1300])
% caxis([-10 4])
% xlabel ('along shore [m]')
% ylabel ('cross shore [m]')
% title ([dates(pairs(1,1),[1:2 6:end]) ' survey' ])
% 
% subplot(3,2,[2 4 ])
% pcolor(xi,yi,altitude_i(:,:,end));  % plot latest survey
% shading interp
% hold on;
% [C,h] = contour(xi,yi,altitude_i(:,:,end),[-9:1:5],'--k');
% clabel(C,h,'LabelSpacing',72,'FontSize',8,'Color','k','Rotation',0) ;
% rectangle('Position', [B_south B_ons_all (B_north-B_south) (B_offs_all-B_ons_all)],'LineWidth',2,'LineStyle','--')
% rectangle('Position', [B_south B_ons_1+2 (B_north-B_south) (B_offs_1-B_ons_1)-2],'EdgeColor','y')
% rectangle('Position', [B_south B_ons_2 (B_north-B_south) (B_offs_2-B_ons_2)-5],'EdgeColor','r')
% rectangle('Position', [B_south B_ons_3 (B_north-B_south) (B_offs_3-B_ons_3)-5],'EdgeColor','m')
% rectangle('Position', [B_south B_ons_4 (B_north-B_south) (B_offs_4-B_ons_4)-5],'EdgeColor','b')
% rectangle('Position', [B_south B_ons_5 (B_north-B_south) (B_offs_5-B_ons_5)-5],'EdgeColor','g')
% text(B_south+10, B_offs_1-30,'Area I','FontAngle','italic','color','w')
% text(B_south+10, B_offs_2-30,'Area II','FontAngle','italic','color','w')
% text(B_south+10, B_offs_3-30,'Area III','FontAngle','italic','color','k')
% text(B_south+10, B_offs_4-30,'Area IV','FontAngle','italic','color','k')
% text(B_south+10, B_offs_5-30,'Area V','FontAngle','italic','color','w')
% axis equal
% ylim([-50 1300])
% caxis([-10 4])
% xlabel ('along shore [m]')
% title ([dates(pairs(end,1),[1:2 6:end]) ' survey' ])
% 
% subplot(3,2,[5 6])
% caxis([-10 4])
% colorbar('location', 'north')
% axis off
% 
% set(gcf,'PaperPosition',[0  0  25 20])
% print Overview_firstandlastsurvey.png -dpng
% 
% 
% 
% %% relate to wave power
%  load('waveheightanalysis/Pperiod.mat')
% 
%  %Pperiod(:,4)=Pperiod(:,1)./(Pperiod(:,3)-Pperiod(:,2));
%  
%  % timeserie
%  figure(34)
%  hold on;
% for ii=2:length(m_p)
%  bar(datenum(dates(pairs(ii,1),:)),(m_p(ii)-m_p(ii-1))*(B_offs_all-B_ons_all) )
% plot(datenum(dates(pairs(ii,1),:)),(m1_p(ii)-m1_p(ii-1))*(B_offs_1-B_ons_1),'*y')
%  plot(datenum(dates(pairs(ii,1),:)),(m2_p(ii)-m2_p(ii-1))*(B_offs_2-B_ons_2),'om')
%  plot(datenum(dates(pairs(ii,1),:)),(m3_p(ii)-m3_p(ii-1))*(B_offs_3-B_ons_3),'vr')
%  plot(datenum(dates(pairs(ii,1),:)),(m4_p(ii)-m4_p(ii-1))*(B_offs_4-B_ons_4),'^b')
%  plot(datenum(dates(pairs(ii,1),:)),(m5_p(ii)-m5_p(ii-1))*(B_offs_5-B_ons_5),'xg')
% end
% 
%  plot( Pperiod(2:end,3), Pperiod(2:end,1)/10^10)
% plot( Pperiod(2:end,3), Pperiod(2:end,4)/10^9)
% 
%  ylim([-120 120])
%  xlim([ datenum(dates(pairs(1,1),:)) datenum(dates(pairs(end,1),:)) ])
%  datetick('x') 
%  grid on
%  ylabel('\Delta V per m alongshore wrt June 2009')
%  legend('Total area','Area I','Area II','Area III','Area IV','Area V','location','northwest')
% 
%  set(gcf,'PaperPosition',[0  0  25 20])
% 
%  % Scatter
%   figure(35)
%  hold on;
% for ii=2:length(m_p)
% scatter(Pperiod(ii,4),(m_p(ii)-m_p(ii-1))*(B_offs_all-B_ons_all),'.k' )
% scatter(Pperiod(ii,4),(m1_p(ii)-m1_p(ii-1))*(B_offs_1-B_ons_1),'*y')
% scatter(Pperiod(ii,4),(m2_p(ii)-m2_p(ii-1))*(B_offs_2-B_ons_2),'om')
% scatter(Pperiod(ii,4),(m3_p(ii)-m3_p(ii-1))*(B_offs_3-B_ons_3),'vr')
% scatter(Pperiod(ii,4),(m4_p(ii)-m4_p(ii-1))*(B_offs_4-B_ons_4),'^b')
% scatter(Pperiod(ii,4),(m5_p(ii)-m5_p(ii-1))*(B_offs_5-B_ons_5),'xg')
% end
% 
% %  plot( Pperiod(2:end,3), Pperiod(2:end,1)/10^10)
% 
% 
% %  ylim([-120 120])
% %  xlim([ datenum(dates(pairs(1,1),:)) datenum(dates(pairs(end,1),:)) ])
% %  datetick('x') 
%  grid on
%  ylabel('\Delta V per m alongshore wrt June 2009')
%  legend('Total area','Area I','Area II','Area III','Area IV','Area V','location','northwest')
% 
%  set(gcf,'PaperPosition',[0  0  25 20])
% 
%  
%  
%  
% 
% %% Figure Volume changes in boxes
%  figure(32)
%  
% for pair_nr=2:19;
%  dV_all(pair_nr)=(m_p(pair_nr)-m_p(pair_nr-1))*(B_offs_all-B_ons_all);
%  dV_1(pair_nr)=(m1_p(pair_nr)-m1_p(pair_nr-1))*(B_offs_1-B_ons_1);
%  dV_2(pair_nr)=(m2_p(pair_nr)-m2_p(pair_nr-1))*(B_offs_2-B_ons_2);
%  dV_3(pair_nr)=(m3_p(pair_nr)-m3_p(pair_nr-1))*(B_offs_3-B_ons_3);
%  dV_4(pair_nr)=(m4_p(pair_nr)-m4_p(pair_nr-1))*(B_offs_4-B_ons_4);
%  dV_5(pair_nr)=(m5_p(pair_nr)-m5_p(pair_nr-1))*(B_offs_5-B_ons_5);
% end
% 
%  bar(datenum(dates(pairs(:,1),:)),dV_all)
% hold on;
% plot(datenum(dates(pairs(:,1),:)),dV_1,'*--y')
%  plot(datenum(dates(pairs(:,1),:)),dV_2,'m')
%  plot(datenum(dates(pairs(:,1),:)),dV_3,'r')
%  plot(datenum(dates(pairs(:,1),:)),dV_4,'b')
%  plot(datenum(dates(pairs(:,1),:)),dV_5,'g')
%  ylim([-250 200])
%  xlim([ datenum(dates(pairs(1,1),:)) datenum(dates(pairs(end,1),:)) ])
%  datetick('x') 
% 
%  grid on
%  ylabel('\Delta V per m alongshore wrt previous survey')
%  legend('Total area','Area I','Area II','Area III','Area IV','Area V','location','northwest')
% 
%  set(gcf,'PaperPosition',[0  0  25 20])
% print Vol_change_total.png -dpng
% 
% %% Figure mean elevation changes in boxes
%  figure
% 
%  bar(datenum(dates(pairs(:,1),:)),m_p-m_p(1))
% hold on;
% plot(datenum(dates(pairs(:,1),:)),m1_p-m1_p(1),'y')
%  plot(datenum(dates(pairs(:,1),:)),m2_p-m2_p(1),'m')
%  plot(datenum(dates(pairs(:,1),:)),m3_p-m3_p(1),'r')
%  plot(datenum(dates(pairs(:,1),:)),m4_p-m4_p(1),'b')
%  plot(datenum(dates(pairs(:,1),:)),m5_p-m5_p(1),'g')
%  ylim([-1.2 1])
%  xlim([ datenum(dates(pairs(1,1),:)) datenum(dates(pairs(end,1),:)) ])
%   datetick('x')
%  grid on
%  ylabel('mean elevation change')
%  legend('Total area','Area I','Area II','Area III','Area IV','Area V','location','northwest')
% 
%  set(gcf,'PaperPosition',[0  0  25 20])
% print Elevation_Areas.png -dpng
% 
% 
% %% figures all surveys
% 
% %% Figure first and last topo and areas of interest
% if 0
% figure(322)
% for ipair=1:length(pairs)
% pcolor(xi,yi,altitude_i(:,:,ipair));  % plot first survey
% shading interp
% hold on;
% [C,h] = contour(xi,yi,altitude_i(:,:,ipair),[-9:1:5],'--k');
% clabel(C,h,'LabelSpacing',72,'FontSize',8,'Color','k','Rotation',0) ;
% 
% rectangle('Position', [B_south B_ons_all (B_north-B_south) (B_offs_all-B_ons_all)],'LineWidth',2,'LineStyle','--')
% rectangle('Position', [B_south B_ons_1+2 (B_north-B_south) (B_offs_1-B_ons_1)-2],'EdgeColor','y')
% rectangle('Position', [B_south B_ons_2 (B_north-B_south) (B_offs_2-B_ons_2)-5],'EdgeColor','r')
% rectangle('Position', [B_south B_ons_3 (B_north-B_south) (B_offs_3-B_ons_3)-5],'EdgeColor','m')
% rectangle('Position', [B_south B_ons_4 (B_north-B_south) (B_offs_4-B_ons_4)-5],'EdgeColor','b')
% rectangle('Position', [B_south B_ons_5 (B_north-B_south) (B_offs_5-B_ons_5)-5],'EdgeColor','g')
% text(B_south+10, B_offs_1-30,'Area I','FontAngle','italic','color','w')
% text(B_south+10, B_offs_2-30,'Area II','FontAngle','italic','color','w')
% text(B_south+10, B_offs_3-30,'Area III','FontAngle','italic','color','k')
% text(B_south+10, B_offs_4-30,'Area IV','FontAngle','italic','color','k')
% text(B_south+10, B_offs_5-30,'Area V','FontAngle','italic','color','w')
% axis equal
% ylim([-50 1300])
% caxis([-10 4])
% xlabel ('along shore [m]')
% ylabel ('cross shore [m]')
% title ([dates(pairs(ipair,1),[4:6 1:3 7:end]) ' Survey' ])
% colorbar
% 
% set(gcf,'Renderer','zbuffer')
% set(gcf,'PaperPosition',[0  0  25 20])
% filenameComb=['CombinedSurvey' dates(pairs(ipair,1),[7:end 1:2 4:5 ]) '.png']
% print (filenameComb, '-dpng')
% hold off
% end
% end
% 

