% plot volumes & differences between years & correct for any suppletions,
% if necessary
clear, clc;

TargetIDs = [...
    5004898 ...
    6002011 ...
    7003975 ...
    8008950 ...
    9011450];

SuppletionURL = 'http://dtvirt5.deltares.nl:8080/thredds/dodsC/opendap/rijkswaterstaat/suppleties/suppleties.nc';

id          = nc_varget('http://dtvirt5.deltares.nl:8080/thredds/dodsC/opendap/rijkswaterstaat/suppleties/suppleties.nc','id');
stretch     = nc_varget('http://dtvirt5.deltares.nl:8080/thredds/dodsC/opendap/rijkswaterstaat/suppleties/suppleties.nc','stretch');
date        = nc_varget('http://dtvirt5.deltares.nl:8080/thredds/dodsC/opendap/rijkswaterstaat/suppleties/suppleties.nc','date')+ datenum('1970-01-01 00:00');
type        = nc_varget('http://dtvirt5.deltares.nl:8080/thredds/dodsC/opendap/rijkswaterstaat/suppleties/suppleties.nc','type');
typeflag    = nc_varget('http://dtvirt5.deltares.nl:8080/thredds/dodsC/opendap/rijkswaterstaat/suppleties/suppleties.nc','type_flag');

NrSuppletions       = size(stretch,1);
RelevantSuppletions = false(NrSuppletions,numel(TargetIDs));

% figure out if suppletions were done at the target locations
for iSuppletion = 1:NrSuppletions
    IDfilter    = id >= stretch(iSuppletion,1) & id <= stretch(iSuppletion,2);
    AffectedIDs = id(IDfilter);
    
    for iTarget = 1:numel(TargetIDs)
        TargetAffectedIdx   = AffectedIDs == TargetIDs(iTarget);
        if any(TargetAffectedIdx)
            RelevantSuppletions(iSuppletion,iTarget) = true;
        end
    end
end

NCfile = 'n:\Projects\1209000\1209436\F. Other information\Temporele Effecten\calculations\volumevar.nc';

id      = nc_varget(NCfile,'id');
Vols    = nc_varget(NCfile,'volumeC');
time    = nc_varget(NCfile,'time') + datenum('1970-01-01 00:00');
Years   = 1965:2013;
Years   = Years';

dVoud       = [ ...
    272 ...
    369 ...
    411 ...
    371 ...
    295];

TTrend = 2;
OTrend = 1;

for iLoc = 1:numel(TargetIDs)
    LocID = find(id == TargetIDs(iLoc));
    
    dVtrend = NaN(numel(time)-1,1);
    dVobs   = NaN(numel(time)-1,1);
    Vtrend  = NaN(numel(time),1);
    for iYear = 2:numel(time)
        fit     = [];
        filter  = [];
        filter  = time >= time(iYear)-TTrend*365 & time <= time(iYear)+TTrend*365; 
        fit     = polyfit(time(filter),Vols(filter,LocID),OTrend);
        dVtrend(iYear-1)    = fit(1)*365;
        Vtrend(iYear)       = polyval(fit,time(iYear));
        
        if iYear == 6
            subplot(1,2,1)
            VtrendL             = polyval(fit,time);
            plot(time,Vols(:,LocID),'.-','DisplayName','Jarkus volumes')
            hold on
            plot(time,VtrendL,'k--','DisplayName','Trend')
            scatter(time(iYear),Vtrend(iYear),'ro','DisplayName','Target year')
            hold off
            legend('show')
            datetick(gca)
            xlabel('Years')
            ylabel('Volume in active dune erosion zone [m^3/m]')
            YLims = get(gca,'YLim');
        end
    end
    
    dVobs   = diff(Vols(:,LocID));
    dVcorr  = dVobs-dVtrend;
    Vdetrend = detrend(Vols(:,LocID));
    
    % filter out the beach and dune suppletions
    SupTimes    = date(RelevantSuppletions(:,iLoc),:);
    SupTypes    = typeflag(RelevantSuppletions(:,iLoc),:);
    
    TypeFilter  = SupTypes == 1 | SupTypes == 3 | SupTypes == 2;
    
    SupTimeIdx  = false(size(time,1)-1,1);
    for iSup = 1:sum(TypeFilter)
        idxStart    = find(time > SupTimes(iSup,1)-14,1,'first')-1;
        idxStop     = find(time < SupTimes(iSup,2)+14,1,'last');
        SupTimeIdx(idxStart:idxStop,1)  = true;
    end
    
    % Disregard reinforcement for 9011450
    if TargetIDs(iLoc) == 9011450
        idxReinforcement    = find(time > datenum('2009-01-01 00:00'));
        SupTimeIdx(idxReinforcement,1)  = true;
    end
    
    [bandwidth,density,xmesh,cdf] = kde(dVcorr(~SupTimeIdx));
    
    NormProbDens = norm_pdf(xmesh,mean(dVcorr(~SupTimeIdx)),std(dVcorr(~SupTimeIdx)));
    
    fprintf('Transect %d mean %f std %f \n',TargetIDs(iLoc),mean(dVcorr(~SupTimeIdx)),std(dVcorr(~SupTimeIdx)));
    
%     subplot(2,2,1)
%     plot(time,Vols(:,LocID),'.-','DisplayName','Jarkus volumes')
%     hold on
%     plot(time,Vtrend,'k--','DisplayName','Trend')
%     SupDates = date(RelevantSuppletions(:,iLoc),:);
%     XLims = get(gca,'XLim');
%     YLims = get(gca,'YLim');
%     patch([SupDates'; flipud(SupDates')],repmat([YLims(1); YLims(1); YLims(2); YLims(2)],1,size(SupDates,1)),'r','DisplayName','Nourishments')
%     hold off
%     ch = get(gca,'children');
%     set(gca,'children', flipud(ch))
%     datetick(gca)
%     legend('show','location','NorthWest')
%     
%     subplot(2,2,2)
%     hist(dVcorr(~SupTimeIdx),-500:50:500)
%     set(gca,'XLim',[-500 500]);
%     
%     subplot(2,2,3)
%     plot(time(2:end)-0.5,dVobs,'b.-','DisplayName','dV_{obs}')
%     hold on
%     plot(time([false; ~SupTimeIdx])-0.5,dVcorr(~SupTimeIdx),'k.-','DisplayName','dV_{corr}')
%     YLims = get(gca,'YLim');
%     
%     patch([SupDates'; flipud(SupDates')],repmat([YLims(1); YLims(1); YLims(2); YLims(2)],1,size(SupDates,1)),'r','DisplayName','Nourishments')
% %     plot(Years,Vdetrend,'m.-','DisplayName','dV_{detrended}')
%     hold off
%     ch = get(gca,'children');
%     set(gca,'children', flipud(ch))
%     set(gca,'XLim',XLims);
%     datetick(gca)
%     grid on
%     legend('show','location','NorthWest')
%     
%     subplot(2,2,4)
%     plot(xmesh,density,'b','DisplayName','ProbDens_{kde}')
%     hold on
%     plot(xmesh,NormProbDens,'k','DisplayName','Normal dist.')
%     plot([-dVoud(iLoc) 0 dVoud(iLoc)],[0 1/dVoud(iLoc) 0],'m','DisplayName','Triangular dist.')
%     hold off
%     set(gca,'XLim',[-500 500]);
%     set(gca,'YLim',[0 8e-3]);
%     legend('show')
    
    subplot(1,2,2)
    plot(time,Vols(:,LocID),'.-','DisplayName','Jarkus volumes')
    hold on
    plot(time,Vtrend,'k--','DisplayName','Trend')
    SupDates = date(RelevantSuppletions(:,iLoc),:);
    XLims = get(gca,'XLim');
%     YLims = get(gca,'YLim');
    patch([SupDates'; flipud(SupDates')],repmat([YLims(1); YLims(1); YLims(2); YLims(2)],1,size(SupDates,1)),'r','DisplayName','Nourishments')
    hold off
    ch = get(gca,'children');
    set(gca,'children', flipud(ch))
    datetick(gca)
    legend('show','location','NorthWest')
    xlabel('Years')
    ylabel('Volume in active dune erosion zone [m^3/m]')
    print('-dpng','-r600','nourishments_volumes')
 end