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

savefigs=0;

%% control zones
% alongshore boundaries
B_north  = 1500;
B_south  = 200;
% Cross shore 
B_offs_all = 1200;
B_ons_all = 000;
B_offs_1   = 1200;
B_ons_1   = 800;
B_offs_2   = 800;
B_ons_2   = 500;
B_offs_3   = 500;
B_ons_3   = 300;
B_offs_4   = 300;
B_ons_4   = 100;
B_offs_5   = 100;
B_ons_5   = 0;

%%  load survey data

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

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

disp('data loaded')
%% Determine index with boxes
i_date=1; % just for determining the boxes


ind_voly=find(Vsurvey(i_date).transects.y  >= B_south & Vsurvey(i_date).transects.y  <= B_north );  % Alonsgore boundaries

ind_volx=find(Vsurvey(i_date).transects.cr_sh  >= B_ons_all & Vsurvey(i_date).transects.cr_sh  <= B_offs_all ); % total
ind_volx1=find(Vsurvey(i_date).transects.cr_sh >= B_ons_1 & Vsurvey(i_date).transects.cr_sh < B_offs_1);
ind_volx2=find(Vsurvey(i_date).transects.cr_sh >= B_ons_2 & Vsurvey(i_date).transects.cr_sh < B_offs_2);
ind_volx3=find(Vsurvey(i_date).transects.cr_sh >= B_ons_3 & Vsurvey(i_date).transects.cr_sh < B_offs_3);
ind_volx4=find(Vsurvey(i_date).transects.cr_sh >= B_ons_4 & Vsurvey(i_date).transects.cr_sh < B_offs_4);
ind_volx5=find(Vsurvey(i_date).transects.cr_sh >= B_ons_5 & Vsurvey(i_date).transects.cr_sh < B_offs_5);

Dx=Vsurvey(i_date).transects.cr_sh(ind_volx(end))-Vsurvey(i_date).transects.cr_sh(ind_volx(1));
Dx1=Vsurvey(i_date).transects.cr_sh(ind_volx1(end))-Vsurvey(i_date).transects.cr_sh(ind_volx1(1));
Dx2=Vsurvey(i_date).transects.cr_sh(ind_volx2(end))-Vsurvey(i_date).transects.cr_sh(ind_volx2(1));
Dx3=Vsurvey(i_date).transects.cr_sh(ind_volx3(end))-Vsurvey(i_date).transects.cr_sh(ind_volx3(1));
Dx4=Vsurvey(i_date).transects.cr_sh(ind_volx4(end))-Vsurvey(i_date).transects.cr_sh(ind_volx4(1));
Dx5=Vsurvey(i_date).transects.cr_sh(ind_volx5(end))-Vsurvey(i_date).transects.cr_sh(ind_volx5(1));



%% Correction
if 0
    %Control area
    ind_volyC=find(Vsurvey(i_date).transects.y  >= 200 & Vsurvey(i_date).transects.y  <= 1500 );  % Alonsgore boundaries
    ind_volxC=find(Vsurvey(i_date).transects.cr_sh  >= 900 & Vsurvey(i_date).transects.cr_sh  <= 1100 ); % total
    
    for i_date=1:length(Vsurvey)
        surveynew=Vsurvey(i_date).transects.elevation(ind_volxC,ind_volyC);
        surveyref=Vsurvey(28).transects.elevation(ind_volxC,ind_volyC);
        CorMat=NaN(size(surveyref));
        CorMat=surveynew./surveyref;
        Cor=nanmean(nanmean(CorMat));
        if 0; figure(55) ; pcolor(CorMat);shading interp; colorbar; caxis([0.9 1.1]);  pause;end
        Cor=nanmean(nanmean(CorMat));
        %sedero=surveynew-surveyref;
        
        [datajetxi,datajetxj]=find(Vsurvey(i_date).transects.elevation <= -1);
        
        %sedero(datajet)=Cor*surveynew(datajet)-surveyref(datajet);
        Vsurvey(i_date).transects.elevation(datajetxi,datajetxj)=Vsurvey(i_date).transects.elevation(datajetxi,datajetxj)./Cor;
        
    end
    disp('Correcting surveys!!')
    
end


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

for i_date=1:length(Vsurvey)
for i_transect=1:22;
xi=Vsurvey(i_date).transects.cr_sh;
zi=Vsurvey(i_date).transects.elevation(:,i_transect);
VolMCL(i_date,i_transect)=jarkus_getVolume(xi, zi, b_up, b_b ,0,1200);
end
end
%%
VolMCL_first=repmat(VolMCL(1,:),length(Vsurvey),1);
DeltaMCL=VolMCL-VolMCL_first;
dMCL=[zeros(1,22); diff(VolMCL)];

%%
figure
subplot(211)
plot(datenum(dates(pairs(:,1),:)),DeltaMCL(:,2:20))
 datetick('x','keeplimits') 
grid on

subplot(212)
plot(datenum(dates(pairs(:,1),:)),dMCL(:,2:20))
 datetick('x','keeplimits') 
grid on

%%
%% Figure Volume changes wrt to first survey in boxes 
% MA method
 figure(32)

  plot(datenum(dates(pairs(1,1),:)),(muV(1)-muV(1)),'or')
  hold on;
  bar(datenum(dates(pairs(:,1),:)),muV-muV(1))

 plot(datenum(dates(pairs(:,1),:)),muV1-muV1(1),'*--g')
 plot(datenum(dates(pairs(:,1),:)),muV2-muV2(1),'o--m')
 plot(datenum(dates(pairs(:,1),:)),muV3-muV3(1),'+--k')
  plot(datenum(dates(pairs(:,1),:)),muV4-muV4(1),'v--b')
  plot(datenum(dates(pairs(:,1),:)),muV5-muV5(1),'^--r')
  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')
 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=[ '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

%% Jarkus method
 figure(33)

  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') 


 
 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=[ '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

%% 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


