% convertBKL_TKL_TND2nc
%
% Generates BKL_TKL.nc which contains the following coastal indicators: 
% BKL   (Basis Kust Lijn)
% TKL   (te Toetsen Kustlijn)
% TND   (Trend Kustlijn)
%                              
%
% BKL, TKL, TND are taken from Deltares repository (.xls format) @
%               https://repos.deltares.nl/repos/OpenEarthRawData/trunk/rijkswaterstaat/kustlijnkaartenboek/raw/
%
%
% The data for the TKL-position in all years, except the year 1990, is the
% formally assigned position (source: Yearly Kustlijnkaartenboek - Rijkswaterstaat).
% The data for the TKL-position in 1990 deviates: for 1990 the data from the 
% nota GWWS-91.006 (RKB 91-19) DE BASISKUSTLIJN een technische/morfologische uitwerking is used. 
% It concerns the original extrapolation of the (lineair) MKL trend over 
% the period 1980-1989 to the year 1990. This TKL-position was used as a basis for the BKL-position, 
% but for several location the position was shifted (verlegd) or released (losgelaten). 
% Shifting or releasing was mainly based on morphological arguments (the applied arguments 
% can be found in the previous reference and DGW-93.035 De Basiskustlijn, Norm voor Dynamisch Handhaven). 
% The direction and magnitude of shifting can be obtained by comparing the TKL1990 position to the BKL-position. 
%
%
%	---------------------------------------
% 
%       Giorgio Santinelli
%
%       Giorgio.Santinelli@deltares.nl
%
%       Deltares
%       P.O. Box 177
%       2600 MH Delft
%       The Netherlands
%
%	---------------------------------------
%
% $Id$
% $Date$
% $Author$
% $Revision$
% $HeadURL$

% clear

%% settings
% BKL, TKL, TND
%OPT.ncpath1      = 'https://repos.deltares.nl/repos/OpenEarthRawData/trunk/rijkswaterstaat/kustlijnkaartenboek/raw/';
OPT.ncpath1      = '..\raw\'; % local
OPT.Filepart1 = dir([OPT.ncpath1 'tkl' '*.xls']);
OPT.Filename1 = {OPT.Filepart1.name};

% new files about BKL, TKL and TND
OPT.ncpath2      = '..\raw\TKL2023\'; % local
OPT.Filepart2 = dir([OPT.ncpath2 '*.csv']);
OPT.Filename2 = {OPT.Filepart2.name};
currentyear = str2num(cell2mat(regexp(OPT.ncpath2 , '[0-9]{4}','match'))); %#ok<ST2NM> 

% jarkus database for transect info
%OPT.ncpath3      = 'd:\JARKUSdata\';
%OPT.ncpath3      = 'C:\Users\santinel\temp\openearthrawdata\jarkus\scripts\';
OPT.ncpath3       = 'C:\Users\fuentesm\Marilu\Deltares\Projects\Coastviewer\jarkus\scripts\';
OPT.Filename3    = 'transect_r20240405.nc';

% general settings
OPT.fillvalue    = NaN;
OPT.refdatenum   = datenum(1970,1,1);
OPT.outputfile   = 'BKL_TKL_TND.nc';
OPT.disp         = true;
OPT.dump         = true;

%% BKL, TKL, TND from Deltares repository
% Read data form jarkus transect.nc
id       = ncread([OPT.ncpath3 OPT.Filename3],'id');
cross    = ncread([OPT.ncpath3 OPT.Filename3],'cross_shore');
along    = ncread([OPT.ncpath3 OPT.Filename3],'alongshore');
areaname = ncread([OPT.ncpath3 OPT.Filename3],'areaname');
areacode = ncread([OPT.ncpath3 OPT.Filename3],'areacode');
time     = ncread([OPT.ncpath3 OPT.Filename3],'time');
time_bounds = ncread([OPT.ncpath3 OPT.Filename3],'time_bounds');

jarkusyear = str2num(datestr(OPT.refdatenum+ncread([OPT.ncpath3 OPT.Filename3],'time'),'yyyy')); %#ok<ST2NM>

x        = ncread([OPT.ncpath3 OPT.Filename3],'x');
y        = ncread([OPT.ncpath3 OPT.Filename3],'y');
rsp_x    = ncread([OPT.ncpath3 OPT.Filename3],'rsp_x');
rsp_y    = ncread([OPT.ncpath3 OPT.Filename3],'rsp_y');

% make time array 5 year longer than the jarkus database, as discussed in
% Feb 2022 with RWS
time(end+1:end+5) = [datenum(['01-Jul-' num2str(jarkusyear(end)+1)],'dd-mmm-yyyy') 
                     datenum(['01-Jul-' num2str(jarkusyear(end)+2)],'dd-mmm-yyyy')
                     datenum(['01-Jul-' num2str(jarkusyear(end)+3)],'dd-mmm-yyyy')
                     datenum(['01-Jul-' num2str(jarkusyear(end)+4)],'dd-mmm-yyyy')
                     datenum(['01-Jul-' num2str(jarkusyear(end)+5)],'dd-mmm-yyyy')] - OPT.refdatenum;
% %% to allow for a year with jarkus transect.nc database not yet available [if necessary !!]

%% Read BKL, (periode), TKL, TND

%BKL
% cycle to read each areas BKL from OLD files
for j = 1:length(OPT.Filename1)
    xlsTKL    = xlsread([OPT.ncpath1, char(OPT.Filename1(j))]);
    xlsBKL(j)  = {(xlsTKL(~isnan(xlsTKL(:,1)),3))'}; % transpose: 2Dmatrix
    xlsTKLy(j) = {(xlsTKL(~isnan(xlsTKL(:,1)),4))'}; % //
    xlsTND(j)  = {(xlsTKL(~isnan(xlsTKL(:,1)),6))'}; % //
    
    xlscode(j) = {xlsTKL(~isnan(xlsTKL(:,1)),1)};
    xlsraai(j) = {xlsTKL(~isnan(xlsTKL(:,1)),2)};
    
    name = char(OPT.Filename1(j));
    year = str2double(name(4:7));
    date_TKL(j) = datenum(['01-Jul-' num2str(year)],'yyyy') - OPT.refdatenum;
    time_TKL(j) = year;
end

% read each areas BKL from new files
% BKL-nota
xlsBKLoverzicht = xlsread([OPT.ncpath1 'BKL-overzicht incl 1990.xls']);
xlscodenota     = {xlsBKLoverzicht(~isnan(xlsBKLoverzicht(:,1)),1)};
xlsraainota     = {xlsBKLoverzicht(~isnan(xlsBKLoverzicht(:,1)),2)};
xlsBKLnota      = {xlsBKLoverzicht(~isnan(xlsBKLoverzicht(:,1)),3)};

%TKL, TND
% now add the 5years streak of TKL and Trend from last year
%find columns. This would probably have to be updated in 2023.
colTrend = []; colTkl = []; colKustvaK = []; colRaai = []; colBkl = [];
xlsTrendfuture = []; xlsTKLfuture = []; xlsKustvaKfuture = []; xlsRaaifuture = []; xlsBKLfuture = [];
for j = 1:length(OPT.Filename2)
    [xlsTKLfuture_all_num, xlsTKLfuture_all_txt, xlsTKLfuture_all_raw]    = xlsread([OPT.ncpath2, char(OPT.Filename2{j})]);
    for jtrend = 1:size(xlsTKLfuture_all_txt,2)
        cKustvaK = strfind(char(xlsTKLfuture_all_txt(1,jtrend)), 'Kustvak');
        if length(cKustvaK) ~= 0
            colKustvaK(length(colKustvaK)+1) = jtrend;
        end
        cRaai = strfind(char(xlsTKLfuture_all_txt(1,jtrend)), 'Raai');
        if length(cRaai) ~= 0
            colRaai(length(colRaai)+1) = jtrend;
        end        
        cTrend = strfind(char(xlsTKLfuture_all_txt(1,jtrend)), 'Trend [m/jaar]');
        if length(cTrend) ~= 0
            colTrend(length(colTrend)+1) = jtrend;
        end
        cTkl = strfind(char(xlsTKLfuture_all_txt(1,jtrend)), 'TKL [m + RSP]');
        if length(cTkl) ~= 0
            colTkl(length(colTkl)+1) = jtrend;
        end
        cBkl = regexp(char(xlsTKLfuture_all_txt(1,jtrend)), '^BKL$');
        if length(cBkl) ~= 0
            colBkl(length(colBkl)+1) = jtrend;
        end
    end
    colTrend_u = unique(colTrend);
    colTkl_u = unique(colTkl);
    colKustvaK_u = unique(colKustvaK);
    colRaai_u = unique(colRaai);
    colBkl_u = unique(colBkl);

    xlsTrendfuture   = [xlsTrendfuture;   cell2mat(nanify(xlsTKLfuture_all_raw(6:end,colTrend_u)))]; 
    xlsTKLfuture     = [xlsTKLfuture;     cell2mat(nanify(xlsTKLfuture_all_raw(6:end,colTkl_u)))]; 
    xlsKustvaKfuture = [xlsKustvaKfuture; cell2mat(nanify(xlsTKLfuture_all_raw(6:end,colKustvaK_u)))]; 
    xlsRaaifuture    = [xlsRaaifuture;    cell2mat(nanify(xlsTKLfuture_all_raw(6:end,colRaai_u)))]; 
    xlsBKLfuture     = [xlsBKLfuture;     cell2mat(nanify(xlsTKLfuture_all_raw(6:end,colBkl_u)))]; 
    
    %xlsTrendfuture = [xlsTrendfuture; str2num(char(xlsTKLfuture_all_raw(6:end,colTrend_u)))]; %#ok<ST2NM> 
    %xlsTKLfuture = [xlsTKLfuture; cell2mat(cellfun( @str2double, xlsTKLfuture_all_raw(6:end,colTkl_u), 'UniformOutput', false))];   
    %xlsKustvaKfuture = [xlsKustvaKfuture; cell2mat(xlsTKLfuture_all_raw(6:end,colKustvaK_u))]; 
    %xlsRaaifuture = [xlsRaaifuture; str2num(char(xlsTKLfuture_all_raw(6:end,colRaai_u)))]; %#ok<ST2NM> 
    %xlsBKLfuture = [xlsBKLfuture; cell2mat(cellfun( @str2double, xlsTKLfuture_all_raw(6:end,colBkl_u), 'UniformOutput', false))]; 
end

%% create the matrix for each parameter, with a structure similar to the jarkus file
% define the final matrix. variables definition
Xbkl = NaN(length(time),length(id));     Xtkl = NaN(length(time),length(id)); 
tnd  = NaN(length(id),1);  yearTKL = NaN(length(time),length(id));

[~, dumt] = ismember(time_TKL, jarkusyear);
for tx = 1:length(time_TKL)
    if (xlscode{tx}*1e6 + xlsraai{tx}) == 13001528; xlsraai{tx}(find(xlscode{tx}*1e6 + xlsraai{tx} == 13001528)) = 1525; end % issues with missing transects: 13001528, 17001242
    if (xlscode{tx}*1e6 + xlsraai{tx}) == 17001242; xlsraai{tx}(find(xlscode{tx}*1e6 + xlsraai{tx} == 17001242)) = 1241; end
    [dum_raai, dum_id] = ismember(id,xlscode{tx}*1e6 + xlsraai{tx});
    for ix = 1:length(id)
        if dum_raai(ix)
            Xbkl(dumt(tx),ix) = xlsBKL{tx}(dum_id(ix));
            % not too correct to insert them, cause boundaries for calculations change 
            %Xtkl(dumt(tx),ix) = xlsTKLy{tx}(dum_id(ix)); 
            %tnd(dumt(tx),ix)  = xlsTND{tx}(dum_id(ix));
        end
    end
end

% do the same for latest year(s), including predictions for the future tnd and TKL
yahead = 5;
dumtnew = currentyear:currentyear+yahead;
for tx = 1:length(dumtnew)
    if (xlsKustvaKfuture*1e6 + xlsRaaifuture) == 13001528; xlsRaaifuture(find(xlsKustvaKfuture*1e6 + xlsRaaifuture == 13001528)) = 1525; end % issues with missing transects: 13001528, 17001242
    if (xlsKustvaKfuture*1e6 + xlsRaaifuture) == 17001242; xlsRaaifuture(find(xlsKustvaKfuture*1e6 + xlsRaaifuture == 17001242)) = 1241; end
    [dum_raai, dum_id] = ismember(id,xlsKustvaKfuture*1e6 + xlsRaaifuture);
    for ix = 1:length(id)
        if dum_raai(ix)
            Xbkl(end-yahead,ix) = xlsBKLfuture(dum_id(ix)); % xbkl confirmed in current year, and kept blank for future years
            tnd(ix)  = xlsTrendfuture(dum_id(ix));
            if tx > 1 % skip the current year
                
                Xtkl(end-yahead+tx-1:end,ix) = xlsTKLfuture(dum_id(ix),tx-1)'; % valid for the last five years of the series
                % Xtkl(end-yahead+tx:end,ix) =
                % xlsTKLfuture(dum_id(ix),tx-1)'; % valid for the last five
                % years of the series % without -1 if jarkus of that year is
                % there already
            end
        end
    end
end


%% insert Bkl - nota
[~, dumt_nota] = ismember(1990, jarkusyear);
if (xlscodenota{1}*1e6 + xlsraainota{1}) == 13001528; xlsraainota(find(xlscodenota{1}*1e6 + xlsraainota{1} == 13001528)) = 1525; end % issues with missing transects: 13001528, 17001242
if (xlscodenota{1}*1e6 + xlsraainota{1}) == 17001242; xlsraainota(find(xlscodenota{1}*1e6 + xlsraainota{1} == 17001242)) = 1241; end
[dum_raai dum_id] = ismember(id, xlscodenota{1}*1e6 + xlsraainota{1});
for ix = 1:length(id) 
	if dum_raai(ix); Xbkl(dumt_nota,ix) = xlsBKLnota{1}(dum_id(ix)); end;
end

%% convert bkl, tkl in RD coordinates and in WGS84
% BKL, TKL in RD coordinates
BKL_x = NaN(length(time),length(id)); BKL_y = NaN(length(time),length(id));
TKL_x = NaN(length(time),length(id)); TKL_y = NaN(length(time),length(id));

for i = 1:length(id)
    disp([num2str(i) '/' num2str(length(id))])
    for t = 1:length(time)
        BKL_x(t,i) = interp1(cross, x(:,i), Xbkl(t,i));
        BKL_y(t,i) = interp1(cross, y(:,i), Xbkl(t,i));
        TKL_x(t,i) = interp1(cross, x(:,i), Xtkl(t,i));
        TKL_y(t,i) = interp1(cross, y(:,i), Xtkl(t,i));
        
    end
end

% BKL, TKL in WGS84 coordinates
[BKL_lon, BKL_lat] = convertCoordinates(BKL_x,BKL_y,'CS1.code',28992,'CS2.code',4326);
[TKL_lon, TKL_lat] = convertCoordinates(TKL_x,TKL_y,'CS1.code',28992,'CS2.code',4326);

%% create structure BKL_TKL_TND
s.BKL_lat   = BKL_lat; 
s.BKL_lon   = BKL_lon;
s.BKL_x     = BKL_x; 
s.BKL_y     = BKL_y; 
s.Xbkl      = Xbkl;
s.TKL_lat   = TKL_lat; 
s.TKL_lon   = TKL_lon;
s.TKL_x     = TKL_x; 
s.TKL_y     = TKL_y; 
s.Xtkl      = Xtkl;
s.tnd       = tnd; 
% include relevant information from Jarkus
s.id   = id;   
s.cross = cross; 
s.along = along;
s.time        = time; 
s.time_bounds = time_bounds;
s.areaname = areaname; 
s.areacode = areacode;
% add for convenience
s.jarkusyear = jarkusyear; 


%% create BKL_TKL_TND
bkl2netcdf(OPT.outputfile, s,'processing_level', 'final');


function output = nanify(input)
% This function replaces string values of 'NaN' in a cell array with the numeric value NaN.
output = input;
% Loop through each element of the input.
for ii = 1:numel(input)
    % Check if the element is a string with the value 'NaN'.
    if ischar(input{ii}) && strcmp(input{ii}, 'NaN')
        % If it is, replace it with the numeric value NaN.
        output{ii} = NaN;
    end
end
end
