% FK_Plots
%
% Coast line indicators and nourishment plots along the Dutch coast:
% - probability of failure (faalkans) - last update 2017
%
%
%	---------------------------------------
%
%       Giorgio Santinelli
%
%       Giorgio.Santinelli@deltares.nl
%
%       Deltares
%       P.O. Box 177
%       2600 MH Delft
%       The Netherlands
%
%	---------------------------------------
%
%
% $Id: FK_Plots.m 9089 2018-01-08 08:41:32Z santinel $
% $Date: 2018-01-08 09:41:32 +0100 (Mon, 08 Jan 2018) $
% $Author: santinel $
% $Revision: 9089 $
% $HeadURL: https://svn.oss.deltares.nl/repos/openearthrawdata/trunk/rijkswaterstaat/coastline_graphs_beheerregister/scripts/FK_Plots.m $
% $Keywords$
%

clear all
close all

%% read .nc files
FKpath         = 'http://opendap.deltares.nl/thredds/dodsC/opendap/rijkswaterstaat/faalkans/';
% FKpath         = 'n:\Projects\1204000\1204421\B. Measurements and calculations\Toestand van de kust\Analysis\probability_of_failure\Netcdf bestanden\Netcdf file\';

NOURpath       = 'http://opendap.deltares.nl/thredds/dodsC/opendap/rijkswaterstaat/suppleties/';

FKurl          = 'faalkans.nc';
NOURurl        = 'nourishments.nc';

id       = nc_varget([FKpath FKurl], 'id');

% read area
areacode = nc_varget([FKpath FKurl], 'areacode');
areaname = nc_varget([FKpath FKurl], 'areaname');

% read time
FKtime   = nc_cf_time([FKpath FKurl]); % year

FK       = nc_varget([FKpath FKurl], 'probability_failure');

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

%% read .nc file of nourishments (m3/m)
NOUR_ntime = str2num(datestr(nc_cf_time([NOURpath NOURurl], 'time'),'yyyy')); % year of beg_date

NOUR_beach = nc_varget([NOURpath NOURurl], 'volume',[0 0 0],[-1 -1 1]); % beach
NOUR_shore = nc_varget([NOURpath NOURurl], 'volume',[0 0 1],[-1 -1 1]); % shoreface
NOUR_dune  = nc_varget([NOURpath NOURurl], 'volume',[0 0 2],[-1 -1 1]); % dune
NOUR_other = nc_varget([NOURpath NOURurl], 'volume',[0 0 3],[-1 -1 1]); % other

TIME_start_beach_n     = nc_varget([NOURpath NOURurl], 'time_start', [0 0 0],[-1 -1 1]) + datenum('1970','yyyy'); % actual date;
TIME_end_beach_n       = nc_varget([NOURpath NOURurl], 'time_end', [0 0 0],[-1 -1 1]) + datenum('1970','yyyy');
TIME_start_shore_n     = nc_varget([NOURpath NOURurl], 'time_start', [0 0 1],[-1 -1 1]) + datenum('1970','yyyy');
TIME_end_shore_n       = nc_varget([NOURpath NOURurl], 'time_end', [0 0 1],[-1 -1 1]) + datenum('1970','yyyy');
TIME_start_dune_n      = nc_varget([NOURpath NOURurl], 'time_start', [0 0 2],[-1 -1 1]) + datenum('1970','yyyy');
TIME_end_dune_n        = nc_varget([NOURpath NOURurl], 'time_end', [0 0 2],[-1 -1 1]) + datenum('1970','yyyy');
TIME_start_other_n     = nc_varget([NOURpath NOURurl], 'time_start', [0 0 3],[-1 -1 1]) + datenum('1970','yyyy');
TIME_end_other_n       = nc_varget([NOURpath NOURurl], 'time_end', [0 0 3],[-1 -1 1]) + datenum('1970','yyyy');

%% Figures...
FKfolder = 'FK_Plots_withtrends';
mkdir(FKfolder);

% nourishment 'fixed' duration set to
nourduration = 0.5; % halfayear

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

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));

for i = 1:length(id)
    % if ismember(i,find(id>=7e6 & id<8e6))
    disp([num2str(i) '/' num2str(length(id))])
    figure('Visible','Off'); hold on;
    
    %% 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_start_beach(i,j) + nourduration, TIME_start_beach(i,j) + nourduration],...
                [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_start_shore(i,j) + nourduration, TIME_start_shore(i,j) + nourduration],...
                [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_start_dune(i,j) + nourduration, TIME_start_dune(i,j) + nourduration],...
                [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_start_other(i,j) + nourduration, TIME_start_other(i,j) + nourduration],...
                [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 = [0 1200]; % 1200 is a nice lcm for visualization purposes
    
    %%MAKE IT LOG
    HandleFK = axes('Position',get(HandleNou,'Position'));
    
    % compute linear trend 1965 - 1990
    fk1 = log10(FK(FKtime<1990,i));
    fk1time = FKtime(FKtime<1990);
    if sum(~isnan(fk1))>6 % ...arbirtary number of points to allow for regression
        b1 = polyfit(fk1time(~isnan(fk1)),fk1(~isnan(fk1)),1); fk1_trend = polyval(b1, fk1time(~isnan(fk1)));
        TR1 = semilogy(fk1time(~isnan(fk1)),10.^(fk1_trend),'-k'); hold on; % params are b1
    end
    fk2 = log10(FK(FKtime>=1990,i));
    fk2time = FKtime(FKtime>=1990);
    if sum(~isnan(fk2))>6 % ...arbirtary number of points to allow for regression
        b2 = polyfit(fk2time(~isnan(fk2)),fk2(~isnan(fk2)),1); fk2_trend = polyval(b2, fk2time(~isnan(fk2)));
        TR2 = semilogy(fk2time(~isnan(fk2)),10.^(fk2_trend),'-k'); hold on; % params are b2
    end
        
    % plot values
    F = semilogy(FKtime, FK(:,i),'ks','LineWidth',2,'MarkerFaceColor','k','MarkerSize',5);
    set(gca,'YAxisLocation','right','Color','none');
    set(gca,'XLimMode','manual','Xlim',[dateini datefin]);
    
    yinc    = (ylimits(2)-ylimits(1))/(length(get(HandleFK,'YTick'))-1);
    set(HandleNou,'YLimMode','manual','YLim',ylimits,'YTick',ylimits(1):yinc:ylimits(2),...
        'YTickLabel', num2str((ylimits(1):yinc:ylimits(2))'));
    set(HandleFK,'YMinorTick','off')
    
    % figure opt
    set(gcf,'Position',[5 35 1040 760]);
    set(gca,'FontSize',12); grid on
    xlabel('Year','FontSize',12); ylabel('Probability','FontSize',12)
    title(['Probability of failure : ' char(cellstr(areaname(i,:))) ' (' num2str(areacode(i)) '), transect '...
        num2str(id(i) - areacode(i)*1000000)],'FontSize',12,'FontWeight','bold')
    legendKL = legend([F], {[' Faalkans']}, 'Location', 'NorthWest');
    if strcmp(char(cellstr(areaname(i,:))),'Maasvlakte / slufter'); areaname(i,1:20) = 'Maasvlakte          '; end;
    if strcmp(char(cellstr(areaname(i,:))),'Oosterschelde / Neeltje Jans'); areaname(i,1:28) = 'Oosterschelde               '; end;
    print('-dpng','-r300',[FKfolder filesep char(cellstr(areaname(i,:))) '_' num2str(id(i))])
    
    close all
    clear('hvol_st','hvol_on','hvol_dn','hvol_ot','F');
    
    %   end
end