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

%% control zones
% alongshore boundaries
B_north  = 1500;
B_south  = 200;
% Cross shore 
B_west_all   = 1200;
B_east_all   = 0;
B_west_1   = 1200;
B_east_1  = 800;
B_west_2   = 800;
B_east_2  = 500;
B_west_3   = 500;
B_east_3  = 300;
B_west_4   = 300;
B_east_4  = 100;
B_west_5   = 100;
B_east_5  = 0;

%%  load survey data
% locate vlugtenburg directory on this computer
curdir=pwd;
ind_v=findstr(curdir, 'vlugtenburg');
pathvlugtenburg=[(curdir(1:ind_v+10)) '\'];  % gives back cometing like F:\check_outs\tudelft\vlugtenburg\
NetCDFfilename=[pathvlugtenburg 'vlugtenburg.nc'];
%nc_dump(NetCDFfilename);  % give a quick overview of the contents

dates=nc_varget(NetCDFfilename,'time');     % survey dates
survey_ind=list_surveys;                    % list with the surveys and the survey pairs
indjet=survey_ind.indjet; indwalk=survey_ind.indwalk; pairs=survey_ind.Survey_pairs;  % rewrite to local variables... bit easier later on

% Data interpolated on surveylines
cross_shore=nc_varget(NetCDFfilename,'cross_shore');
dist_alongshore = nc_varget(NetCDFfilename,'dist_alongshore');
% 'raw' data on path
path=nc_varget(NetCDFfilename,'survey_path');


%% make rectified top view of the surveys.
[xi,yi] = meshgrid(dist_alongshore,cross_shore);

ind_vol=find(xi >= B_south & xi <= B_north & yi >= B_east_all & yi < B_west_all);
ind_vol1=find(xi >= B_south & xi <= B_north & yi >= B_east_1 & yi < B_west_1);
ind_vol2=find(xi >= B_south & xi <= B_north & yi >= B_east_2 & yi < B_west_2);
ind_vol3=find(xi >= B_south & xi <= B_north & yi >= B_east_3 & yi < B_west_3);
ind_vol4=find(xi >= B_south & xi <= B_north & yi >= B_east_4 & yi < B_west_4);
ind_vol5=find(xi >= B_south & xi <= B_north & yi >= B_east_5 & yi < B_west_5);


for pair_nr=1:length(pairs);
[cross_shore_i,altitude_i(:,:,pair_nr)] = combine_surveys_nr(survey_ind,pair_nr);
sedero=altitude_i(:,:,pair_nr)-altitude_i(:,:,1);

% calculate Volumes above -10 and mean elevation changes per area
Vol_p_10(pair_nr)=sum(altitude_i(ind_vol)+10); %Volume area above 10 m
m_p(pair_nr)=nanmean(sedero(ind_vol));
Vol1_p_10(pair_nr)=sum(altitude_i(ind_vol1)+10); %Volume area above 10 m
m1_p(pair_nr)=nanmean(sedero(ind_vol1));
Vol2_p_10(pair_nr)=sum(altitude_i(ind_vol2)+10); %Volume area above 10 m
m2_p(pair_nr)=nanmean(sedero(ind_vol2));
Vol3_p_10(pair_nr)=sum(altitude_i(ind_vol3)+10); %Volume area above 10 m
m3_p(pair_nr)=nanmean(sedero(ind_vol3));
Vol4_p_10(pair_nr)=sum(altitude_i(ind_vol4)+10); %Volume area above 10 m
m4_p(pair_nr)=nanmean(sedero(ind_vol4));
Vol5_p_10(pair_nr)=sum(altitude_i(ind_vol5)+10); %Volume area above 10 m
m5_p(pair_nr)=nanmean(sedero(ind_vol5));
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_east_all (B_north-B_south) (B_west_all-B_east_all)],'LineWidth',2,'LineStyle','--')
rectangle('Position', [B_south B_east_1+2 (B_north-B_south) (B_west_1-B_east_1)-2],'EdgeColor','y')
rectangle('Position', [B_south B_east_2 (B_north-B_south) (B_west_2-B_east_2)-5],'EdgeColor','r')
rectangle('Position', [B_south B_east_3 (B_north-B_south) (B_west_3-B_east_3)-5],'EdgeColor','m')
rectangle('Position', [B_south B_east_4 (B_north-B_south) (B_west_4-B_east_4)-5],'EdgeColor','b')
rectangle('Position', [B_south B_east_5 (B_north-B_south) (B_west_5-B_east_5)-5],'EdgeColor','g')
text(B_south+10, B_west_1-30,'Area I','FontAngle','italic','color','w')
text(B_south+10, B_west_2-30,'Area II','FontAngle','italic','color','w')
text(B_south+10, B_west_3-30,'Area III','FontAngle','italic','color','k')
text(B_south+10, B_west_4-30,'Area IV','FontAngle','italic','color','k')
text(B_south+10, B_west_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_east_all (B_north-B_south) (B_west_all-B_east_all)],'LineWidth',2,'LineStyle','--')
rectangle('Position', [B_south B_east_1+2 (B_north-B_south) (B_west_1-B_east_1)-2],'EdgeColor','y')
rectangle('Position', [B_south B_east_2 (B_north-B_south) (B_west_2-B_east_2)-5],'EdgeColor','r')
rectangle('Position', [B_south B_east_3 (B_north-B_south) (B_west_3-B_east_3)-5],'EdgeColor','m')
rectangle('Position', [B_south B_east_4 (B_north-B_south) (B_west_4-B_east_4)-5],'EdgeColor','b')
rectangle('Position', [B_south B_east_5 (B_north-B_south) (B_west_5-B_east_5)-5],'EdgeColor','g')
text(B_south+10, B_west_1-30,'Area I','FontAngle','italic','color','w')
text(B_south+10, B_west_2-30,'Area II','FontAngle','italic','color','w')
text(B_south+10, B_west_3-30,'Area III','FontAngle','italic','color','k')
text(B_south+10, B_west_4-30,'Area IV','FontAngle','italic','color','k')
text(B_south+10, B_west_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

%% Figure Volume changes wrt to first survey in boxes
 figure

 bar(datenum(dates(pairs(:,1),:)),(m_p-m_p(1))*(B_west_all-B_east_all) )
hold on;
plot(datenum(dates(pairs(:,1),:)),(m1_p-m1_p(1))*(B_west_1-B_east_1),'*--y')
 plot(datenum(dates(pairs(:,1),:)),(m2_p-m2_p(1))*(B_west_2-B_east_2),'m')
 plot(datenum(dates(pairs(:,1),:)),(m3_p-m3_p(1))*(B_west_3-B_east_3),'r')
 plot(datenum(dates(pairs(:,1),:)),(m4_p-m4_p(1))*(B_west_4-B_east_4),'b')
 plot(datenum(dates(pairs(:,1),:)),(m5_p-m5_p(1))*(B_west_5-B_east_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 June 2009')
 legend('Total area','Area I','Area II','Area III','Area IV','Area V','location','northwest')

 set(gcf,'PaperPosition',[0  0  25 20])
print Vol_Areas.png -dpng

%% Figure Volume changes in boxes
 figure
 
for pair_nr=2:length(pairs);
 dV_all(pair_nr)=(m_p(pair_nr)-m_p(pair_nr-1))*(B_west_all-B_east_all);
 dV_1(pair_nr)=(m1_p(pair_nr)-m1_p(pair_nr-1))*(B_west_1-B_east_1);
 dV_2(pair_nr)=(m2_p(pair_nr)-m2_p(pair_nr-1))*(B_west_2-B_east_2);
 dV_3(pair_nr)=(m3_p(pair_nr)-m3_p(pair_nr-1))*(B_west_3-B_east_3);
 dV_4(pair_nr)=(m4_p(pair_nr)-m4_p(pair_nr-1))*(B_west_4-B_east_4);
 dV_5(pair_nr)=(m5_p(pair_nr)-m5_p(pair_nr-1))*(B_west_5-B_east_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 1
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_east_all (B_north-B_south) (B_west_all-B_east_all)],'LineWidth',2,'LineStyle','--')
rectangle('Position', [B_south B_east_1+2 (B_north-B_south) (B_west_1-B_east_1)-2],'EdgeColor','y')
rectangle('Position', [B_south B_east_2 (B_north-B_south) (B_west_2-B_east_2)-5],'EdgeColor','r')
rectangle('Position', [B_south B_east_3 (B_north-B_south) (B_west_3-B_east_3)-5],'EdgeColor','m')
rectangle('Position', [B_south B_east_4 (B_north-B_south) (B_west_4-B_east_4)-5],'EdgeColor','b')
rectangle('Position', [B_south B_east_5 (B_north-B_south) (B_west_5-B_east_5)-5],'EdgeColor','g')
text(B_south+10, B_west_1-30,'Area I','FontAngle','italic','color','w')
text(B_south+10, B_west_2-30,'Area II','FontAngle','italic','color','w')
text(B_south+10, B_west_3-30,'Area III','FontAngle','italic','color','k')
text(B_south+10, B_west_4-30,'Area IV','FontAngle','italic','color','k')
text(B_south+10, B_west_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

%% Figure Sed erosion 
if 1
figure(39)
for ipair=2:length(pairs)
pcolor(xi,yi,altitude_i(:,:,ipair)-altitude_i(:,:,1));  % 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_east_all (B_north-B_south) (B_west_all-B_east_all)],'LineWidth',2,'LineStyle','--')
rectangle('Position', [B_south B_east_1+2 (B_north-B_south) (B_west_1-B_east_1)-2],'EdgeColor','y')
rectangle('Position', [B_south B_east_2 (B_north-B_south) (B_west_2-B_east_2)-5],'EdgeColor','r')
rectangle('Position', [B_south B_east_3 (B_north-B_south) (B_west_3-B_east_3)-5],'EdgeColor','m')
rectangle('Position', [B_south B_east_4 (B_north-B_south) (B_west_4-B_east_4)-5],'EdgeColor','b')
rectangle('Position', [B_south B_east_5 (B_north-B_south) (B_west_5-B_east_5)-5],'EdgeColor','g')
text(B_south+10, B_west_1-30,'Area I','FontAngle','italic','color','w')
text(B_south+10, B_west_2-30,'Area II','FontAngle','italic','color','w')
text(B_south+10, B_west_3-30,'Area III','FontAngle','italic','color','k')
text(B_south+10, B_west_4-30,'Area IV','FontAngle','italic','color','k')
text(B_south+10, B_west_5-30,'Area V','FontAngle','italic','color','w')
axis equal
ylim([-50 1300])
caxis([-2 2])
xlabel ('along shore [m]')
ylabel ('cross shore [m]')
title ([dates(pairs(ipair,1),[4:6 1:3 7:end]) ' Survey' ])
colorbar
cmap=colormap(jet)
cmap(30:34,:)=1
colormap(cmap)
set(gcf,'Renderer','zbuffer')
set(gcf,'PaperPosition',[0  0  25 20])
filenameComb=['CombinedSedEro' dates(pairs(ipair,1),[7:end 1:2 4:5 ]) '.png']
print (filenameComb, '-dpng')
hold off
end
end