NCfile='http://dtvirt5.deltares.nl:8080/thredds/dodsC/opendap/rijkswaterstaat/jarkus/profiles/transect.nc';

transectid=11560; % ten noorden van de opgang

d07 = jarkus_readTransectDataNetcdf (NCfile,'Delfland',transectid, 2007);
d08 = jarkus_readTransectDataNetcdf (NCfile,'Delfland',transectid, 2008);
d09 = jarkus_readTransectDataNetcdf (NCfile,'Delfland',transectid, 2009);
d10 = jarkus_readTransectDataNetcdf (NCfile,'Delfland',transectid, 2010);


%% Figure profiles 2007to2010
figure
hold all
set(gca,'XDir','reverse')
xlim([00 1400])
ylim([-10 8])
xlabel(' cross shore distance [m]' )
ylabel(' elevation [m NAP]' )
ind07=find(~isnan(d07.zi)==1);
ind08=find(~isnan(d08.zi)==1);
ind09=find(~isnan(d09.zi)==1);
ind10=find(~isnan(d10.zi)==1);
plot(d07.xi(ind07),d07.zi(ind07))
plot(d08.xi(ind08),d08.zi(ind08))
plot(d09.xi(ind09),d09.zi(ind09))
plot(d10.xi(ind10),d10.zi(ind10))
grid on;   
legend('2007','2008','2009','2010','location','NorthWest')
title(['Delfland transect 11560' ])
set(gcf,'paperPosition',[0 0 15 10])
print 'Transect11560_2007to2010.jpg' -djpeg
 
%% Figure Total Nourished volume
b_off=1200;
b_on=0;
b_b=-10;
b_up=10;

ind_vol07=find(d07.xi(ind07) >= b_on & d07.xi(ind07) < b_off);
ind_vol10=find(d10.xi(ind10) >= b_on & d10.xi(ind10) < b_off);

[Vol_07] =jarkus_getVolume(d07.xi, d07.zi, b_up, b_b,b_on,b_off)
[Vol_10] =jarkus_getVolume(d10.xi, d10.zi, b_up, b_b,b_on,b_off)

figure
hold all
set(gca,'XDir','reverse')
xlim([00 1400])
ylim([-10 8])
xlabel(' cross shore distance [m]' )
ylabel(' elevation [m NAP]' )
plot(d07.xi(ind07),d07.zi(ind07),'LineWidth',2)
plot(d10.xi(ind10),d10.zi(ind10),'LineWidth',2)

grid on;   
legend('2007','2010','location','NorthWest')

% arcering
% for indh=ind10(ind_vol10(1)) : 5: ind10(ind_vol10(end))
%     plot([d10.xi(indh) d10.xi(indh)],[d07.zi(indh) d10.zi(indh)],'k')
% end
a=flipud(d10.xi(ind10(ind_vol10)))
xp=[d10.xi(ind10(ind_vol10)) ; a([1:67 69:2:107 109:end ]) ];
b=flipud(d07.zi(ind10(ind_vol10)))
yp=[d10.zi(ind10(ind_vol10 )); b([1:67 69:2:107 109:end ])];
zdata=ones(size(xp))
patch(xp,yp,zdata,'facecolor',1/255*[255 222 173],'FaceAlpha',0.85)

plot(tud1.xi(indtu1),tud1.zi(indtu1),'linewidth',2)  

LLWS=-0.93; % getij.nl
HHWS=LLWS+2.2; 
plot([320 1400],[320 1400]*0+LLWS,'--b','LineWidth',1)
plot([0 1400],[0 1400]*0,'b','LineWidth',1)
plot([210 1400],[210 1400]*0+HHWS,'--b','LineWidth',1)
grid on;   

text(1390,LLWS-0.4,'LLWS')
text(1390,HHWS+0.4,'HHWS')

%plot([b_off b_off],[-10 8],'--')
%plot([b_on b_on],[-10 8],'--')
%plot([0 1400],[3 3],'--')
%plot([0 1400],[b_b b_b],'--')

title(['Delfland transect 11560' ])
set(gcf,'paperPosition',[0 0 15 10])

set(gca,'ytick',[-10:2:10])
print 'Transect11560_NourishedVolume.jpg' -djpeg
 
dV_total_nourishment=Vol_10-Vol_07



%% Beun in met TUD orrsprong
% Figure Total Nourished volume
shift=300
figure
hold all
set(gca,'XDir','reverse')
xlim([00 1400]-shift)
ylim([-10 8])
xlabel(' cross shore distance [m]' )
ylabel(' elevation [m NAP]' )
plot(d07.xi(ind07)-shift,d07.zi(ind07),'LineWidth',2)
plot(d10.xi(ind10)-shift,d10.zi(ind10),'LineWidth',2)

grid on;   
legend('2007','2010','location','NorthWest')

% arcering
% for indh=ind10(ind_vol10(1)) : 5: ind10(ind_vol10(end))
%     plot([d10.xi(indh) d10.xi(indh)],[d07.zi(indh) d10.zi(indh)],'k')
% end
a=flipud(d10.xi(ind10(ind_vol10)))
xp=[d10.xi(ind10(ind_vol10)) ; a([1:67 69:2:107 109:end ]) ];
b=flipud(d07.zi(ind10(ind_vol10)))
yp=[d10.zi(ind10(ind_vol10 )); b([1:67 69:2:107 109:end ])];
zdata=ones(size(xp))
patch(xp-shift,yp,zdata,'facecolor',1/255*[255 222 173],'FaceAlpha',0.65)

a=flipud(d10.xi(ind10(ind_vol10)))
xp=[ a([1:67 69:2:107 109:end ]) ;0 ;1400; 1400];
b=flipud(d07.zi(ind10(ind_vol10)))
yp=[ b([1:67 69:2:107 109:end ]) ;-10 ;-10 ;-8.64];
zdata=ones(size(xp))
patch(xp-shift,yp,zdata,'facecolor',1/255*[139	69	19],'FaceAlpha',0.65)



%plot(tud1.xi(indtu1)-shift,tud1.zi(indtu1),'linewidth',2)  

LLWS=-0.93; % getij.nl
HHWS=LLWS+2.2; 
plot([320 1400]-shift,[320 1400]*0+LLWS,'--b','LineWidth',1)
plot([0 1400]-shift,[0 1400]*0,'b','LineWidth',1)
plot([210 1400]-shift,[210 1400]*0+HHWS,'--b','LineWidth',1)
grid on;   

text(1390-shift,LLWS-0.4,'LLWS')
text(1390-shift,HHWS+0.4,'HHWS')

%plot([b_off b_off],[-10 8],'--')
%plot([b_on b_on],[-10 8],'--')
%plot([0 1400],[3 3],'--')
%plot([0 1400],[b_b b_b],'--')


set(gcf,'paperPosition',[0 0 15 10])

set(gca,'ytick',[-10:2:10])








%% Total nourished volume in the beach and shoreface between -8 and 3m
b_off=1170;
b_on=110;
b_b=-8;
b_up=3;

ind_vol07=find(d07.xi(ind07) >= b_on & d07.xi(ind07) < b_off);
ind_vol10=find(d10.xi(ind10) >= b_on & d10.xi(ind10) < b_off);

[Vol_07] =jarkus_getVolume(d07.xi, d07.zi, b_up, b_b,b_on,b_off)
[Vol_10] =jarkus_getVolume(d10.xi, d10.zi, b_up, b_b,b_on,b_off)

figure
hold all
set(gca,'XDir','reverse')
xlim([00 1400])
ylim([-10 8])
xlabel(' cross shore distance [m]' )
ylabel(' elevation [m NAP]' )
plot(d07.xi(ind07),d07.zi(ind07))
% plot(d08.xi(ind08),d08.zi(ind08))
% plot(d09.xi(ind09),d09.zi(ind09))
plot(d10.xi(ind10),d10.zi(ind10))


grid on;   
legend('2007','2010','location','NorthWest')

for indh=ind10(ind_vol10(1)) : 5: ind10(ind_vol10(end))
    plot([d10.xi(indh) d10.xi(indh)],[d07.zi(indh) min(d10.zi(indh),3)],'k')
end

plot([b_off b_off],[-10 8],'--')
plot([b_on b_on],[-10 8],'--')
plot([0 1400],[3 3],'--')
plot([0 1400],[-8 -8],'--')

title(['Delfland transect 11560' ])
set(gcf,'paperPosition',[0 0 15 10])
print 'Transect11560_NourishedVolumeBeach.jpg' -djpeg
 
dV_min8tot3=Vol_10-Vol_07

%% Total nourished volume in the beach and shoreface between -8 and 3m
%% zeewaatrs van de nieuwe +3
b_off=1170;
b_on=375;
b_b=-8;
b_up=3;

ind_vol07=find(d07.xi(ind07) >= b_on & d07.xi(ind07) < b_off);
ind_vol10=find(d10.xi(ind10) >= b_on & d10.xi(ind10) < b_off);

[Vol_07] =jarkus_getVolume(d07.xi, d07.zi, b_up, b_b,b_on,b_off)
[Vol_10] =jarkus_getVolume(d10.xi, d10.zi, b_up, b_b,b_on,b_off)

figure
hold all
set(gca,'XDir','reverse')
xlim([00 1400])
ylim([-10 8])
xlabel(' cross shore distance [m]' )
ylabel(' elevation [m NAP]' )
plot(d07.xi(ind07),d07.zi(ind07))
plot(d10.xi(ind10),d10.zi(ind10))



grid on;   
legend('2007','2010','location','NorthWest')

for indh=ind10(ind_vol10(1)) : 5: ind10(ind_vol10(end))
    plot([d10.xi(indh) d10.xi(indh)],[d07.zi(indh) min(d10.zi(indh),3)],'k')
end

plot([b_off b_off],[-10 8],'--')
plot([b_on b_on],[-10 8],'--')
plot([0 1400],[3 3],'--')
plot([0 1400],[-8 -8],'--')

title(['Delfland transect 11560' ])
set(gcf,'paperPosition',[0 0 15 10])
print 'Transect11560_NourishedVolumeBeach2.jpg' -djpeg
 
dV_min8tot3zee=Vol_10-Vol_07


%% plot TUD data on Jarkus 11560 line
% load TUD data 
pathvlugtenburg='F:\check_outs\tudelft\vlugtenburg\'
NetCDFfilename=[pathvlugtenburg 'vlugtenburg.nc'];
%%
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');
X_RD= nc_varget(NetCDFfilename,'X_RD');
Y_RD= nc_varget(NetCDFfilename,'Y_RD');
%%
pair_nr=1;
[cross_shore_i,altitude_i] = combine_surveys_nr(survey_ind,pair_nr);
tud1.zi= griddata(X_RD,Y_RD,altitude_i',d10.xRD,d10.yRD);
tud1.xi=d10.xi;


pair_nr=11;
[cross_shore_i,altitude_i] = combine_surveys_nr(survey_ind,pair_nr);
tud10.zi= griddata(X_RD,Y_RD,altitude_i',d10.xRD,d10.yRD);
tud10.xi=d10.xi;

pair_nr=19;
[cross_shore_i,altitude_i] = combine_surveys_nr(survey_ind,pair_nr);
tud19.zi= griddata(X_RD,Y_RD,altitude_i',d10.xRD,d10.yRD);
tud19.xi=d10.xi;

pair_nr=23;
[cross_shore_i,altitude_i] = combine_surveys_nr(survey_ind,pair_nr);
tud23.zi= griddata(X_RD,Y_RD,altitude_i',d10.xRD,d10.yRD);
tud23.xi=d10.xi;

pair_nr=length(pairs);
[cross_shore_i,altitude_i] = combine_surveys_nr(survey_ind,pair_nr);
tudlast.zi= griddata(X_RD,Y_RD,altitude_i',d10.xRD,d10.yRD);
tudlast.xi=d10.xi;

%% Figure Vergelijk Jarkus and Jetski profiles
figure
hold all
set(gca,'XDir','reverse')
xlim([00 1400])
ylim([-10 8])
xlabel(' cross shore distance [m]' )
ylabel(' elevation [m NAP]' )
plot(d09.xi(ind09),d09.zi(ind09))
plot(d10.xi(ind10),d10.zi(ind10))
grid on; 
indtu1=find(~isnan(tud1.zi)==1);
plot(tud1.xi(indtu1),tud1.zi(indtu1),'linewidth',2)  
indtu10=find(~isnan(tud10.zi)==1);
plot(tud10.xi(indtu10),tud10.zi(indtu10),'linewidth',2)  

indtudlast=find(~isnan(tudlast.zi)==1);
plot(tudlast.xi(indtudlast),tudlast.zi(indtudlast),'linewidth',2)  
%legend('2007','2008','2009','2010','TUD june 2009','TUD may 2010','TUD feb 2011','location','NorthWest')
legend('Jarkus 2009','Jarkus 2010','TUD june 2009','TUD may 2010','TUD last','location','NorthWest')

title(['Delfland transect 11560, comparison Jarkus and TUD surveys interpolated on Jarkus grid' ])
set(gcf,'paperPosition',[0 0 15 10])
print 'Transect11560_TUD_and_Jarkus.jpg' -djpeg

%% Figure afname volume ter hoogte van de Jarkus Raai.
b_off=1210;
b_on=shift;
b_b=-10;
b_up=3;

% ind_vol07=find(d07.xi(ind07) >= b_on & d07.xi(ind07) < b_off);
ind_vol01=find(tud1.xi >= b_on & tud1.xi < b_off);

[Vol_t01] =jarkus_getVolume(tud1.xi,tud1.zi, b_up, b_b,b_on,b_off)
[Vol_t19] =jarkus_getVolume(tud19.xi,tud19.zi, b_up, b_b,b_on,b_off)

figure
hold all
set(gca,'XDir','reverse')
xlim([00 1400])
ylim([-10 8])
xlabel(' cross shore distance [m]' )
ylabel(' elevation [m NAP]' )
plot(tud1.xi,tud1.zi)
plot(tud19.xi,tud19.zi)

grid on;   
legend('June 2009','Feb 2011','location','NorthWest')

for indh=ind_vol01(1): 5: ind_vol01(end)
    if tud1.zi(indh)> tud19.zi(indh)
    plot([tud1.xi(indh) tud1.xi(indh)],[tud1.zi(indh) tud19.zi(indh)],'r')
    else
       plot([tud1.xi(indh) tud1.xi(indh)],[tud1.zi(indh) tud19.zi(indh)],'g') 
    end
end

plot([b_off b_off],[-10 8],'--')
plot([b_on b_on],[-10 8],'--')
plot([0 1400],[b_up b_up],'--')
plot([0 1400],[b_b b_b],'--')

title(['Delfland transect 11560, Meting TUD' ])
set(gcf,'paperPosition',[0 0 15 10])
print 'Transect11560_TUDJune09-Feb11.jpg' -djpeg
 
dV_Tud=Vol_t19-Vol_t01

%% Figure afname volume ter hoogte met raai 7
b_off=1210;
b_on=-150;
b_b=-10;
b_up=5;

transect_nr=7;

pair_nr=2;
[cross_shore_i,altitude_first] = combine_surveys_nr(survey_ind,pair_nr);
pair_nr=length(pairs)-1;
[cross_shore_i,altitude_last] = combine_surveys_nr(survey_ind,pair_nr);





% ind_vol07=find(d07.xi(ind07) >= b_on & d07.xi(ind07) < b_off);
ind_vol01=find(cross_shore_i >= b_on & cross_shore_i< b_off);

[Vol_t01] =jarkus_getVolume(cross_shore_i,altitude_first(:,transect_nr), b_up, b_b,b_on,b_off)
[Vol_t19] =jarkus_getVolume(cross_shore_i,altitude_last(:,transect_nr), b_up, b_b,b_on,b_off)

figure
hold all
set(gca,'XDir','reverse')
xlim([-200 1260])
ylim([-10 8])
xlabel(' cross shore distance [m]' )
ylabel(' elevation [m NAP]' )
plot(cross_shore_i,altitude_first(:,transect_nr),'-m','LineWidth',3)
plot(cross_shore_i,altitude_last(:,transect_nr),'-k','LineWidth',3)
%ind=find(altitude_first(:,transect_nr) <= 0);
%plot(cross_shore_i(ind),cross_shore_i(ind)*0,'-b','LineWidth',2)
LLWS=-0.93; % getij.nl
HHWS=LLWS+2.2; 
ind=find(altitude_first(:,transect_nr) <= LLWS);
plot(cross_shore_i(ind),cross_shore_i(ind)*0+LLWS,'--b','LineWidth',2)
ind=find(altitude_first(:,transect_nr) <= HHWS);
plot(cross_shore_i(ind),cross_shore_i(ind)*0+HHWS,'--b','LineWidth',2)
grid on;   

text(1190,LLWS-0.4,'Low Water line')
text(1190,HHWS-0.4,'High Water line')
legend('June 2009','March 2011','location','NorthWest')

for indh=ind_vol01(1): 2: ind_vol01(end)
    if altitude_first(indh,transect_nr) >= altitude_last(indh,transect_nr)
    plot([cross_shore_i(indh) cross_shore_i(indh)],[altitude_first(indh,transect_nr) altitude_last(indh,transect_nr)],'r')
    else
       plot([cross_shore_i(indh) cross_shore_i(indh)],[altitude_first(indh,transect_nr) altitude_last(indh,transect_nr)],'g') 
    end
end

 %plot([b_off b_off],[-10 8],'--')
 %plot([b_on b_on],[-10 8],'--')
 %plot([0 1400],[b_up b_up],'--')
 %plot([0 1400],[b_b b_b],'--')

title([' TUD transect nr. 7' ])
set(gcf,'paperPosition',[0 0 15 10])
print 'Transectnr7_TUDJune09-Feb11.jpg' -djpeg
 
dV_RaaiTUD=Vol_t19-Vol_t01




%% Afname volume met de TUD raaien
b_off=1210;
b_on=000;
b_b=-10;
b_up=3;

pair_nr=1;
[cross_shore_i,altitude_1] = combine_surveys_nr(survey_ind,pair_nr);

pair_nr=length(pairs);
[cross_shore_i,altitude_last] = combine_surveys_nr(survey_ind,pair_nr);

for raai_nr=1:22;
[Vol_t01] =jarkus_getVolume(cross_shore_i,altitude_1(:,raai_nr), b_up, b_b,b_on,b_off);
[Vol_tlast] =jarkus_getVolume(cross_shore_i,altitude_last(:,raai_nr), b_up, b_b,b_on,b_off);

dV_Tud_last(raai_nr)=Vol_tlast-Vol_t01;

end

pair_nr=11;
[cross_shore_i,altitude_firstyr] = combine_surveys_nr(survey_ind,pair_nr);

for raai_nr=1:22;
[Vol_t01] =jarkus_getVolume(cross_shore_i,altitude_1(:,raai_nr), b_up, b_b,b_on,b_off);
[Vol_firstyr] =jarkus_getVolume(cross_shore_i,altitude_firstyr(:,raai_nr), b_up, b_b,b_on,b_off);

dV_Tud_firstyr(raai_nr)=Vol_firstyr-Vol_t01;

end


%% Figure
figure
plot(dV_Tud_last(1,1:22));figure(gcf)
hold all;
plot(dV_Tud_firstyr,'DisplayName','dV_Tud_11','YDataSource','dV_Tud_11');figure(gcf)
xlabel('transect nummer (1 is Noord; 22 Zuid)')
ylabel('Berekend Volume verlies tussen 3 en -8 m NAP')
ylim([-400 100])
xlim([1 22])
grid on;
legend( 'Present-June 2009','May 2010-June 2009')

%% 
if 0
figure
ind=find(cross_shore_i >= 0 & cross_shore_i <= 1200);
surf(dist_alongshore,cross_shore_i(ind),altitude_19(ind,:))
hold on
plot3(cross_shore_i(ind)*0+dist_alongshore(10),cross_shore_i(ind),altitude_19(ind,10),'k','LineWidth',3)
plot3(cross_shore_i(ind)*0+dist_alongshore(10),cross_shore_i(ind),altitude_19(ind,10)*0,'-b','LineWidth',2)

[C,h] = contour3(dist_alongshore,cross_shore_i(ind),altitude_19(ind,:),[-9:.5:0 1.5 3],'--k');
clabel(C,h,'LabelSpacing',72,'FontSize',10,'Color','k','Rotation',0) ;

axis on

view([-110 40])

zlim([-10 5])
clim([-10 5])
shading interp

end