% convertDF2nc
%
% Generates Dune Foot position NetCDF file in
% http://opendap.deltares.nl/thredds/dodsC/opendap/rijkswaterstaat/DuneFoot/DF.nc
%
% DF.nc contains:
% Dune Foot(1) defined as the crossing between the upper horizontal boundary
% (not always at +3.00m nap) for the calculation of the MKL and the coastline profile:
% Dune Foot position and related variables following this definition are set as DF_upperMKL*
% 
% Dune Foot(2) defined as the crossing between the water level at +3.00m NAP 
% and the JarKus coastline profile. 
% Dune Foot position and related variables following this definition are set as DF_threeNAP*
% 
% NB. 
% a) Dune Foot is taken from STRANDLIJNEN dataset (1843 - 1998), already in .nc format @
% https://repos.deltares.nl/repos/OpenEarthRawData/trunk/rijkswaterstaat/strandlijnen/processed/strandlijnen.nc
% b) Dune Foot is calculated using JarKus dataset (1965 - ...) @
% http://opendap.deltares.nl/thredds/dodsC/opendap/rijkswaterstaat/jarkus/profiles/transect.nc
% c) STRANDLIJNEN dataset is now MERGED WITH JarKus profiles intersections. 
% d) Dune foot in 1843 - 1964 is kept same for both Dune Foot (1 & 2) forementioned definitions
%
% $Id$
% $Date$
% $Author$
% $Revision$
% $HeadURL$
% $Keywords$ 
%
%%

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

stringsize = 100;

% strandlijnen database - dataset 1843 - 1998 (1843-1964 taken into account)
%OPT.ncpath1      = 'c:\Users\santinel\Documents\openearthrawdata\rijkswaterstaat\strandlijnen\';
OPT.ncpath1      = 'http://opendap.deltares.nl/thredds/dodsC/opendap/rijkswaterstaat/strandlijnen/';
OPT.Filename1    = 'strandlijnen.nc';
% JarKus database - dataset 1965 - ...
OPT.ncpath2      = 'http://opendap.deltares.nl/thredds/dodsC/opendap/rijkswaterstaat/jarkus/profiles//';
OPT.Filename2    = 'transect_r20240405.nc';
% % MKL database for the upper horizontal limit
% OPT.ncpath3      =
% 'http://opendap.deltares.nl/thredds/dodsC/opendap/rijkswaterstaat/BKL_TKL_MKL/';%
% % local
OPT.ncpath3      = 'http://opendap.deltares.nl/thredds/dodsC/opendap/rijkswaterstaat/BKL_TKL_MKL/';
OPT.Filename3    = 'MKL.nc';

%% Read data (1/3)
% JarKus database - dataset 1965 - ...
id       = nc_varget([OPT.ncpath2 OPT.Filename2],'id');
timej    = nc_varget([OPT.ncpath2 OPT.Filename2],'time');
timebj   = nc_varget([OPT.ncpath2 OPT.Filename2],'time_bounds');
lat      = nc_varget([OPT.ncpath2 OPT.Filename2],'lat');
lon      = nc_varget([OPT.ncpath2 OPT.Filename2],'lon');
rsp_x    = nc_varget([OPT.ncpath2 OPT.Filename2],'rsp_x');
rsp_y    = nc_varget([OPT.ncpath2 OPT.Filename2],'rsp_y');
areacode = nc_varget([OPT.ncpath2 OPT.Filename2],'areacode');
areaname = nc_varget([OPT.ncpath2 OPT.Filename2],'areaname');
along    = nc_varget([OPT.ncpath2 OPT.Filename2],'alongshore');
cross    = nc_varget([OPT.ncpath2 OPT.Filename2],'cross_shore');

time_year = str2num(datestr(timej+datenum('1970','yyyy'),'yyyy')); %#ok<ST2NM>
transectid = id - areacode*1e6;

%% Read data (2/3)
% % MKL database [for the upper horizontal limit]
id_MKL     = nc_varget([OPT.ncpath3, OPT.Filename3],'id');
time_MKL = nc_varget([OPT.ncpath3, OPT.Filename3],'time');
HibndMKL_oldindex = nc_varget([OPT.ncpath3, OPT.Filename3],'high_boundary_MKL');
ThreeNAP_oldindex = 3*ones(size(HibndMKL_oldindex));
time_year_MKL = str2num(datestr(nc_cf_time([OPT.ncpath3, OPT.Filename3]),'yyyy')); %#ok<ST2NM>

% %% 2.5/3 !!!
% % Fake MKL data for the last year while waiting for MKL_TKL_BKL data
if length(time_year_MKL) < length(time_year)
    addY = length(time_year) - length(time_year_MKL);
    time_MKL(end+addY) = time_MKL(end)+366*addY; % very quick..
    HibndMKL_oldindex(end+addY,:) = NaN(addY,length(id_MKL));
    ThreeNAP_oldindex(end+addY,:) = 3*ones(addY,length(id_MKL));
    time_year_MKL(end+addY) = time_year_MKL(end) + addY;
end


%% Read data (3/3)
% strandlijnen database - dataset 1843 - 1964 [from 1843-1998 strandlijnen.nc file]
id_STRAND     = nc_varget([OPT.ncpath1 OPT.Filename1],'id');
time_STRAND   = nc_varget([OPT.ncpath1 OPT.Filename1],'time'); % year
DF_STRAND     = nc_varget([OPT.ncpath1 OPT.Filename1],'DF');
DF_STRANDx    = nc_varget([OPT.ncpath1 OPT.Filename1],'DF_x');
DF_STRANDy    = nc_varget([OPT.ncpath1 OPT.Filename1],'DF_y');
% LAT and LON are reversed in strandlijnen.nc
DF_STRANDlat  = nc_varget([OPT.ncpath1 OPT.Filename1],'DF_lon');
DF_STRANDlon  = nc_varget([OPT.ncpath1 OPT.Filename1],'DF_lat');

time_STRAND_year  = str2num(datestr(time_STRAND+datenum('1970','yyyy'),'yyyy')); %#ok<ST2NM>
time_total        = (time_STRAND_year(1):1:time_year_MKL(end))';
time_totalstring  = arrayfun(@(x) ['01-Jul-' num2str(x)], time_total, 'UniformOutput', false); % set at 01-Jul-yyyy
time_totaldatenum = datenum(time_totalstring);

time_totalstring_beg  = arrayfun(@(x) ['01-Jan-' num2str(x)], time_total, 'UniformOutput', false); % set at 01-Jan-yyyy
time_totalstring_end  = arrayfun(@(x) ['31-Dec-' num2str(x)], time_total, 'UniformOutput', false); % set at 31-Dec-yyyy
time_totaldatenum_beg = datenum(time_totalstring_beg);
time_totaldatenum_end = datenum(time_totalstring_end);

[idID_STa, idID_STb]       = ismember(id,id_STRAND); % 324 common ones of 334 from strandlijnen.nc
[timeIDa, timeIDb]   = ismember(time_year_MKL, time_STRAND_year); % 34 common years 

% leave out the part after 1964 (1964 included)
DF_STRAND(:,timeIDb(1):end) = [];
DF_STRANDx(:,timeIDb(1):end) = [];
DF_STRANDy(:,timeIDb(1):end) = [];
DF_STRANDlat(:,timeIDb(1):end) = [];
DF_STRANDlon(:,timeIDb(1):end) = [];

% Preallocation and definition
DF_upperMKL_tot     = NaN(length(id),length(time_totaldatenum));
DF_upperMKL_x_tot   = NaN(length(id),length(time_totaldatenum));
DF_upperMKL_y_tot   = NaN(length(id),length(time_totaldatenum));
DF_upperMKL_lat_tot = NaN(length(id),length(time_totaldatenum));
DF_upperMKL_lon_tot = NaN(length(id),length(time_totaldatenum));
DF_upperMKL         = NaN(length(id),length(time_totaldatenum));

DF_threeNAP_tot     = NaN(length(id),length(time_totaldatenum));
DF_threeNAP_x_tot   = NaN(length(id),length(time_totaldatenum));
DF_threeNAP_y_tot   = NaN(length(id),length(time_totaldatenum));
DF_threeNAP_lat_tot = NaN(length(id),length(time_totaldatenum));
DF_threeNAP_lon_tot = NaN(length(id),length(time_totaldatenum));
DF_threeNAP         = NaN(length(id),length(time_totaldatenum));

% fill the variables with old values from strandlijnen.nc
for ii = 1:length(id)
        if idID_STa(ii)
            DF_upperMKL_tot(ii,1:(timeIDb(1)-1))     = DF_STRAND(idID_STb(ii),:);
            DF_upperMKL_x_tot(ii,1:(timeIDb(1)-1))   = DF_STRANDx(idID_STb(ii),:);
            DF_upperMKL_y_tot(ii,1:(timeIDb(1)-1))   = DF_STRANDy(idID_STb(ii),:);
            DF_upperMKL_lat_tot(ii,1:(timeIDb(1)-1)) = DF_STRANDlat(idID_STb(ii),:);
            DF_upperMKL_lon_tot(ii,1:(timeIDb(1)-1)) = DF_STRANDlon(idID_STb(ii),:);
            
            DF_threeNAP_tot(ii,1:(timeIDb(1)-1))     = DF_STRAND(idID_STb(ii),:);
            DF_threeNAP_x_tot(ii,1:(timeIDb(1)-1))   = DF_STRANDx(idID_STb(ii),:);
            DF_threeNAP_y_tot(ii,1:(timeIDb(1)-1))   = DF_STRANDy(idID_STb(ii),:);
            DF_threeNAP_lat_tot(ii,1:(timeIDb(1)-1)) = DF_STRANDlat(idID_STb(ii),:);
            DF_threeNAP_lon_tot(ii,1:(timeIDb(1)-1)) = DF_STRANDlon(idID_STb(ii),:);
        end
end

% do the same for elevation boundaries in the MKL dataset
[idID_MKLa, idID_MKLb]       = ismember(id,id_MKL);

HibndMKL     = NaN(length(time_MKL),length(id));
ThreeNAP     = NaN(length(time_MKL),length(id));

% fill the variables
for ii = 1:length(id)
        if idID_MKLa(ii)
            HibndMKL(:,ii)     = HibndMKL_oldindex(:,idID_MKLb(ii));
            ThreeNAP(:,ii)     = ThreeNAP_oldindex(:,idID_MKLb(ii));
        end
end


%% Evaluation of DF
% counters
[kv, m, dummy] = unique(areacode);
% n          = vertcat(1, m(1:end-1)+1);
n          = vertcat(m(2:end)-1,length(id));

% Dune Foot (1) definition
% DuneFootm is the same as 'high_boundary_MKL' in MKL.nc
dune_foot_upperMKL_cross = NaN(length(areacode),length(time_MKL));
dune_foot_upperMKL_lat   = NaN(length(areacode),length(time_MKL)); 
dune_foot_upperMKL_lon   = NaN(length(areacode),length(time_MKL));

% Dune Foot (2) definition
% most seaward intersection with +3m NAP
dune_foot_threeNAP_cross = NaN(length(areacode),length(time_MKL));
dune_foot_threeNAP_lat   = NaN(length(areacode),length(time_MKL)); 
dune_foot_threeNAP_lon   = NaN(length(areacode),length(time_MKL));

%% dune_foot_upperMKL: intersection with the upper boundary of MKL definition
for i = 1:length(kv) % loop for area
    num_tr  = n(i)-m(i)+1;

    for j = 1:num_tr % loop for transect
        altitude = squeeze(nc_varget([OPT.ncpath2 OPT.Filename2],'altitude',[0 m(i)+j-2 0],[-1 1 -1]));
        disp(['DF1: transect ' num2str(j) ' of area ' num2str(kv(i)) ...
            ' (' num2str(m(i)+j-1) '/' num2str(length(id)) ')']);
        
        for k = 1:length(time_MKL) % loop for year (...in MKL, because the upper bnd is required)
            if ~isempty(altitude(k,~isnan(altitude(k,:)))) % altitude measurement exists that year
                % disp(num2str(datestr(timenum(k))))
                 % intersection, crossdir=2 is descending; *1_new_out is
                 % the new array with intersection points
                [x_int z_int x1_new_out z1_new_out dum1_new_out dum2_new_out crossdir] = ...
                    findCrossings((cross(~isnan(altitude(k,:))))',altitude(k,~isnan(altitude(k,:))),...
                    [cross(1) cross(end)],[HibndMKL(k,m(i)+j-1) HibndMKL(k,m(i)+j-1)]);

                if any(crossdir == 2 | crossdir == 1)
                    dune_foot_upperMKL_xintersection = max(x_int(crossdir == 2 | crossdir == 1));
                else
                    dune_foot_upperMKL_xintersection = NaN;
                end
                
                dune_foot_upperMKL_cross(m(i)+j-1,k) = dune_foot_upperMKL_xintersection;
                
                if ~isnan(dune_foot_upperMKL_xintersection) % if it is non-NaN, calculate WGS84 coords.
                    dune_foot_upperMKL_latdum = interp1(cross',lat(m(i)+j-1,:),dune_foot_upperMKL_xintersection);
                    dune_foot_upperMKL_londum = interp1(cross',lon(m(i)+j-1,:),dune_foot_upperMKL_xintersection);
                    
                    dune_foot_upperMKL_lat(m(i)+j-1,k) = dune_foot_upperMKL_latdum; 
                    dune_foot_upperMKL_lon(m(i)+j-1,k) = dune_foot_upperMKL_londum;
                else
                    dune_foot_upperMKL_latdum = NaN; 
                    dune_foot_upperMKL_londum = NaN;
                    
                    dune_foot_upperMKL_lat(m(i)+j-1,k) = dune_foot_upperMKL_latdum; 
                    dune_foot_upperMKL_lon(m(i)+j-1,k) = dune_foot_upperMKL_londum;
                end
                
                clear dune_foot_xintersection
            end
        end
    end

end

% dune_foot in RD coordinates
[dune_foot_upperMKL_x, dune_foot_upperMKL_y] = convertCoordinates(dune_foot_upperMKL_lon,dune_foot_upperMKL_lat,'CS1.code',4326,'CS2.code',28992);

% remove old files
if exist('..\raw\dune_foot_upperMKL.mat', 'file')==2
  delete('..\raw\dune_foot_upperMKL.mat');
end

% save
save('..\raw\dune_foot_upperMKL.mat','dune_foot_upperMKL_cross',...
    'dune_foot_upperMKL_lat', 'dune_foot_upperMKL_lon', 'dune_foot_upperMKL_x', 'dune_foot_upperMKL_y');


%% dune_foot_threeNAP: intersection with +3m NAP
for i = 1:length(kv) % loop for area
%     num_tr  = m(i)-n(i)+1;
    num_tr  = n(i)-m(i)+1;

    for j = 1:num_tr % loop for transect
        altitude = squeeze(nc_varget([OPT.ncpath2 OPT.Filename2],'altitude',[0 m(i)+j-2 0],[-1 1 -1]));
        disp(['DF2: transect ' num2str(j) ' of area ' num2str(kv(i)) ...
            ' (' num2str(m(i)+j-1) '/' num2str(length(id)) ')']);
        
        for k = 1:length(time_MKL) % loop for year (...in MKL, because the upper bnd is required)
            if ~isempty(altitude(k,~isnan(altitude(k,:)))) % altitude measurement exists that year
                % disp(num2str(datestr(timenum(k))))
                 % intersection, crossdir=2 is descending; *1_new_out is
                 % the new array with intersection points
                [x_int z_int x1_new_out z1_new_out dum1_new_out dum2_new_out crossdir] = ...
                    findCrossings((cross(~isnan(altitude(k,:))))',altitude(k,~isnan(altitude(k,:))),...
                    [cross(1) cross(end)],[ThreeNAP(k,m(i)+j-1) ThreeNAP(k,m(i)+j-1)]);

                if any(crossdir == 2 | crossdir == 1)
                    dune_foot_threeNAP_xintersection = max(x_int(crossdir == 2 | crossdir == 1));
                else
                    dune_foot_threeNAP_xintersection = NaN;
                end
                
                dune_foot_threeNAP_cross(m(i)+j-1,k) = dune_foot_threeNAP_xintersection;
                
                if ~isnan(dune_foot_threeNAP_xintersection) % is it is non-NaN, calculate WGS84 coords.
                    dune_foot_threeNAP_latdum = interp1(cross',lat(m(i)+j-1,:),dune_foot_threeNAP_xintersection);
                    dune_foot_threeNAP_londum = interp1(cross',lon(m(i)+j-1,:),dune_foot_threeNAP_xintersection);
                    
                    dune_foot_threeNAP_lat(m(i)+j-1,k) = dune_foot_threeNAP_latdum; 
                    dune_foot_threeNAP_lon(m(i)+j-1,k) = dune_foot_threeNAP_londum;
                else
                    dune_foot_threeNAP_latdum = NaN; 
                    dune_foot_threeNAP_londum = NaN;
                    
                    dune_foot_threeNAP_lat(m(i)+j-1,k) = dune_foot_threeNAP_latdum; 
                    dune_foot_threeNAP_lon(m(i)+j-1,k) = dune_foot_threeNAP_londum;
                end
                
                clear dune_foot_threeNAP_xintersection
            end
        end
    end

end

% dune_foot in RD coordinates
[dune_foot_threeNAP_x, dune_foot_threeNAP_y] = convertCoordinates(dune_foot_threeNAP_lon,dune_foot_threeNAP_lat,'CS1.code',4326,'CS2.code',28992);

% remove old files
if exist('..\raw\dune_foot_threeNAP.mat', 'file')==2
  delete('..\raw\dune_foot_threeNAP.mat');
end

% save
save('..\raw\dune_foot_threeNAP.mat','dune_foot_threeNAP_cross',...
    'dune_foot_threeNAP_lat', 'dune_foot_threeNAP_lon', 'dune_foot_threeNAP_x', 'dune_foot_threeNAP_y');

% read from .mat file [with no use unless a division in shorter scripts is found convenient]
% load('dune_foot*.mat')

% remove old nc file
if exist(OPT.outputfile, 'file')==2
  delete(OPT.outputfile);
end

% fill the variables with new values (intersection from transect.nc profiles)
DF_upperMKL_tot(:,timeIDb(1):length(time_totaldatenum))       = dune_foot_upperMKL_cross;
DF_upperMKL_x_tot(:,timeIDb(1):length(time_totaldatenum))     = dune_foot_upperMKL_x;
DF_upperMKL_y_tot(:,timeIDb(1):length(time_totaldatenum))     = dune_foot_upperMKL_y;
DF_upperMKL_lat_tot(:,timeIDb(1):length(time_totaldatenum))   = dune_foot_upperMKL_lat;
DF_upperMKL_lon_tot(:,timeIDb(1):length(time_totaldatenum))   = dune_foot_upperMKL_lon;

DF_threeNAP_tot(:,timeIDb(1):length(time_totaldatenum))       = dune_foot_threeNAP_cross;
DF_threeNAP_x_tot(:,timeIDb(1):length(time_totaldatenum))     = dune_foot_threeNAP_x;
DF_threeNAP_y_tot(:,timeIDb(1):length(time_totaldatenum))     = dune_foot_threeNAP_y;
DF_threeNAP_lat_tot(:,timeIDb(1):length(time_totaldatenum))   = dune_foot_threeNAP_lat;
DF_threeNAP_lon_tot(:,timeIDb(1):length(time_totaldatenum))   = dune_foot_threeNAP_lon;

DF_threeNAP(:,timeIDb(1):length(time_totaldatenum))       = ThreeNAP';
DF_upperMKL(:,timeIDb(1):length(time_totaldatenum))       = HibndMKL';

%% assembly variables from strandlijnen.nc and .mat files all together
% epsg      = epsg
% id        = id
% areacode  = areacode
% areaname  = areaname
alongshore  = along;
cross_shore = cross;
time        = (time_totaldatenum - datenum('1970','yyyy'))';
time_bounds = [time_totaldatenum_beg - datenum('1970','yyyy'), time_totaldatenum_end - datenum('1970','yyyy')];

dune_foot_upperMKL       = DF_upperMKL'; % upper bnd for MKL calculation
dune_foot_upperMKL_cross = DF_upperMKL_tot';
dune_foot_upperMKL_x     = DF_upperMKL_x_tot';
dune_foot_upperMKL_y     = DF_upperMKL_y_tot';
dune_foot_upperMKL_lat   = DF_upperMKL_lat_tot';
dune_foot_upperMKL_lon   = DF_upperMKL_lon_tot';

dune_foot_threeNAP       = DF_threeNAP'; % +3m NAP
dune_foot_threeNAP_cross = DF_threeNAP_tot';
dune_foot_threeNAP_x     = DF_threeNAP_x_tot';
dune_foot_threeNAP_y     = DF_threeNAP_y_tot';
dune_foot_threeNAP_lat   = DF_threeNAP_lat_tot';
dune_foot_threeNAP_lon   = DF_threeNAP_lon_tot';

%% create empty file
nc_create_empty(OPT.outputfile)

%% Add overall meta info

attval =  sprintf('created on %s', datestr(now));
nc_attput(OPT.outputfile, nc_global, 'naming_authority'         , 'deltares.nl');
nc_attput(OPT.outputfile, nc_global, 'title'                    , 'dune foot position based on JarKus data and 1843-1998 database');
nc_attput(OPT.outputfile, nc_global, 'summary'                  , 'dune foot position along the Dutch coast since 1843, following +3m NAP definition (dune_foot_threeNAP) and the upper boundary of MKL definition (dune_foot_upperMKL)');
nc_attput(OPT.outputfile, nc_global, 'keywords'                 , 'Dune foot, Dutch coast');
nc_attput(OPT.outputfile, nc_global, 'standard_name_vocabulary' , 'http://cf-pcmdi.llnl.gov/documents/cf-standard-names/');
nc_attput(OPT.outputfile, nc_global, 'history'                  , sprintf('Data received from Rijkswaterstaat\nNetCDF created on %s by %s on computer %s\\%s with script %s',...
    datestr(now,31), getenv('USERNAME'), getenv('USERDOMAIN'), getenv('COMPUTERNAME'), 'convertDF2nc.m'));
nc_attput(OPT.outputfile, nc_global, 'institution'              , 'Rijkswaterstaat');
nc_attput(OPT.outputfile, nc_global, 'source'                   , 'https://svn.oss.deltares.nl/repos/openearthrawdata/trunk/rijkswaterstaat/ (Rijkswaterstaat)');
nc_attput(OPT.outputfile, nc_global, 'references'               , 'original source: http://www.rijkswaterstaat.nl , Deltares storage: http://opendap.deltares.nl/thredds/dodsC/opendap/rijkswaterstaat/DuneFoot/ , converted with script convertDF2nc.m');
nc_attput(OPT.outputfile, nc_global, 'creator_name'             , 'Rijkswaterstaat');
nc_attput(OPT.outputfile, nc_global, 'creator_url'              , 'www.rijkswaterstaat.nl');
nc_attput(OPT.outputfile, nc_global, 'creator_email'            , 'info@rijkswaterstaat.nl');
nc_attput(OPT.outputfile, nc_global, 'publisher_name'           , sprintf('%s (%s)', getenv('USERNAME'), getenv('USERID')));
nc_attput(OPT.outputfile, nc_global, 'publisher_url'            , 'www.deltares.nl');
nc_attput(OPT.outputfile, nc_global, 'publisher_email'          , 'Giorgio.Santinelli@deltares.nl');
% Other attributes
nc_attput(OPT.outputfile, nc_global, 'processing_level'         , 'preliminary');        
nc_attput(OPT.outputfile, nc_global, 'comment'                  , 'Most seaward point with descending JarKus transect.nc profile is taken as intersection for both dune_foot deifnitions, following +3m NAP definition (dune_foot_threeMKL) and the upper boundary of MKL definition (dune_foot_upperMKL)');
nc_attput(OPT.outputfile, nc_global, 'version'                  , '0.2');
nc_attput(OPT.outputfile, nc_global, 'conventions'              , 'CF-1.6');
nc_attput(OPT.outputfile, nc_global, 'license'                  ,[sprintf('These data can be used freely for research purposes provided that the following source is acknowledged: %s. ', 'RIJKSWATERSTAAT')...
            'disclaimer: This data is made available in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.']);

%% create dimensions
% add dimensions
nc_add_dimension(OPT.outputfile, 'time' , 0);
nc_add_dimension(OPT.outputfile,  'alongshore' , length(along));
nc_add_dimension(OPT.outputfile,  'cross_shore', length(cross));
nc_add_dimension(OPT.outputfile,  'bounds2'    , 2);
nc_add_dimension(OPT.outputfile,  'stringsize' , stringsize); % comments make CHARs' bigger in size

%% add variables
clear nc
i = 0;

%% add epsg

% smaller variables first
[lon2,lat2,OPT2] = convertCoordinates(0,0,'CS1.code',28992,'CS2.code',4326);
%[CoordinateSystems, Operations , CoordSysCart ,CoordSysGeo] = GetCoordinateSystems();
epsg        = 28992;

i = i + 1;
nc(i).Name      = 'epsg';
nc(i).Nctype    = 'int';
nc(i).Dimension = {};
nc(i).Attribute = struct('Name', ...
   {'grid_mapping_name', ...
    'semi_major_axis', ...
    'semi_minor_axis', ...
    'inverse_flattening', ...
    'latitude_of_projection_origin', ...
    'longitude_of_projection_origin', ...
    'false_easting', ...
    'false_northing', ...
    'scale_factor_at_projection_origin',...
    'comment'}, ...
    'Value', ...
    {OPT2.proj_conv1.method.name,    ...
    OPT2.CS1.ellips.semi_major_axis, ...
    OPT2.CS1.ellips.semi_minor_axis, ...
    OPT2.CS1.ellips.inv_flattening,  ...
    OPT2.proj_conv1.param.value(strcmp(OPT2.proj_conv1.param.name,'Latitude of natural origin'    )), ...
    OPT2.proj_conv1.param.value(strcmp(OPT2.proj_conv1.param.name,'Longitude of natural origin'   )),...
    OPT2.proj_conv1.param.value(strcmp(OPT2.proj_conv1.param.name,'False easting'                 )),...
    OPT2.proj_conv1.param.value(strcmp(OPT2.proj_conv1.param.name,'False northing'                )),...
    OPT2.proj_conv1.param.value(strcmp(OPT2.proj_conv1.param.name,'Scale factor at natural origin')),...
    'value is equal to EPSG code'});

%% same structure as in jarkus/profiles/transect.nc

% id
i = i + 1;
nc(i).Name         = 'id';
nc(i).Nctype       = nc_int;
nc(i).Dimension    = {'alongshore'};
nc(i).Attribute(1) = struct('Name', 'long_name' ,'Value', 'identifier');
nc(i).Attribute(2) = struct('Name', 'comment'   ,'Value', 'sum of area code (x1000000) and alongshore coordinate');

% areacode
i = i + 1;
nc(i).Name         = 'areacode';
nc(i).Nctype       = nc_double;
nc(i).Dimension    = {'alongshore'};
nc(i).Attribute(1) = struct('Name', 'long_name' ,'Value', 'area code');
nc(i).Attribute(2) = struct('Name', 'comment'   ,'Value', 'codes for the 15 coastal areas as defined by rijkswaterstaat');

% areaname
i = i + 1;
nc(i).Name         = 'areaname';
nc(i).Nctype       = nc_char;
nc(i).Dimension    = {'alongshore','stringsize'};
nc(i).Attribute(1) = struct('Name', 'long_name' ,'Value', 'area name');
nc(i).Attribute(2) = struct('Name', 'comment'   ,'Value', 'names for the 15 coastal areas as defined by rijkswaterstaat');

% alongshore
i = i + 1;
nc(i).Name         = 'alongshore';
nc(i).Nctype       = nc_double; 
nc(i).Dimension    = {'alongshore'};
nc(i).Attribute(1) = struct('Name', 'long_name' ,'Value', 'alongshore coordinate');
nc(i).Attribute(2) = struct('Name', 'units'     ,'Value', 'm');
nc(i).Attribute(3) = struct('Name', 'comment'   ,'Value', 'alongshore coordinate relative to the rsp (rijks strand paal)');

% cross_shore
i = i + 1;
nc(i).Name         = 'cross_shore';
nc(i).Nctype       = nc_double; 
nc(i).Dimension    = {'cross_shore'};
nc(i).Attribute(1) = struct('Name', 'long_name' ,'Value', 'cross-shore coordinate');
nc(i).Attribute(2) = struct('Name', 'units'     ,'Value', 'm');
nc(i).Attribute(3) = struct('Name', 'comment'   ,'Value', 'cross-shore coordinate relative to the rsp (rijks strand paal)');

% time
i = i + 1;
nc(i).Name         = 'time';
nc(i).Nctype       = nc_double;
nc(i).Dimension    = {'time'};
nc(i).Attribute(1) = struct('Name', 'standard_name' ,'Value', 'time');
nc(i).Attribute(2) = struct('Name', 'axis'          ,'Value', 'T');
nc(i).Attribute(3) = struct('Name', 'units'         ,'Value', 'days since 1970-01-01');
nc(i).Attribute(4) = struct('Name', 'cell_methods'  ,'Value', 'mean');
nc(i).Attribute(5) = struct('Name', 'bounds'        ,'Value', 'time_bounds');

% time bounds
i = i + 1;
nc(i).Name         = 'time_bounds';
nc(i).Nctype       = nc_double;
nc(i).Dimension    = {'time', 'bounds2'};
nc(i).Attribute    = struct('Name', 'standard_name', 'Value', 'time');
nc(i).Attribute(2) = struct('Name', 'units'        , 'Value', 'days since 1970-01-01');

% DF_upperMKL
i = i + 1;
nc(i).Name         = 'dune_foot_upperMKL';
nc(i).Nctype       = nc_double;
nc(i).Dimension    = {'time','alongshore'};
nc(i).Attribute(1) = struct('Name', 'long_name' ,'Value', 'upper boundary of MKL calculation');
nc(i).Attribute(2) = struct('Name', 'units'     ,'Value', 'm');
nc(i).Attribute(3) = struct('Name', 'comment'   ,'Value', 'dune foot height: upper boundary for calculation of the MKL');  

% DF_upperMKL_cross
i = i + 1;
nc(i).Name         = 'dune_foot_upperMKL_cross';
nc(i).Nctype       = nc_double;
nc(i).Dimension    = {'time','alongshore'};
nc(i).Attribute(1) = struct('Name', 'long_name' ,'Value', 'dune foot cross-shore intersection with upper boundary for calculation of the MKL');
nc(i).Attribute(2) = struct('Name', 'units'     ,'Value', 'm');
nc(i).Attribute(3) = struct('Name', 'comment'   ,'Value', 'dune foot intersection with transect profile altitude. Most seaward point with descending profile is taken');

% DF_upperMKL_x
i = i + 1;
nc(i).Name         = 'dune_foot_upperMKL_x';
nc(i).Nctype       = nc_double;
nc(i).Dimension    = {'time','alongshore'};
nc(i).Attribute(1) = struct('Name', 'long_name' ,'Value', 'dune foot cross-shore intersection in x-coordinate');
nc(i).Attribute(2) = struct('Name', 'units'     ,'Value', 'm');
nc(i).Attribute(3) = struct('Name', 'axis'      ,'Value', 'X');

% DF_upperMKL_y
i = i + 1;
nc(i).Name         = 'dune_foot_upperMKL_y';
nc(i).Nctype       = nc_double;
nc(i).Dimension    = {'time','alongshore'};
nc(i).Attribute(1) = struct('Name', 'long_name' ,'Value', 'dune foot cross-shore intersection in y-coordinate');
nc(i).Attribute(2) = struct('Name', 'units'     ,'Value', 'm');
nc(i).Attribute(3) = struct('Name', 'axis'      ,'Value', 'Y');

% DF_upperMKL_lat
i = i + 1;
nc(i).Name         = 'dune_foot_upperMKL_lat';
nc(i).Nctype       = nc_double;
nc(i).Dimension    = {'time','alongshore'};
nc(i).Attribute(1) = struct('Name', 'standard_name' ,'Value', 'latitude');
nc(i).Attribute(2) = struct('Name', 'units'         ,'Value', 'degree_north');
nc(i).Attribute(3) = struct('Name', 'axis'          ,'Value', 'X');

% DF_upperMKL_lon
i = i + 1;
nc(i).Name         = 'dune_foot_upperMKL_lon';
nc(i).Nctype       = nc_double;
nc(i).Dimension    = {'time','alongshore'};
nc(i).Attribute(1) = struct('Name', 'standard_name' ,'Value', 'longitude');
nc(i).Attribute(2) = struct('Name', 'units'         ,'Value', 'degree_east');
nc(i).Attribute(3) = struct('Name', 'axis'          ,'Value', 'Y');

% DF_threeNAP
i = i + 1;
nc(i).Name         = 'dune_foot_threeNAP';
nc(i).Nctype       = nc_double;
nc(i).Dimension    = {'time','alongshore'};
nc(i).Attribute(1) = struct('Name', 'long_name' ,'Value', '+3m NAP');
nc(i).Attribute(2) = struct('Name', 'units'     ,'Value', 'm');
nc(i).Attribute(3) = struct('Name', 'comment'   ,'Value', 'dune foot height: +3m NAP');  

% DF_threeNAP_cross
i = i + 1;
nc(i).Name         = 'dune_foot_threeNAP_cross';
nc(i).Nctype       = nc_double;
nc(i).Dimension    = {'time','alongshore'};
nc(i).Attribute(1) = struct('Name', 'long_name' ,'Value', 'dune foot cross-shore intersection with +3m NAP');
nc(i).Attribute(2) = struct('Name', 'units'     ,'Value', 'm');
nc(i).Attribute(3) = struct('Name', 'comment'   ,'Value', 'dune foot intersection with transect profile altitude. Most seaward point with descending profile is taken');

% DF_threeNAP_x
i = i + 1;
nc(i).Name         = 'dune_foot_threeNAP_x';
nc(i).Nctype       = nc_double;
nc(i).Dimension    = {'time','alongshore'};
nc(i).Attribute(1) = struct('Name', 'long_name' ,'Value', 'dune foot cross-shore intersection in x-coordinate');
nc(i).Attribute(2) = struct('Name', 'units'     ,'Value', 'm');
nc(i).Attribute(3) = struct('Name', 'axis'      ,'Value', 'X');

% DF_threeNAP_y
i = i + 1;
nc(i).Name         = 'dune_foot_threeNAP_y';
nc(i).Nctype       = nc_double;
nc(i).Dimension    = {'time','alongshore'};
nc(i).Attribute(1) = struct('Name', 'long_name' ,'Value', 'dune foot cross-shore intersection in y-coordinate');
nc(i).Attribute(2) = struct('Name', 'units'     ,'Value', 'm');
nc(i).Attribute(3) = struct('Name', 'axis'      ,'Value', 'Y');

% DF_threeNAP_lat
i = i + 1;
nc(i).Name         = 'dune_foot_threeNAP_lat';
nc(i).Nctype       = nc_double;
nc(i).Dimension    = {'time','alongshore'};
nc(i).Attribute(1) = struct('Name', 'standard_name' ,'Value', 'latitude');
nc(i).Attribute(2) = struct('Name', 'units'         ,'Value', 'degree_north');
nc(i).Attribute(3) = struct('Name', 'axis'          ,'Value', 'X');

% DF_threeNAP_lon
i = i + 1;
nc(i).Name         = 'dune_foot_threeNAP_lon';
nc(i).Nctype       = nc_double;
nc(i).Dimension    = {'time','alongshore'};
nc(i).Attribute(1) = struct('Name', 'standard_name' ,'Value', 'longitude');
nc(i).Attribute(2) = struct('Name', 'units'         ,'Value', 'degree_east');
nc(i).Attribute(3) = struct('Name', 'axis'          ,'Value', 'Y');

%% create variables with attributes
for i = 1:length(nc)
    if OPT.disp;disp(['adding ',num2str(i),' ',nc(i).Name]);end
    nc_addvar(OPT.outputfile, nc(i));
end

%% store variables
nc_varput(OPT.outputfile, 'id'            , id);
nc_varput(OPT.outputfile, 'areacode'      , areacode);
nc_varput(OPT.outputfile, 'areaname'      , areaname);
nc_varput(OPT.outputfile, 'alongshore'    , along);
nc_varput(OPT.outputfile, 'cross_shore'   , cross);
nc_varput(OPT.outputfile, 'time'          , time);
nc_varput(OPT.outputfile, 'time_bounds'   , time_bounds);

nc_varput(OPT.outputfile, 'dune_foot_upperMKL'       , dune_foot_upperMKL);
nc_varput(OPT.outputfile, 'dune_foot_upperMKL_cross' , dune_foot_upperMKL_cross);
nc_varput(OPT.outputfile, 'dune_foot_upperMKL_x'     , dune_foot_upperMKL_x);
nc_varput(OPT.outputfile, 'dune_foot_upperMKL_y'     , dune_foot_upperMKL_y);
nc_varput(OPT.outputfile, 'dune_foot_upperMKL_lat'   , dune_foot_upperMKL_lat);
nc_varput(OPT.outputfile, 'dune_foot_upperMKL_lon'   , dune_foot_upperMKL_lon);

nc_varput(OPT.outputfile, 'dune_foot_threeNAP'       , dune_foot_threeNAP);
nc_varput(OPT.outputfile, 'dune_foot_threeNAP_cross' , dune_foot_threeNAP_cross);
nc_varput(OPT.outputfile, 'dune_foot_threeNAP_x'     , dune_foot_threeNAP_x);
nc_varput(OPT.outputfile, 'dune_foot_threeNAP_y'     , dune_foot_threeNAP_y);
nc_varput(OPT.outputfile, 'dune_foot_threeNAP_lat'   , dune_foot_threeNAP_lat);
nc_varput(OPT.outputfile, 'dune_foot_threeNAP_lon'   , dune_foot_threeNAP_lon);

%% check
if OPT.dump
    nc_dump(OPT.outputfile);
end

%% end
