% DW_Plots 
%
% Deep Water Indicators and nourishment plots along the Dutch coast
%
% NB. indicators by tommer.vermaas@wur.nl
% Reference file DWI
% 
% giorgio.santinelli@deltares.nl

clear all
close all

%% read .nc files
% DWpath         = 'http://opendap.deltares.nl/thredds/dodsC/opendap/rijkswaterstaat/DW/';
DWpath         = 'd:\santinel\Giorgio\Kustlijnzorg\Tommer\';
NOURpath       = 'http://opendap.deltares.nl/thredds/dodsC/opendap/rijkswaterstaat/suppleties/';

DWurl          = 'DWI.nc';
NOURurl        = 'suppleties.nc';

id       = nc_varget([DWpath DWurl], 'id');

% read area
areacode = nc_varget([DWpath DWurl], 'areacode');
areaname = nc_varget([DWpath DWurl], 'areaname');

% read boundaries
low    = nc_varget([DWpath DWurl], 'lower_boundary');
up     = nc_varget([DWpath DWurl], 'upper_boundary');

% read time
DWyear   = nc_cf_time([DWpath DWurl]); % year

% read volume for different boundaries TOBECHANGED as preference
% DWvv2    = nc_varget([DWpath DWurl], 'vertical_volume', [0 0 2], [-1 -1 1]);
% DWvv4    = nc_varget([DWpath DWurl], 'vertical_volume', [0 0 4], [-1 -1 1]);
% DWvv5    = nc_varget([DWpath DWurl], 'vertical_volume', [0 0 5], [-1 -1 1]);
DWvv0    = nc_varget([DWpath DWurl], 'vertical_volume', [0 0 0], [-1 -1 1]);
DWvv1    = nc_varget([DWpath DWurl], 'vertical_volume', [0 0 1], [-1 -1 1]);
DWvv3    = nc_varget([DWpath DWurl], 'vertical_volume', [0 0 3], [-1 -1 1]);

% time re-definition to yyyy.mm...
DWtime  = str2num(datestr(DWyear,'yyyy'));

%% read .nc file of nourishments (m3/m)
NOUR_ntime = nc_varget([NOURpath NOURurl], 'n_time') + datenum('1970','yyyy'); % year of beg_date

NOUR_beach = nc_varget([NOURpath NOURurl], 'vol_beach');
NOUR_shore = nc_varget([NOURpath NOURurl], 'vol_shoreface');
NOUR_dune  = nc_varget([NOURpath NOURurl], 'vol_dune');
NOUR_other = nc_varget([NOURpath NOURurl], 'vol_other');

TIME_start_beach_n     = nc_varget([NOURpath NOURurl], 'time_start_beach') + datenum('1970','yyyy'); % actual date;
TIME_end_beach_n       = nc_varget([NOURpath NOURurl], 'time_end_beach') + datenum('1970','yyyy');
TIME_start_shore_n     = nc_varget([NOURpath NOURurl], 'time_start_shoreface') + datenum('1970','yyyy');
TIME_end_shore_n       = nc_varget([NOURpath NOURurl], 'time_end_shoreface') + datenum('1970','yyyy');
TIME_start_dune_n      = nc_varget([NOURpath NOURurl], 'time_start_dune') + datenum('1970','yyyy');
TIME_end_dune_n        = nc_varget([NOURpath NOURurl], 'time_end_dune') + datenum('1970','yyyy');
TIME_start_other_n     = nc_varget([NOURpath NOURurl], 'time_start_other') + datenum('1970','yyyy');
TIME_end_other_n       = nc_varget([NOURpath NOURurl], 'time_end_other') + datenum('1970','yyyy');

% definition in yyyy to avoid issues for datetick command
dateini = 1960;       
datefin = 2015;

TIME_start_beach = NaN(size(TIME_start_beach_n,1),size(TIME_start_beach_n,2)); 
TIME_end_beach   = NaN(size(TIME_end_beach_n,1),size(TIME_end_beach_n,2));
TIME_start_shore = NaN(size(TIME_start_shore_n,1),size(TIME_start_shore_n,2)); 
TIME_end_shore   = NaN(size(TIME_end_shore_n,1),size(TIME_end_shore_n,2)); 
TIME_start_dune  = NaN(size(TIME_start_dune_n,1),size(TIME_start_dune_n,2));   
TIME_end_dune    = NaN(size(TIME_end_dune_n,1),size(TIME_end_dune_n,2)); 
TIME_start_other = NaN(size(TIME_start_other_n,1),size(TIME_start_other_n,2)); 
TIME_end_other   = NaN(size(TIME_end_other_n,1),size(TIME_end_other_n,2)); 

% Indicators to inspect TOBECHANGED as preference
% DWI      = ['DWvv2';'DWvv4';'DWvv5'];
% boundLim = [num2str(up(3)) ' to ' num2str(low(3)); num2str(up(5)) ' to ' num2str(low(5)); num2str(up(6)) ' to ' num2str(low(6))];
DWI      = ['DWvv0';'DWvv1';'DWvv3'];
boundLim = [num2str(up(1)) ' to ' num2str(low(1)); num2str(up(2)) ' to ' num2str(low(2)); num2str(up(4)) ' to ' num2str(low(4))];



for ii = 1:size(DWI,1)
    DWvv = eval(DWI(ii,:));
    
    % define coefficient of linear regression + error for 0.95 conf. in DW
    p1      = NaN(length(id),1); p2      = NaN(length(id),1); p3      = NaN(length(id),1);
    uncert1 = NaN(length(id),1); uncert2 = NaN(length(id),1); uncert3 = NaN(length(id),1);
    
    DW_mean1 = NaN(length(id),1); DW_mean2 = NaN(length(id),1); DW_mean3 = NaN(length(id),1);
    
    %% Figures...
    DWfolder = [DWI(ii,:) '_Plots']; mkdir(DWfolder);

for i = 1:length(id)
    i
    if id(i)<8e6 && id(i)>7e6
    figure('Visible','Off')

    %% Nourishments + add axis

    hvol_st = []; hvol_on = []; hvol_dn = []; hvol_ot = [];
    
    for j = 1:length(NOUR_ntime)
        if ~isnan(NOUR_beach(i,j)) % beach nourishment
            if ~isnan(TIME_start_beach_n(i,j))
                TIME_start_beach(i,j) = str2double(datestr(TIME_start_beach_n(i,j),'yyyy')) + (str2double(datestr(TIME_start_beach_n(i,j),'mm'))-1)/12;
                TIME_end_beach(i,j) = str2double(datestr(TIME_end_beach_n(i,j),'yyyy')) + (str2double(datestr(TIME_end_beach_n(i,j),'mm'))-1)/12;
            end
            hvol_st = fill([TIME_start_beach(i,j), TIME_start_beach(i,j),...
                TIME_end_beach(i,j), TIME_end_beach(i,j)],...
                [0, NOUR_beach(i,j), NOUR_beach(i,j), 0],[1 0.69 0.39]); % fill the plot
            set(hvol_st,'EdgeColor',[1 0.69 0.39],'EdgeAlpha',0.2);
            set(gca,'XLimMode','manual','Xlim',[dateini datefin]);
            set(gca,'XTick', [],'XTickLabel', []);   
            hold on        
        end
        if ~isnan(NOUR_shore(i,j)) % shoreface nourishment
            if ~isnan(TIME_start_shore_n(i,j))
                TIME_start_shore(i,j) = str2double(datestr(TIME_start_shore_n(i,j),'yyyy')) + (str2double(datestr(TIME_start_shore_n(i,j),'mm'))-1)/12;
                TIME_end_shore(i,j) = str2double(datestr(TIME_end_shore_n(i,j),'yyyy')) + (str2double(datestr(TIME_end_shore_n(i,j),'mm'))-1)/12;
            end
            hvol_on = fill([TIME_start_shore(i,j), TIME_start_shore(i,j),...
                TIME_end_shore(i,j), TIME_end_shore(i,j)],...
                [0, NOUR_shore(i,j), NOUR_shore(i,j), 0],[0.04 0.52 0.78]);
            set(hvol_on,'EdgeColor',[0.04 0.52 0.78],'EdgeAlpha',0.2);
            set(gca,'XLimMode','manual','Xlim',[dateini datefin]);
            set(gca,'XTick', [],'XTickLabel', []);        
            hold on          
        end
        if ~isnan(NOUR_dune(i,j))  % dune nourishment
            if ~isnan(TIME_start_dune_n(i,j))
                TIME_start_dune(i,j) = str2double(datestr(TIME_start_dune_n(i,j),'yyyy')) + (str2double(datestr(TIME_start_dune_n(i,j),'mm'))-1)/12;
                TIME_end_dune(i,j) = str2double(datestr(TIME_end_dune_n(i,j),'yyyy')) + (str2double(datestr(TIME_end_dune_n(i,j),'mm'))-1)/12;
            end            
            hvol_dn = fill([TIME_start_dune(i,j), TIME_start_dune(i,j),...
                TIME_end_dune(i,j), TIME_end_dune(i,j)],...
                [0, NOUR_dune(i,j), NOUR_dune(i,j), 0],[0.5 0.5 0]);
            set(hvol_dn,'EdgeColor',[0.5 0.5 0],'EdgeAlpha',0.2);
            set(gca,'XLimMode','manual','Xlim',[dateini datefin]);
            set(gca,'XTick', [],'XTickLabel', []);        
            hold on           
        end
        if ~isnan(NOUR_other(i,j)) % other type of nourishment
            if ~isnan(TIME_start_other_n(i,j))
                TIME_start_other(i,j) = str2double(datestr(TIME_start_other_n(i,j),'yyyy')) + (str2double(datestr(TIME_start_other_n(i,j),'mm'))-1)/12;
                TIME_start_other(i,j) = str2double(datestr(TIME_start_other_n(i,j),'yyyy')) + (str2double(datestr(TIME_start_other_n(i,j),'mm'))-1)/12;
            end            
            hvol_ot = fill([TIME_start_other(i,j), TIME_start_other(i,j),...
                TIME_end_other(i,j), TIME_end_other(i,j)],...
                [0, NOUR_other(i,j), NOUR_other(i,j), 0],[1 0 0]);
            set(hvol_ot,'EdgeColor',[1 0 0],'EdgeAlpha',0.2);
            set(gca,'XLimMode','manual','Xlim',[dateini datefin]);
            set(gca,'XTick', [],'XTickLabel', []);
            hold on            
        end    
    end
    ylabel('Nourishment volume [m  ^3/m]','FontSize',12)
    set(gca,'XTick', [],'XTickLabel', []);    

    % routine to display legend [dumb...]
    Mhvol     = [' Beach    ';' Shoreface';' Dune     ';' Other    '];
    if isempty(hvol_st); Mhvol(find(ismember(Mhvol,' Beach    ', 'rows')),:) = []; end
    if isempty(hvol_on); Mhvol(find(ismember(Mhvol,' Shoreface', 'rows')),:) = []; end
    if isempty(hvol_dn); Mhvol(find(ismember(Mhvol,' Dune     ', 'rows')),:) = []; end
    if isempty(hvol_ot); Mhvol(find(ismember(Mhvol,' Other    ', 'rows')),:) = []; end

    % do not display ticks and labels if no nourishment occurred
    if isempty(hvol_st) && isempty(hvol_on) && isempty(hvol_dn) && isempty(hvol_ot)
        set(gca,'YTick', [],'YTickLabel', []);
    end

    hvol      = ([hvol_st; hvol_on; hvol_dn; hvol_ot]); 
    legend(hvol, char(Mhvol), 'Location', 'SouthWest');
    HandleNou = gca;

    % set # of ticks to 10
    ylimits = get(gca,'YLim');
    yinc    = (ylimits(2)-ylimits(1))/10;
    set(gca,'YTick', ylimits(1):yinc:ylimits(2));
    
    %% DW

    HandleDW = axes('Position',get(HandleNou,'Position'));
    
    D = plot(DWtime, DWvv(:,i),'bs','LineWidth',2,'MarkerSize',7);
    hold on
    
    set(HandleDW,'YAxisLocation','right','Color','none');
    set(HandleDW,'XLimMode','manual','Xlim',[dateini datefin]);

    % set # of ticks to 10
    ylimits = get(gca,'YLim');
    yinc    = (ylimits(2)-ylimits(1))/10;
    set(gca,'YTick', ylimits(1):yinc:ylimits(2));

    % figure opt
    set(gcf,'Position',[5 35 1040 760]); 
    set(gca,'FontSize',12); grid on
    xlabel('Year','FontSize',12); ylabel('Vertical volume [m  ^3/m]','FontSize',12)
    title(['Volume Indicator  ' boundLim(ii,:) ' m : ' char(cellstr(areaname(i,:))) ' (' num2str(areacode(i)) '), transect '...
        num2str(id(i) - areacode(i)*1000000)],'FontSize',12,'FontWeight','bold')
    legendDW = legend(D, {' Volume'}, 'Location', 'NorthWest');  
    
    print('-dpng','-r300',[DWfolder filesep char(cellstr(areaname(i,:))) '_' num2str(id(i))])

    close all
    clear('hvol_st','hvol_on','hvol_dn','hvol_ot','D');
    
    end

end

clear('DWvv');
    
end
