% WL_Plots 
%
% Waterline and nourishment plots along the Dutch coast
%
% Two different primary sources are used:
%		Jarkus database (transects.nc)		1965 - today
%		[old] database  (strandlijnen.nc) 	    1843 - 1998
%
% MHW, MLW, DF have been created using the [old] database strandlijnen.nc for period 1843 - 1965
% Available in http://opendap.deltares.nl/thredds/catalog/opendap/rijkswaterstaat/
% 
% DF in the plots is defined as the intersection of the profile with +3 NAP
% 
%	---------------------------------------
% 
%       Giorgio Santinelli
%
%       Giorgio.Santinelli@deltares.nl
%
%       Deltares
%       P.O. Box 177
%       2600 MH Delft
%       The Netherlands
%
%	---------------------------------------
%
% $Id: WL_Plots.m 9013 2016-12-21 08:00:07Z santinel $
% $Date: 2016-12-21 09:00:07 +0100 (Wed, 21 Dec 2016) $
% $Author: santinel $
% $Revision: 9013 $
% $HeadURL: https://svn.oss.deltares.nl/repos/openearthrawdata/trunk/rijkswaterstaat/coastline_graphs_beheerregister/scripts/WL_Plots.m $

clear
close all

%% read .nc files
DFpath           = 'http://opendap.deltares.nl/thredds/dodsC/opendap/rijkswaterstaat/DuneFoot/';
% local
%DFpath           = 'z:\D\OpenEarthRawData\rijkswaterstaat\DuneFoot\processed\';
MHW_MLWpath      = 'http://opendap.deltares.nl/thredds/dodsC/opendap/rijkswaterstaat/MHW_MLW/';
% local
%MHW_MLWpath      = 'z:\d\OpenEarthRawData\rijkswaterstaat\MHW_MLW\processed\';
% online nc throw the following error: 'NC_UNLIMITED in the wrong index (NC_EUNLIMPOS)'
% nourishments.nc online dataset has to be fixed.
NOURpath       = 'http://opendap.deltares.nl/thredds/dodsC/opendap/rijkswaterstaat/suppleties/'; 

MHW_MLWurl        = 'MHW_MLW.nc';
DFurl            = 'DF.nc';
% NOURurl        = 'suppleties.nc';
NOURurl        = 'nourishments.nc';

id       = nc_varget([MHW_MLWpath MHW_MLWurl], 'id'); % using DF for this and further vars doesn't make difference

% read area
areacode = nc_varget([MHW_MLWpath MHW_MLWurl], 'areacode');
areaname = nc_varget([MHW_MLWpath MHW_MLWurl], 'areaname');

% read time (1 - Jan - XXXX represents year XXXX on the plots)
time_JARKUS = nc_cf_time([MHW_MLWpath MHW_MLWurl]); % year

% DF as profile intersection with +3NAP
MHW      = nc_varget([MHW_MLWpath MHW_MLWurl], 'mean_high_water_cross');
MLW      = nc_varget([MHW_MLWpath MHW_MLWurl], 'mean_low_water_cross');
DF       = nc_varget([DFpath DFurl], 'dune_foot_threeNAP_cross');

MHW_JARKUS  = MHW';
MLW_JARKUS  = MLW';
DF_JARKUS   = DF';

% time re-definition to yyyy.mm...
time_JARKUS = str2num(datestr(time_JARKUS,'yyyy')); % + (str2num(datestr(time_JARKUS,'mm'))-1)/12;

% distinguish between new and old database
commonID = false(size(MHW_JARKUS,1),1);
for ii = 1:size(MHW_JARKUS,1)
    if any(~isnan(MHW_JARKUS(ii,time_JARKUS<1965))) ||...
        any(~isnan(MLW_JARKUS(ii,time_JARKUS<1965))) ||...
        any(~isnan(DF_JARKUS(ii,time_JARKUS<1965)));
        commonID(ii) = true;
    end
end

%% 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...
WLfolder = 'WL_Plots_withtrends';
mkdir(WLfolder);

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

% definition in yyyy to avoid issues for datetick command
dateini     = 1960;       datefin = 2020;
dateini_old = 1840;   datefin_old = 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)

    disp([num2str(i) '/' num2str(length(id))])
    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_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);
            if commonID(i) % both old and new database
                set(gca,'XLimMode','manual','Xlim',[dateini_old datefin_old]);
            else
                set(gca,'XLimMode','manual','Xlim',[dateini datefin]);
            end
            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);
            if commonID(i)
                set(gca,'XLimMode','manual','Xlim',[dateini_old datefin_old]);
            else
                set(gca,'XLimMode','manual','Xlim',[dateini datefin]);
            end
            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);
            if commonID(i)
                set(gca,'XLimMode','manual','Xlim',[dateini_old datefin_old]);
            else
                set(gca,'XLimMode','manual','Xlim',[dateini datefin]);
            end
            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);
            if commonID(i)
                set(gca,'XLimMode','manual','Xlim',[dateini_old datefin_old]);
            else
                set(gca,'XLimMode','manual','Xlim',[dateini datefin]);
            end
            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));

    
    % % % Not sure if it makes sense to implement this, because of the timespan.     
    %% 	MHW, MLW, DF
    HandleWL = axes('Position',get(HandleNou,'Position')); % use part of the same handle

    % compute linear trend 1965 - 1990
    mlw1 = MLW_JARKUS(i,time_JARKUS<1990 & time_JARKUS>=1965); mlw1time = time_JARKUS(time_JARKUS<1990 & time_JARKUS>=1965)';
    mhw1 = MHW_JARKUS(i,time_JARKUS<1990 & time_JARKUS>=1965); mhw1time = time_JARKUS(time_JARKUS<1990 & time_JARKUS>=1965)';
    df1 = DF_JARKUS(i,time_JARKUS<1990 & time_JARKUS>=1965); df1time = time_JARKUS(time_JARKUS<1990 & time_JARKUS>=1965)';
    if sum(~isnan(mlw1))>6 % ...arbirtary number of points to allow for regression
        Lb1 = polyfit(mlw1time(~isnan(mlw1)),mlw1(~isnan(mlw1)),1); mlw1_trend = polyval(Lb1, mlw1time(~isnan(mlw1)));
        LTR1 = plot(mlw1time(~isnan(mlw1)),mlw1_trend,'-b'); hold on; % params are b1
    end
    if sum(~isnan(mhw1))>6 % ...arbirtary number of points to allow for regression
        Hb1 = polyfit(mhw1time(~isnan(mhw1)),mhw1(~isnan(mhw1)),1); mhw1_trend = polyval(Hb1, mhw1time(~isnan(mhw1)));
        HTR1 = plot(mhw1time(~isnan(mhw1)),mhw1_trend,'-r'); hold on; % params are b1
    end
    if sum(~isnan(df1))>6 % ...arbirtary number of points to allow for regression
        Db1 = polyfit(df1time(~isnan(df1)),df1(~isnan(df1)),1); df1_trend = polyval(Db1, df1time(~isnan(df1)));
        DTR1 = plot(df1time(~isnan(df1)),df1_trend,'-','Color',[0 0.5 0]); hold on; % params are b1
    end
    mlw2 = MLW_JARKUS(i,time_JARKUS>=1990); mlw2time = time_JARKUS(time_JARKUS>=1990)';
    mhw2 = MHW_JARKUS(i,time_JARKUS>=1990); mhw2time = time_JARKUS(time_JARKUS>=1990)';
    df2 = DF_JARKUS(i,time_JARKUS>=1990); df2time = time_JARKUS(time_JARKUS>=1990)';
    if sum(~isnan(mlw2))>6 % ...arbirtary number of points to allow for regression
        Lb2 = polyfit(mlw2time(~isnan(mlw2)),mlw2(~isnan(mlw2)),1); mlw2_trend = polyval(Lb2, mlw2time(~isnan(mlw2)));
        LTR2 = plot(mlw2time(~isnan(mlw2)),mlw2_trend,'-b'); hold on; % params are b2
    end
    if sum(~isnan(mhw2))>6 % ...arbirtary number of points to allow for regression
        Hb2 = polyfit(mhw2time(~isnan(mhw2)),mhw2(~isnan(mhw2)),1); mhw2_trend = polyval(Hb2, mhw2time(~isnan(mhw2)));
        HTR2 = plot(mhw2time(~isnan(mhw2)),mhw2_trend,'-r'); hold on; % params are b2
    end
    if sum(~isnan(df2))>6 % ...arbirtary number of points to allow for regression
        Db2 = polyfit(df2time(~isnan(df2)),df2(~isnan(df2)),1); df2_trend = polyval(Db2, df2time(~isnan(df2)));
        DTR2 = plot(df2time(~isnan(df2)),df2_trend,'-','Color',[0 0.5 0]); hold on; % params are b2
    end
    
    % plot parameters
    L = plot(time_JARKUS, MLW_JARKUS(i,:),'bx','LineWidth',2,'MarkerSize',6);    
    hold on
    H = plot(time_JARKUS, MHW_JARKUS(i,:),'rx','LineWidth',2,'MarkerSize',6);
    D = plot(time_JARKUS, DF_JARKUS(i,:),'x','LineWidth',2,'MarkerSize',6,'Color',[0 0.5 0]);
 
    set(HandleWL,'YAxisLocation','right','Color','none');
    if commonID(i)
        set(gca,'XLimMode','manual','Xlim',[dateini_old datefin_old]);
    else
        set(gca,'XLimMode','manual','Xlim',[dateini datefin]);
    end

    % 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('Distance wrt RSP [m]','FontSize',12)
    title(['MLW - MHW - DF : ' char(cellstr(areaname(i,:))) ' (' num2str(areacode(i)) '), transect '...
        num2str(id(i) - areacode(i)*1000000)],'FontSize',12,'FontWeight','bold')
    leghandle = legend([L; H; D],[' MLW'; ' MHW'; ' DF ']); set(leghandle,'Position',[0.155 0.81 0.0971 0.0877]);
% 	if any(id_old == id(i)); leghandle = legend([L; H; D; Lo; Ho; Do;],[' MLW    '; ' MHW    '; ' DF     '; ' MLW old'; ' MHW old'; ' DF old ']); set(leghandle,'Position',[0.172 0.756 0.0971 0.0877]); end
    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',[WLfolder filesep char(cellstr(areaname(i,:))) '_' num2str(id(i))])

    close all
    clear('hvol_st','hvol_on','hvol_dn','hvol_ot','L','H','D','Ho','Lo','Do');

end
