% BKL_TKL_TND2kml
%
% Just like the Kustlijnkaartenboek in Google Earth
%
% Generates BKL_TKL_TND.kml which shows the following coastal indicators in GE: 
% BKL   (Basis Kust Lijn)
% TKL   (te Toetsen Kustlijn)
% TND   (Trend Kustlijn)
%                              
% in the standard format (from OPENDAP)
% BKL, TKL, TND are taken from opendap.deltares.n (.nc format) @
%               http://opendap.deltares.nl/thredds/catalog/opendap/rijkswat
%               erstaat/BKL_TKL_MKL/catalog.html
%
% raw data (from RWS)
% BKL, TKL, TND are taken from Deltares repository (.xls format) @
%               https://repos.deltares.nl/repos/OpenEarthRawData/trunk/rijk
%               swaterstaat/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

% BKL_TKL_TND.nc source
OPT.ncpath      = 'http://opendap.deltares.nl/thredds/dodsC/opendap/rijkswaterstaat/BKL_TKL_MKL/';
OPT.ncpath      = 'd:\checkouts\kustlijnkaartenboek\scripts\';
OPT.Filename    = 'BKL_TKL_TND.nc';

% transects angle source
OPT.ncpath1      = 'http://opendap.deltares.nl/thredds/dodsC/opendap/rijkswaterstaat/jarkus/profiles/';
% OPT.ncpath1      = 'F:\JARKUSdata\';
OPT.Filename1    = 'transect.nc';

% .kml storage
OPT.kmlpath      = 'd:\checkouts\kustlijnkaartenboek\kml\'; %'\\kml.deltares.nl\kml\kml\rijkswaterstaat\kustlijnkaartenboek\'; %'//kml.deltares.nl/kml/rijkswaterstaat/kustlijnkaartenboek/';

% additional options
OPT.z     = 'clampToGround'; % 50;
TNDfactor = 25;

%% Read data form BKL_TKL_TND.nc
id       = nc_varget([OPT.ncpath OPT.Filename],'id');
areaname = nc_varget([OPT.ncpath OPT.Filename],'areaname');
areacode = nc_varget([OPT.ncpath OPT.Filename],'areacode');
time     = nc_cf_time([OPT.ncpath OPT.Filename],'time');

angle    = nc_varget([OPT.ncpath1 OPT.Filename1],'angle');

BKL      = nc_varget([OPT.ncpath OPT.Filename],'basal_coastline');
TKL      = nc_varget([OPT.ncpath OPT.Filename],'testing_coastline');
TND      = nc_varget([OPT.ncpath OPT.Filename],'trend');

BKL_lat  = nc_varget([OPT.ncpath OPT.Filename],'basal_coastline_lat');
BKL_lon  = nc_varget([OPT.ncpath OPT.Filename],'basal_coastline_lon');
TKL_lat  = nc_varget([OPT.ncpath OPT.Filename],'testing_coastline_lat');
TKL_lon  = nc_varget([OPT.ncpath OPT.Filename],'testing_coastline_lon');

BKL_x  = nc_varget([OPT.ncpath OPT.Filename],'basal_coastline_x');
BKL_y  = nc_varget([OPT.ncpath OPT.Filename],'basal_coastline_y');
TKL_x  = nc_varget([OPT.ncpath OPT.Filename],'testing_coastline_x');
TKL_y  = nc_varget([OPT.ncpath OPT.Filename],'testing_coastline_y');

jarkusyear = str2num(datestr(nc_cf_time([OPT.ncpath OPT.Filename]),'yyyy')); %#ok<ST2NM>

timeIn   = datenum(num2str(jarkusyear),'yyyy');
timeOut  = datenum(num2str(jarkusyear+1),'yyyy'); 

%time: starting from 1990, j = 26
yearID   = find(jarkusyear==1990);
timeIn   = timeIn(26:end);
timeOut  = timeOut(26:end);

[a1, a2]   = unique(areacode);
a3        = a2 + 1;
a3(2:end) = a3(1:end-1);
a3(1)     = 1;
% add Z-Vlaanderen
a3(end+1) = a2(end)+1;
a2(end+1) = 2268;

% improved...
a_ini = a2;
a_fin = a_ini-1;
a_fin(1:end-1) = a_fin(2:end);
a_fin(end) = size(id,1);

%% BKL
for i = 1:length(a1)+1
%     % File changes to make it more readable: remove a certain number of
%     % NaN's between notNaN's

    for jj = yearID:size(BKL_lat,1) %time: starting from 1990, j = 26
        for ii = 1:size(BKL_lat(:,a3(i):a2(i)),2)
            iii  = ii + 2;
            if iii-1 < a2(i) - a3(i) + 1
                if isnan(BKL_lat(jj,a3(i)+ii-1)) && ~isnan(BKL_lat(jj,a3(i)+ii))
                    BKL_lat(jj,a3(i)+ii-1)   = BKL_lat(jj,a3(i)+ii);
                    BKL_lon(jj,a3(i)+ii-1)   = BKL_lon(jj,a3(i)+ii);
                elseif isnan(BKL_lat(jj,a3(i)+ii-1)) && isnan(BKL_lat(jj,a3(i)+ii)) && ~isnan(BKL_lat(jj,a3(i)+iii-1))
                    BKL_lat(jj,a3(i)+ii-1)   = BKL_lat(jj,a3(i)+iii-1);
                    BKL_lat(jj,a3(i)+ii) = BKL_lat(jj,a3(i)+iii-1);
                    BKL_lon(jj,a3(i)+ii-1)   = BKL_lon(jj,a3(i)+iii-1);
                    BKL_lon(jj,a3(i)+ii) = BKL_lon(jj,a3(i)+iii-1);
                end
            end
        end
    end


    % NaN to not defined BKL
    % Scheveningen harbour in Delfland
    BKL_lat(:,id==9010193) = NaN(size(BKL_lat,1),1); BKL_lon(:,id==9010193) = NaN(size(BKL_lon,1),1);
    % harbour in Breskens harbour Zeeuws-Vlaanderen
    BKL_lat(:,id==17000138) = NaN(size(BKL_lat,1),1); BKL_lon(:,id==17000138) = NaN(size(BKL_lon,1),1);
    
    if i > 1
        if ~isnumeric(OPT.z)
            KMLline(BKL_lat(yearID:end,a3(i):a2(i))',BKL_lon(yearID:end,a3(i):a2(i))',... % 26 is year1990
                OPT.z,...
                'fileName',[OPT.kmlpath, 'BKLline_', num2str(a1(i-1)),'.kml'],...
                'kmlName',['BKL: ' char(cellstr(areaname(a2(i-1),:)))],...
                'timeIn',timeIn,...
                'timeOut',timeOut,...
                'lineWidth',3,...
                'lineColor',[0.4 0 0.4],...
                'lineAlpha',.7,...
                'fillColor',[0.4 0 0.4],...
                'fillAlpha',.7);        
        else
            KMLline(BKL_lat(yearID:end,a3(i):a2(i))',BKL_lon(yearID:end,a3(i):a2(i))',... % 26 is year1990
                OPT.z*ones(size(BKL_lat(yearID:end,a3(i):a2(i)),2),size(BKL_lat(yearID:end,a3(i):a2(i)),1)),...
                'fileName',[OPT.kmlpath, 'BKLline_', num2str(a1(i-1)),'.kml'],...
                'kmlName',['BKL: ' char(cellstr(areaname(a2(i-1),:)))],...
                'timeIn',timeIn,...
                'timeOut',timeOut,...
                'lineWidth',3,...
                'lineColor',[0.4 0 0.4],...
                'lineAlpha',.7,...
                'fillColor',[0.4 0 0.4],...
                'fillAlpha',.7);
        end
    end
        
end

BKLname = [];
for i = 1:length(a1)
    BKLname = [BKLname; {[OPT.kmlpath, 'BKLline_', num2str(a1(i)),'.kml']}];
end

KMLmerge_files('fileName',[OPT.kmlpath 'BKLline.kml'],...
            'description','BKL along the Dutch coast',...
            'sourceFiles',BKLname,...
                'kmlName','BKL',...
      'deleteSourceFiles',0);

%% TKL - BKL

% get TND lat,lon on the map
TND_x   = BKL_x + TNDfactor*TND.*repmat((sin(angle/180*pi))',length(time),1);
TND_y   = BKL_y + TNDfactor*TND.*repmat((cos(angle/180*pi))',length(time),1);

[TND_lon, TND_lat] = convertCoordinates(TND_x,TND_y,'CS1.code',28992,'CS2.code',4326);
    
for i = 1:length(a1)+1
    for k = 1:(a2(i) - a3(i) + 1)
        Color = ones(size(TKL,1)-yearID+1,3);
        for j = yearID:size(TKL,1) %time: starting from 1990, j = 26
            % TKL - BKL
            buffer = TKL(j,k + a3(i) - 1) - BKL(j,k + a3(i) - 1);
            if isnan(buffer); TND_lat(j,k + a3(i) - 1) = NaN; TND_lon(j,k + a3(i) - 1) = NaN; end % allows for missing TKL but present TND
            if buffer >= 0
                if TND(j,k + a3(i) - 1) >= 0
                    Color(j-yearID+1,:) = [0 0.6 0];
                elseif TND(j,k + a3(i) - 1) < 0
                    Color(j-yearID+1,:) = [0.6 1 0.6];
                end
            elseif buffer < 0
                if TND(j,k + a3(i) - 1) >= 0
                    Color(j-yearID+1,:) = [1 0.6 0.6];
                elseif TND(j,k + a3(i) - 1) < 0
                    Color(j-yearID+1,:) = [1 0 0];
                end        
            end
        end
% % line with TKL in the absolute position
% if i > 1
%         KMLline([TKL_lat(yearID:end,k + a3(i) - 1) BKL_lat(yearID:end,k + a3(i) - 1)]',...
%             [TKL_lon(yearID:end,k + a3(i) - 1) BKL_lon(yearID:end,k + a3(i) - 1)]',...
%             OPT.z*ones(2,size(TKL_lat,1)),...
%              'fileName',['TKL_BKL', num2str(id(k + a3(i) - 1)),'.kml'],...
%               'kmlName',[char(cellstr(areaname(a2(i-1),:))) '_' num2str(id(k + a3(i) - 1))],...
%                'timeIn',timeIn,...
%               'timeOut',timeOut,...
%             'lineWidth',3,...
%             'lineColor',Color,...
%             'lineAlpha',.7,...
%             'fillColor',Color,...
%             'fillAlpha',1);
% end

    % % line with trend value, as in the Kustlijnkaartenboek
    if i > 1
        if ~isnumeric(OPT.z)
            KMLline([TND_lat(yearID:end,k + a3(i) - 1) BKL_lat(yearID:end,k + a3(i) - 1)]',... % 26 is year1990
                [TND_lon(yearID:end,k + a3(i) - 1) BKL_lon(yearID:end,k + a3(i) - 1)]',...
                OPT.z,...
                'fileName',[OPT.kmlpath, 'TKL_BKL', num2str(id(k + a3(i) - 1)),'.kml'],...
                'kmlName',[char(cellstr(areaname(a2(i-1),:))) '_' num2str(id(k + a3(i) - 1))],...
                'timeIn',timeIn,...
                'timeOut',timeOut,...
                'lineWidth',3,...
                'lineColor',Color,...
                'lineAlpha',.7,...
                'fillColor',Color,...
                'fillAlpha',1);
        else
            KMLline([TND_lat(yearID:end,k + a3(i) - 1) BKL_lat(yearID:end,k + a3(i) - 1)]',... % 26 is year1990
                [TND_lon(yearID:end,k + a3(i) - 1) BKL_lon(yearID:end,k + a3(i) - 1)]',...
                OPT.z*ones(2,size(TKL_lat,1)),...
                'fileName',[OPT.kmlpath, 'TKL_BKL', num2str(id(k + a3(i) - 1)),'.kml'],...
                'kmlName',[char(cellstr(areaname(a2(i-1),:))) '_' num2str(id(k + a3(i) - 1))],...
                'timeIn',timeIn,...
                'timeOut',timeOut,...
                'lineWidth',3,...
                'lineColor',Color,...
                'lineAlpha',.7,...
                'fillColor',Color,...
                'fillAlpha',1);
        end
    end

    end
    
    TKLnamearea = [];
    for k = 1:(a2(i) - a3(i) + 1)
        TKLnamearea = [TKLnamearea; {[OPT.kmlpath, 'TKL_BKL', num2str(id(k + a3(i) - 1)),'.kml']}];
    end    
    if i > 1
        KMLmerge_files('fileName',[OPT.kmlpath, 'TKLline_'  num2str(a1(i-1)) '.kml'],...
                'sourceFiles',TKLnamearea,...
                    'kmlName',['TKL and TND: ' char(cellstr(areaname(a2(i-1),:)))],...
          'deleteSourceFiles',1);
     end

end

TKLname = [];
for i = 1:length(a1)
    TKLname = [TKLname; {[OPT.kmlpath, 'TKLline_', num2str(a1(i)),'.kml']}];
end

KMLmerge_files('fileName',[OPT.kmlpath, 'TKLline.kml'],...
            'description','TKL along the Dutch coast',...
            'sourceFiles',TKLname,...
                'kmlName','TKL',...
      'deleteSourceFiles',0);
  
% merge TKL and BKL by area
for i = 1:length(a1)
    BKLTKLarea = cellstr([OPT.kmlpath, 'BKLline_' num2str(a1(i)) '.kml'; OPT.kmlpath, 'TKLline_' num2str(a1(i)) '.kml']);
    areanamekml = strrep(char(cellstr(areaname(a2(i),:))),' ','');
    areanamekml = strrep(areanamekml,'/','_');
    
    KMLmerge_files('fileName',[OPT.kmlpath 'Kustvak' num2str(a1(i)) '_' areanamekml '_noLogo' '.kml'],...
                'description',['BKL and TKL in ' areanamekml],...
                'sourceFiles',BKLTKLarea,...
                    'kmlName',['BKL, TKL and TND: ' areanamekml],...
          'deleteSourceFiles',1);

    KMLmerge_files('fileName',[OPT.kmlpath 'Kustvak' num2str(a1(i)) '_' areanamekml '.kml'],...
                'description',['BKL and TKL in ' areanamekml],...
                'sourceFiles',{[OPT.kmlpath 'Kustvak' num2str(a1(i)) '_' areanamekml '_noLogo' '.kml']; [OPT.kmlpath, 'logos.kml']},...
                    'kmlName',['BKL, TKL and TND: ' areanamekml],...
          'deleteSourceFiles',0);
      
      delete([OPT.kmlpath 'Kustvak' num2str(a1(i)) '_' areanamekml '_noLogo' '.kml'])
  
end



% merge TKL and BKL total  
KMLmerge_files('fileName',[OPT.kmlpath 'Kustlijnkaartenboek.kml'],...
            'description','BKL and TKL along the Dutch coast',...
            'sourceFiles',{[OPT.kmlpath 'BKLline.kml']; [OPT.kmlpath, 'TKLline.kml']; [OPT.kmlpath, 'logos.kml']},...
                'kmlName','BKL, TKL and TND',...
      'deleteSourceFiles',0);

% Dutch version _nl
KML2kmz([OPT.kmlpath 'Kustlijnkaartenboek.kml'], ...
    [OPT.kmlpath 'logos' filesep 'oet4GE.png'], ...
    [OPT.kmlpath 'logos' filesep 'RWS_transparent_small.png'], ...
    [OPT.kmlpath 'logos' filesep 'Leg_nl.png'])

% English version _en
% not yet available...
  
  
  