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

%% control zones
% alongshore boundaries
B_north  = 1500;
B_south  = 1050;
% 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:19;
    [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 Volume changes 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([-1.2 1])
%  xlim([ datenum(dates(pairs(1,1),:)) datenum(dates(pairs(end,1),:)) ])
%   datetick('x')
%  grid on
%  ylabel('\Delta V per m alongshore')
%  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 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 [m]')
%  xlabel('Time')
%  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


%% Figure first and last topo and areas of interest
figure

subplot(1,2,[1])
pcolor(xi,yi,altitude_i(:,:,end));  % 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])
colorbar('location', 'eastOutside')
xlabel ('along shore [m]')
ylabel ('cross shore [m]')
title ([dates(pairs(end,1),[1:2 6:end]) ' survey' ])

subplot(1,2,[2 ])
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')
legend('Total area','Area I','Area II','Area III','Area IV','Area V','location','Best')

set(gcf,'position',[181         669        1049         253]) 
set(gcf,'PaperPosition',[0  0  25 6])
print Overview_Vol_Areas_North.png -dpng