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

%%% Options
savefigs=1; % 1 to save figures
addNorthMiddleSouth=0 % to add the partitions of north middle and south.
usecorrectedSurveys=1 % to use the corrected Volumes


%% control zones for volume
% Cross shore bpundaries
B_offs_all= 1200;
B_ons_all = 0;
% V vertical (planview boxes)
b_b=-10;
b_up=10;

%%  load survey data

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

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

load([Mpath '\AllVbeachsurveys.mat'])
if usecorrectedSurveys==1
    load([Mpath '\AllVbeachsurveys_corr.mat'])
end

disp('data loaded')

%% Jarkus method
if 1
    disp('Calculating Volumes survey:')
    %   dy=-1*diff(Vsurvey(1,1).transects.y);
    for i_date=1:length(Vsurvey)
        disp([' survey ' num2str(i_date) ])
        for i_transect=1:22;
            xi=Vsurvey(i_date).transects.cr_sh;
            zi=Vsurvey(i_date).transects.elevation2(:,i_transect);
            if i_transect>=10
                % pause
                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);
            if isnan(Vsurvey(i_date).VolV(i_transect))
                error('er gaat iets mis')
            end
        end
        % Remove transect 11 (vaak geen droog)
        Vsurvey(i_date).VolV(11)=NaN;
        
        % sommeer over alle transects, volum change per m alongshore
        Vsurvey(i_date).Vtotal=nansum(Vsurvey(i_date).VolV)/21;
        
        % idem, sommeer over alle transects   maar handiger array om verder mee te werken
        Vtotal(i_date)=nansum(Vsurvey(i_date).VolV)/21;
        Vtotal_std(i_date)=nanstd(Vsurvey(i_date).VolV)/21;
        
        if addNorthMiddleSouth==1
            % Sommeer volumes North, middle, South
            Vtotal_n(i_date)=nansum(Vsurvey(i_date).VolV(1:7))/7;
            Vtotal_m(i_date)=nansum(Vsurvey(i_date).VolV(8:15))/7; % 1 extra for the nan 11
            Vtotal_s(i_date)=nansum(Vsurvey(i_date).VolV(16:22))/7;
        end
    end
    
    %% Moving average filter
    windowSize = 3;
    Vtotal_f=filtfilt(ones(1,windowSize)/windowSize,1,(Vtotal-Vtotal(1)));

    %% remove nr 4 survey.. heeft Bas verklooid
    Vtotal(4)=NaN; Vtotal_f(4)=NaN;
    
    %% North middle south filtered
     if addNorthMiddleSouth==1
        Vtotal_n_f=filtfilt(ones(1,windowSize)/windowSize,1,(Vtotal_n-Vtotal_n(1)));
        Vtotal_m_f=filtfilt(ones(1,windowSize)/windowSize,1,(Vtotal_m-Vtotal_m(1)));
        Vtotal_s_f=filtfilt(ones(1,windowSize)/windowSize,1,(Vtotal_s-Vtotal_s(1)));
    Vtotal_n_f(4)=NaN; Vtotal_m_f(4)=NaN; Vtotal_s_f(4)=NaN
    end
       
    %% Bulk volume losses
    disp('Bulk volume losses:')
    disp([ 'after app 1 yr:' num2str(Vtotal_f(13))])
    disp([ 'after app 2 yrs:' num2str(Vtotal_f(24)) ', volume loss 2nd yr:' num2str(Vtotal_f(24)-Vtotal_f(13)) ])
    disp([ 'after app 3 yrs:' num2str(Vtotal_f(35)) ', volume loss 3nd yr:' num2str(Vtotal_f(35)-Vtotal_f(24)) ])
    
    %% Figure Volume changes wrt to first survey in boxes
    % Planfom boxes!
    figure
    
    plot(datenum(dates(pairs(:,1),:)),(Vtotal-Vtotal(1)),'.-k')
    hold on;
    %plot(datenum(dates(pairs(:,1),:)),(Vtotal_std),'color',0.5*[1 1 1])
    plot(datenum(dates(pairs(:,1),:)),(Vtotal_f),'color','r')
    plot(datenum(dates(pairs(13,1),:)),(Vtotal_f(13)),'or')
    plot(datenum(dates(pairs(24,1),:)),(Vtotal_f(24)),'or')
    plot(datenum(dates(pairs(35,1),:)),(Vtotal_f(35)),'or')
    
    if addNorthMiddleSouth==1
        %%% Plot Partitions
        plot(datenum(dates(pairs(:,1),:)),(Vtotal_n_f),'color',[0 0 0.5])
        plot(datenum(dates(pairs(:,1),:)),(Vtotal_m_f),'color',[0 0.3 0.3])
        plot(datenum(dates(pairs(:,1),:)),(Vtotal_s_f),'color',[0 0.5 0 ])
        text(datenum(dates(pairs(end,1),:))-60,(Vtotal_n_f(end))-5,'North','color',[0 0 0.5])
        text(datenum(dates(pairs(end,1),:))-70,(Vtotal_m_f(end))-5,'Middle','color',[0 0.3 0.3])
        text(datenum(dates(pairs(end,1),:))-70,(Vtotal_s_f(end))-5,'South','color',[0 0.5 0])
        text(datenum(dates(pairs(end,1),:))-60,(Vtotal_f(end))-5,'Total','color','r')
    end
    
    % Check which surveys are close to T0+n years
    % plot(datenum(dates(pairs(1,1),:))+365,0,'+r') % zero point
    % plot(datenum(dates(pairs(13,1),:)),0,'+k') % zero point
    %
    % plot(datenum(dates(pairs(1,1),:))+2*365,0,'+r') % zero point
    % plot(datenum(dates(pairs(24,1),:)),0,'+k') % zero point
    %
    % plot(datenum(dates(pairs(1,1),:))+3*365,0,'+r') % zero point
    % plot(datenum(dates(pairs(35,1),:)),0,'+k') % zero point
    
    ylim([-250 100])
    xlim([ datenum(dates(pairs(1,1),:)) datenum(dates(pairs(end,1),:)) ])
    datetick('x')
    grid on
    ylabel({'\Delta V per m alongshore wrt July 2009';'Erosion          Accretion' })
    
    plot(datenum(dates(pairs(1,1),:)),0,'or') % zero point
    
    
    set(gcf,'position',[ 189         372        1094         321])
    set(gcf,'PaperPosition',[0 0  35 10]*0.7)
    hold off
    
    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;-300; -300];
    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;;-300; -300];
    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(['04/01/2009'],'mm/dd/yyyy') datenum(['12/31/2012'],'mm/dd/yyyy')])
    
    text(datenum(['06/01/2010'],'mm/dd/yyyy'),-220,'2010','color','w','fontsize',18,'FontWeight','bold')
    text(datenum(['04/10/2012'],'mm/dd/yyyy'),-220,'2012','color','w','fontsize',18,'FontWeight','bold')
    text(datenum(['08/01/2009'],'mm/dd/yyyy'),-220,'2009','color',0.7*[1 1 1],'fontsize',18,'FontWeight','bold')
    text(datenum(['06/01/2011'],'mm/dd/yyyy'),-220,'2011','color',0.7*[1 1 1],'fontsize',18,'FontWeight','bold')
    
    
    xlabel('Date')
 
    
    
    if savefigs==1
        filename=[ 'VolumeAreaPerM']
        print ([ Mpath 'scripts\Matthieu\VolumeAnalysis\images\' filename '.jpg'],'-djpeg')
    end
      
end


%% 
NetCDFfilename=[Mpath '\Vlugtenburg.nc'];
line=7 % number of transect line(s) 2 be extracted (value(s) between 1 and 22)
i_date=22

survey_ind = list_surveys;

[cross_shore,altitude] = combine_surveys_nr(survey_ind,i_date,line)

% x=nc_varget(NetCDFfilename,'cross_shore',[0],[-1]);
% alt=nc_varget(NetCDFfilename,'altitude',[date-1 line-1 0],[1 1 -1]);   % index -1 because netcdf starts counting at zero.. Weird Yanks.

%%

figure
set(gcf,'position',[ 428   388   660   305])
    set(gcf,'PaperPosition',[0 0  25 10]*0.7)
hold all 
    plot(cross_shore,altitude,'color',1/256*[133 87 35],'linewidth',2)
    plot([135 1480],[0 0],'b')
    plot([-300 135],[0 0],':b')
    plot([0 0],[4.35 6],'k','linewidth',2)
    plot([0 0 0],[5 5.5 6],'xk')

    indp=[61:301];
    pro=altitude;
    xp=[0 0 cross_shore(indp)' cross_shore(indp(end))'];
    yp=[-10 4.5 pro(indp)' -10];
    tcolor(1,1,1:3) = 0.6*[1 1 1];
    patch(xp,yp,tcolor,'facealpha',0.4,'EdgeColor','k')

    grid on

set(gca,'XDir','reverse');
xlim([-300 1500])
ylim([-10.5 12])
% 
%    set(gca,'xtick',[0:100:1200])
%     set(gca,'ytick',[-10:2:5])
%     
xlabel('Cross-shore distance [m]')
ylabel('Profile elevation wrt NAP [m]')
  if savefigs==1
        filename=[ 'VolumeCalcMethod']
        print ([ Mpath 'scripts\Matthieu\VolumeAnalysis\images\' filename '.jpg'],'-djpeg')
    end




%%
break
%% Jarkus method
% V horizontal (discs )

b_b=-10;
b_up=10;

b_bA=2.5;
b_upA=10;

b_bB =-1.5;
b_upB=2.5;

b_bC=-5;
b_upC=-1.5;

b_bD=-7;
b_upD=-5;

b_bE=-10;
b_upE=-7;


for i_date=1:length(Vsurvey)
    disp(['' num2str(i_date) ])
    for i_transect=[1:10 12: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
            b_upA=4;
        end
        Vsurvey(i_date).VolV(i_transect)=jarkus_getVolume(xi, zi, b_up, b_b ,B_ons_all,B_offs_all);
        
        Vsurvey(i_date).VolVA(i_transect)=jarkus_getVolume(xi, zi, b_upA , b_bA ,B_ons_all,B_offs_all);
        Vsurvey(i_date).VolVB(i_transect)=jarkus_getVolume(xi, zi, b_upB, b_bB ,B_ons_all,B_offs_all);
        Vsurvey(i_date).VolVC(i_transect)=jarkus_getVolume(xi, zi, b_upC, b_bC ,B_ons_all,B_offs_all);
        Vsurvey(i_date).VolVD(i_transect)=jarkus_getVolume(xi, zi, b_upD, b_bD ,B_ons_all,B_offs_all);
        Vsurvey(i_date).VolVE(i_transect)=jarkus_getVolume(xi, zi, b_upE, b_bE  ,B_ons_all,B_offs_all);
    end
    
    
    % Remove transect 11 (vaak geen droog)
    Vsurvey(i_date).VolV(11)=NaN; Vsurvey(i_date).VolVA(11)=NaN; Vsurvey(i_date).VolVB(11)=NaN; Vsurvey(i_date).VolVC(11)=NaN; Vsurvey(i_date).VolVD(11)=NaN; Vsurvey(i_date).VolVE(11)=NaN;
    
    % sommeer over alle transects, volum change per m alongshore
    Vsurvey(i_date).Vtotal=nansum(Vsurvey(i_date).VolV)/21;
    Vsurvey(i_date).VtotalA=nansum(Vsurvey(i_date).VolVA)/21;
    Vsurvey(i_date).VtotalB=nansum(Vsurvey(i_date).VolVB)/21;
    Vsurvey(i_date).VtotalC=nansum(Vsurvey(i_date).VolVC)/21;
    Vsurvey(i_date).VtotalD=nansum(Vsurvey(i_date).VolVD)/21;
    Vsurvey(i_date).VtotalE=nansum(Vsurvey(i_date).VolVE)/21;
    
    % idem, sommeer over alle transects   maar handiger array om verder mee te werken
    Vtotal(i_date)=nansum(Vsurvey(i_date).VolV)/21;
    VtotalA(i_date)=nansum(Vsurvey(i_date).VolVA)/21;
    VtotalB(i_date)=nansum(Vsurvey(i_date).VolVB)/21;
    VtotalC(i_date)=nansum(Vsurvey(i_date).VolVC)/21;
    VtotalD(i_date)=nansum(Vsurvey(i_date).VolVD)/21;
    VtotalE(i_date)=nansum(Vsurvey(i_date).VolVE)/21;
end
%% Check is volumes add up
% for i_date=1:length(Vsurvey)
%
%     difference=Vtotal(i_date) - ( VtotalA(i_date)+VtotalB(i_date)+VtotalC(i_date)+VtotalD(i_date)+VtotalE(i_date));
%     if difference >=1
%         i_date
%         disp('Volumes do not add up!')
%     end
% end

Vtotal(4)=NaN;
VtotalB(4)=NaN;
VtotalC(4)=NaN;



%% Figure Volume changes wrt to first survey in boxes
% elevation2 boxes!
figure
bar(datenum(dates(pairs(:,1),:)),(Vtotal-Vtotal(1)),'facecolor',0.5*[1 1 1])
hold on;
plot(datenum(dates(pairs(:,1),:)),(VtotalA-VtotalA(1)),'*--g')
plot(datenum(dates(pairs(:,1),:)),(VtotalB-VtotalB(1)),'o--m')
plot(datenum(dates(pairs(:,1),:)),(VtotalC-VtotalC(1)),'+--k')
plot(datenum(dates(pairs(:,1),:)),(VtotalD-VtotalD(1)),'v--b')
plot(datenum(dates(pairs(:,1),:)),(VtotalE-VtotalE(1)),'^--r')
windowSize = 3;
Va_f=filtfilt(ones(1,windowSize)/windowSize,1,(VtotalA-VtotalA(1)))
Vb_f=filtfilt(ones(1,windowSize)/windowSize,1,(VtotalB-VtotalB(1)))
Vc_f=filtfilt(ones(1,windowSize)/windowSize,1,(VtotalC-VtotalC(1)))
Vd_f=filtfilt(ones(1,windowSize)/windowSize,1,(VtotalD-VtotalD(1)))
Ve_f=filtfilt(ones(1,windowSize)/windowSize,1,(VtotalE-VtotalE(1)))
plot(datenum(dates(pairs(:,1),:)),(Va_f),'g')
plot(datenum(dates(pairs(:,1),:)),(Vb_f),'m')
plot(datenum(dates(pairs(:,1),:)),(Vc_f),'k')
plot(datenum(dates(pairs(:,1),:)),(Vd_f),'b')
plot(datenum(dates(pairs(:,1),:)),(Ve_f),'r')
plot([datenum(dates(pairs(1,1),:)) datenum(dates(pairs(1,1),:))],[-100 100],'--k') % zero point
plot([datenum(dates(pairs(1,1),:)) datenum(dates(pairs(1,1),:))]+365,[-100 100],'--k') % zero point
plot([datenum(dates(pairs(1,1),:)) datenum(dates(pairs(1,1),:))]+365*2,[-100 100],'--k') % zero point
plot([datenum(dates(pairs(1,1),:)) datenum(dates(pairs(1,1),:))]+365*3,[-100 100],'--k') % zero point

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 June 2009';'Erosion          Accretion' })
legend('Total profile','Area A (z_b> 2.5 m)','Area B (2.5>z_b>-1.5 m)','Area C (-1.5>z_b>-5 m)','Area D (-5>z_b>-7 m)','Area E (z_b< -7 m)','location','northwest')
plot(datenum(dates(pairs(1,1),:)),0,'or') % zero point
set(gcf,'PaperPosition',[0  0  18 15])
hold off

title('Cumulative changes since June 2009. Horizontal slices')

if savefigs==1
    filename=[ 'VolumezonesHorizontal']
    print ([ 'F:\check_outs\tudelft\vlugtenburg\scripts\Matthieu\VolumeAnalysis\images\' filename '.jpg'],'-djpeg')
end


if 0
    %% Show profiles boxes
    figure
    plot(Vsurvey(1).transects.cr_sh,nanmean(Vsurvey(1).transects.elevation2(:,:)'),'r','linewidth',2)
    hold all
    plot(Vsurvey(1).transects.cr_sh,nanmean(Vsurvey(length(Vsurvey)).transects.elevation2(:,:)'),'k','linewidth',2)
    
    legend('June 2009','Latest survey','location','NorthWest')
    plot([205 1200],[0 0],'b')
    plot([0 205],[0 0],':b')
    xlim([0 1200])
    grid on
    
    set(gca,'xtick',[0:100:1200])
    set(gca,'ytick',[-10:2:5])
    
    
    
    % set(gca,'xtick',[B_ons_all B_offs_1 B_offs_2 B_offs_3 B_offs_4 B_offs_5])
    % set(gca,'ytick',[b_b b_upE b_upD b_upC b_upB b_upA])
    
    xlabel ('Cross shore distance [m]')
    ylabel ('elevation2 [m NAP]')
    
    set(gcf,'position',[ 560   662   557   286])
    
    set(gca,'xdir','reverse')
    
    
    %% Show profiles slices
    figure
    plot(Vsurvey(1).transects.cr_sh,nanmean(Vsurvey(length(Vsurvey)).transects.elevation2(:,:)'),'k','linewidth',2)
    hold all
    %legend('June 2009','Nov 2011')
    plot([135 1200],[0 0],'b')
    plot([0 135],[0 0],':b')
    xlim([0 1200])
    grid on
    
    indp=[61:72];
    pro=nanmean(Vsurvey(length(Vsurvey)).transects.elevation2(:,:)');
    xp=[0 0 Vsurvey(1).transects.cr_sh(indp)' Vsurvey(1).transects.cr_sh(indp(end))'];
    yp=[2.5 10 pro(indp) 2.5];
    zdata=ones(size(xp))
    patch(xp,yp,zdata,'g','FaceAlpha',0.85)
    
    indp=[73:104];
    pro=nanmean(Vsurvey(length(Vsurvey)).transects.elevation2(:,:)');
    xp=[0 0 55 Vsurvey(1).transects.cr_sh(indp)' 222];
    yp=[-1.5 2.5 2.5 pro(indp) -1.5];
    zdata=ones(size(xp))
    patch(xp,yp,zdata,'m','FaceAlpha',0.85)
    
    indp=[106:165];
    pro=nanmean(Vsurvey(length(Vsurvey)).transects.elevation2(:,:)');
    xp=[0 0 Vsurvey(1).transects.cr_sh(indp)' 530];
    yp=[-5 -1.5 pro(indp) -5];
    zdata=ones(size(xp))
    patch(xp,yp,zdata,'k','FaceAlpha',0.85)
    
    indp=[167:199];
    pro=nanmean(Vsurvey(length(Vsurvey)).transects.elevation2(:,:)');
    xp=[0 0 Vsurvey(1).transects.cr_sh(indp)' 690];
    yp=[-7 -5 pro(indp) -7];
    zdata=ones(size(xp))
    patch(xp,yp,zdata,'b','FaceAlpha',0.85)
    
    indp=[200:301];
    pro=nanmean(Vsurvey(length(Vsurvey)).transects.elevation2(:,:)');
    xp=[0 0 Vsurvey(1).transects.cr_sh(indp)' 1200];
    yp=[-10 -7 pro(indp) -10];
    zdata=ones(size(xp))
    patch(xp,yp,zdata,'r','FaceAlpha',0.85)
    
    ylim([-9 +5])
    grid on
    
    % set(gca,'xtick',[B_ons_all B_offs_1 B_offs_2  B_offs_3 B_offs_4 B_offs_5])
    % set(gca,'ytick',[b_b b_upE b_upD b_upC b_upB b_upA])
    
    xlabel ('Cross shore distance [m]')
    ylabel ('elevation [m NAP]')
    
    set(gca,'xtick',[0:100:1200])
    set(gca,'ytick',[-10:2:5])
    
    
    set(gca,'xdir','reverse')
    
    set(gcf,'position',[ 560   808   394   140])
end
%% Show profiles boxes
figure
plot(Vsurvey(1).transects.cr_sh,nanmean(Vsurvey(length(Vsurvey)).transects.elevation2(:,:)'),'k','linewidth',2)
hold all
%legend('June 2009','Nov 2011')
plot([135 1200],[0 0],'b')
plot([0 135],[0 0],':b')
xlim([0 1200])
grid on

indp=[61:81];
pro=nanmean(Vsurvey(length(Vsurvey)).transects.elevation2(:,:)');
xp=[0 0 Vsurvey(1).transects.cr_sh(indp)' 100];
yp=[-10 10 pro(indp) -10];
zdata=ones(size(xp))
patch(xp,yp,zdata,'g','FaceAlpha',0.85)

indp=[81:121];
pro=nanmean(Vsurvey(length(Vsurvey)).transects.elevation2(:,:)');
xp=[100 100 Vsurvey(1).transects.cr_sh(indp)' 300];
yp=[-10 1.55 pro(indp) -10];
zdata=ones(size(xp))
patch(xp,yp,zdata,'m','FaceAlpha',0.85)

indp=[121:161];
pro=nanmean(Vsurvey(length(Vsurvey)).transects.elevation2(:,:)');
xp=[300 300 Vsurvey(1).transects.cr_sh(indp)' 500];
yp=[-10 -3.6 pro(indp) -10];
zdata=ones(size(xp))
patch(xp,yp,zdata,'k','FaceAlpha',0.85)

indp=[161:221];
pro=nanmean(Vsurvey(length(Vsurvey)).transects.elevation2(:,:)');
xp=[500 500 Vsurvey(1).transects.cr_sh(indp)' 800];
yp=[-10 -5.2 pro(indp) -10];
zdata=ones(size(xp))
patch(xp,yp,zdata,'b','FaceAlpha',0.85)

indp=[221:301];
pro=nanmean(Vsurvey(length(Vsurvey)).transects.elevation2(:,:)');
xp=[800 800 Vsurvey(1).transects.cr_sh(indp)' 1200];
yp=[-10 -7.7 pro(indp) -10];
zdata=ones(size(xp))
patch(xp,yp,zdata,'r','FaceAlpha',0.85)

ylim([-9 +5])
grid on

% set(gca,'xtick',[B_ons_all B_offs_1 B_offs_2  B_offs_3 B_offs_4 B_offs_5])
% set(gca,'ytick',[b_b b_upE b_upD b_upC b_upB b_upA])

xlabel ('Cross shore distance [m]')
ylabel ('elevation2 [m NAP]')

set(gcf,'position',[ 560   808   394   140])

%plot(Vsurvey(1).transects.cr_sh,nanmean(Vsurvey(1).transects.elevation2(:,:)'),'r','linewidth',2)
%% Figure first and last topo and areas of interest
figure
% alongshore boundaries
B_north  = 1745;
B_south  = 0;

subplot(3,2,[1 3])
i_date=1;
pcolor(Vsurvey(i_date).transects.y,Vsurvey(i_date).transects.cr_sh,Vsurvey(i_date).transects.elevation2);  % plot first survey
shading interp
hold on;
[C,h] = contour(Vsurvey(i_date).transects.y,Vsurvey(i_date).transects.cr_sh,Vsurvey(i_date).transects.elevation2,[-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','g')
rectangle('Position', [B_south B_ons_2 (B_north-B_south) (B_offs_2-B_ons_2)-5],'EdgeColor','m')
rectangle('Position', [B_south B_ons_3 (B_north-B_south) (B_offs_3-B_ons_3)-5],'EdgeColor','k')
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','r')
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([-10 1210])
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 ])
i_date=length(Vsurvey);
pcolor(Vsurvey(i_date).transects.y,Vsurvey(i_date).transects.cr_sh,Vsurvey(i_date).transects.elevation2);  % plot latest survey
shading interp
hold on;
[C,h] = contour(Vsurvey(i_date).transects.y,Vsurvey(i_date).transects.cr_sh,Vsurvey(i_date).transects.elevation2,[-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_firstandlastsurveyPlanformboxes.png -dpng




break


%% relate to wave power
load('waveheightanalysis/Pperiod.mat')

%Pperiod(:,4)=Pperiod(:,1)./(Pperiod(:,3)-Pperiod(:,2));

% timeserie
figure(34)
hold on;
plot(datenum(dates(pairs(:,1),:)),(VtotalB-VtotalB(1))/1745,'o--m')
% 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=6:min(length(Vsurvey),length(Pperiod))
    % 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),(VtotalB(ii)-VtotalB(ii-1))/1745 ,'om')
    scatter(Pperiod(ii,4),(VtotalC(ii)-VtotalC(ii-1))/1745 ,'+k')
    % 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 elevation2 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 elevation2 change')
legend('Total area','Area I','Area II','Area III','Area IV','Area V','location','northwest')

set(gcf,'PaperPosition',[0  0  25 20])
print elevation2_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


