%% settings
% general settings
OPT.outputfile   = 'DF_2nd_deriv.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      = 'http://opendap.deltares.nl/thredds/dodsC/opendap/rijkswaterstaat/strandlijnen/';
OPT.Filename1    = 'strandlijnen.nc';
% % JarKus database - dataset 1965 - ...
OPT.ncpath2      = 'C:\Users\fuentesm\Marilu\Deltares\Projects\Coastviewer\jarkus\scripts\';
OPT.Filename2    = 'transect_r20240405.nc';

% OPT.ncpath2      = 'http://opendap.deltares.nl/thredds/dodsC/opendap/rijkswaterstaat/jarkus/profiles/';
% OPT.Filename2    = 'transect.nc';

% % MKL database for the upper horizontal limit
% OPT.ncpath3      =
% 'http://opendap.deltares.nl/thredds/dodsC/opendap/rijkswaterstaat/BKL_TKL_MKL/';%
% % local
OPT.ncpath3      = 'C:\Users\fuentesm\Marilu\Deltares\Projects\Coastviewer\MHW_MLW\processed\';
OPT.Filename3    = 'MHW_MLW_r20240414.nc';

%% Read data 
% 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;

id_MKL     = nc_varget([OPT.ncpath3 OPT.Filename3],'id');
[idID_STa, idID_STb] = ismember(id,id_MKL);

df= load('..\Scripts for 2nd derivative method\DF_2nd_deriv.mat');

% Preallocation and definition
DF_2nd_deriv_tot     = NaN(length(id_MKL),length(timej));
DF_2nd_deriv_x_tot   = NaN(length(id_MKL),length(timej));
DF_2nd_deriv_y_tot   = NaN(length(id_MKL),length(timej));
DF_2nd_deriv_lat_tot = NaN(length(id_MKL),length(timej));
DF_2nd_deriv_lon_tot = NaN(length(id_MKL),length(timej));
DF_2nd_deriv         = NaN(length(id_MKL),length(timej));

%% 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 definition
dune_foot_2nd_deriv_cross = NaN(length(timej),length(areacode));
dune_foot_2nd_deriv_lat   = NaN(length(timej),length(areacode));
dune_foot_2nd_deriv_lon   = NaN(length(timej),length(areacode));

disp(['check if ' num2str(length(df.id)) ' is equal to ' num2str(length(id)) '. If not, get worried.'])
for t=1:size(df.xDF,1)
    for trans=1:size(df.xDF,2)
        [c, index]=min(abs(cross-df.xDF(t,trans)));
        
        if isnan(c)
            dune_foot_2nd_deriv_lat(t,trans)    = NaN;
            dune_foot_2nd_deriv_lon(t,trans)    = NaN;
        
        else
            dune_foot_2nd_deriv_lat(t,trans)    = lat(trans,index);
            dune_foot_2nd_deriv_lon(t,trans)    = lon(trans,index);
            dune_foot_2nd_deriv_cross(t,trans) = cross(index);
        end
    end
end
dune_foot_2nd_deriv = df.zDF;
% dune_foot in RD coordinates
[dune_foot_2nd_deriv_x, dune_foot_2nd_deriv_y] = convertCoordinates(dune_foot_2nd_deriv_lon,dune_foot_2nd_deriv_lat,'CS1.code',4326,'CS2.code',28992);

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

% save
save('dune_foot_2nd_deriv.mat','dune_foot_2nd_deriv_cross',...
    'dune_foot_2nd_deriv_lat', 'dune_foot_2nd_deriv_lon', 'dune_foot_2nd_deriv_x', 'dune_foot_2nd_deriv_y');

% remove old files in raw folder
if exist('../../../raw/dune_foot_2nd_deriv.mat', 'file')==2
  delete('../../../raw/dune_foot_2nd_deriv.mat');
end

% save in raw folder too
save('../../../raw/dune_foot_2nd_deriv.mat','dune_foot_2nd_deriv_cross',...
    'dune_foot_2nd_deriv_lat', 'dune_foot_2nd_deriv_lon', 'dune_foot_2nd_deriv_x', 'dune_foot_2nd_deriv_y');


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

% remove old also from processed foledr
if exist(['..', filesep, '..', filesep, '..', filesep, 'processed', filesep, OPT.outputfile], 'file')==2
  delete(['..', filesep, '..', filesep, '..', filesep, 'processed', filesep, OPT.outputfile]);
end

%% make netcdf [quick]
% 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');
nc_attput(OPT.outputfile, nc_global, 'summary'                  , 'dune foot position along the Dutch coast since 1843, following the second derivative method');
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'          , 'diamantidou_el@hotmail.com');
% Other attributes
nc_attput(OPT.outputfile, nc_global, 'processing_level'         , 'preliminary');        
nc_attput(OPT.outputfile, nc_global, 'comment'                  , 'Dune foot position based on the second derivative method');
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_adddim(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_2nd_deriv
i = i + 1;
nc(i).Name         = 'dune_foot_2nd_deriv';
nc(i).Nctype       = nc_double;
nc(i).Dimension    = {'time','alongshore'};
nc(i).Attribute(1) = struct('Name', 'long_name' ,'Value', 'dune foot based on second derivative method');
nc(i).Attribute(2) = struct('Name', 'units'     ,'Value', 'm');
nc(i).Attribute(3) = struct('Name', 'comment'   ,'Value', 'dune foot height: second derivative menthod');  

% DF_2nd_deriv_cross
i = i + 1;
nc(i).Name         = 'dune_foot_2nd_deriv_cross';
nc(i).Nctype       = nc_double;
nc(i).Dimension    = {'time','alongshore'};
nc(i).Attribute(1) = struct('Name', 'long_name' ,'Value', 'dune foot based on the second derivative method');
nc(i).Attribute(2) = struct('Name', 'units'     ,'Value', 'm');
nc(i).Attribute(3) = struct('Name', 'comment'   ,'Value', 'dune foot based on the second derivative method');

% DF_2nd_deriv_x
i = i + 1;
nc(i).Name         = 'dune_foot_2nd_deriv_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_2nd_deriv_y
i = i + 1;
nc(i).Name         = 'dune_foot_2nd_deriv_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_2nd_deriv_lat
i = i + 1;
nc(i).Name         = 'dune_foot_2nd_deriv_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_2nd_deriv_lon
i = i + 1;
nc(i).Name         = 'dune_foot_2nd_deriv_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'          , timej);
nc_varput(OPT.outputfile, 'time_bounds'   , timebj);

nc_varput(OPT.outputfile, 'dune_foot_2nd_deriv'       , dune_foot_2nd_deriv);
nc_varput(OPT.outputfile, 'dune_foot_2nd_deriv_cross' , dune_foot_2nd_deriv_cross);
nc_varput(OPT.outputfile, 'dune_foot_2nd_deriv_x'     , dune_foot_2nd_deriv_x);
nc_varput(OPT.outputfile, 'dune_foot_2nd_deriv_y'     , dune_foot_2nd_deriv_y);
nc_varput(OPT.outputfile, 'dune_foot_2nd_deriv_lat'   , dune_foot_2nd_deriv_lat);
nc_varput(OPT.outputfile, 'dune_foot_2nd_deriv_lon'   , dune_foot_2nd_deriv_lon);

%% copy in processed folder, too
copyfile(OPT.outputfile, ['..', filesep, '..', filesep, '..', filesep, 'processed', filesep, OPT.outputfile])

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

%% end

