% convertDF2nc
%
% Generates DF.nc (for MorphAn) which contains
%
% Dune Foot  defined as the crossing between the water level at +3.00m nap
%            and the coastline profile
% 
% Data are taken from Jarkus dataset (1965 - 2010) @
% http://opendap.deltares.nl/thredds/dodsC/opendap/rijkswaterstaat/jarkus/profiles/transect.nc
%
% NB. strandlijnen dataset (1843 - 1998) is already in a .nc format @
% https://repos.deltares.nl/repos/OpenEarthRawData/trunk/rijkswaterstaat/strandlijnen/processed/strandlijnen.nc
% 
% STRANDLIJNEN dataset is NOT MERGED WITH JARKUS YET! 
% - Decision on how combining them is required
%
% giorgio.santinelli@deltares.nl
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% settings
% strandlijnen database - dataset 1843 - 1998
OPT.ncpath1      = 'https://repos.deltares.nl/repos/OpenEarthRawData/trunk/rijkswaterstaat/strandlijnen/processed/';
OPT.ncpath1      = 'F:\checkout\OpenEarthRawData\trunk\rijkswaterstaat\strandlijnen\processed\';
OPT.Filename1    = 'strandlijnen.nc';
% jarkus database - dataset 1965 - 2010
OPT.ncpath2      = 'http://opendap.deltares.nl/thredds/dodsC/opendap/rijkswaterstaat/jarkus/profiles/';
OPT.ncpath2      = 'F:\JARKUSdata\';
OPT.Filename2    = 'transect.nc';
% general settings
OPT.outputfile   = '..\processed\DF.nc';
OPT.disp         = true;
OPT.dump         = true;

%% jarkus database - dataset 1965 - 2010
%% Read data
id       = nc_varget([OPT.ncpath2 OPT.Filename2],'id');
MHW      = nc_varget([OPT.ncpath2 OPT.Filename2],'mean_high_water');
MLW      = nc_varget([OPT.ncpath2 OPT.Filename2],'mean_low_water');
lat      = nc_varget([OPT.ncpath2 OPT.Filename2],'lat');
lon      = nc_varget([OPT.ncpath2 OPT.Filename2],'lon');
cross    = nc_varget([OPT.ncpath2 OPT.Filename2],'cross_shore');
along    = nc_varget([OPT.ncpath2 OPT.Filename2],'alongshore');
areaname = nc_varget([OPT.ncpath2 OPT.Filename2],'areaname');
areacode = nc_varget([OPT.ncpath2 OPT.Filename2],'areacode');
time     = nc_varget([OPT.ncpath2 OPT.Filename2],'time');

rsp_x    = nc_varget([OPT.ncpath2 OPT.Filename2],'rsp_x');
rsp_y    = nc_varget([OPT.ncpath2 OPT.Filename2],'rsp_y');
x        = nc_varget([OPT.ncpath2 OPT.Filename2],'x');
y        = nc_varget([OPT.ncpath2 OPT.Filename2],'y');
angle    = nc_varget([OPT.ncpath2 OPT.Filename2],'angle');

dune_foot_name = 'dune_foot';

%% Evaluation of DF

% counters
[kv, m, ~] = unique(areacode);
n          = vertcat(1, m(1:end-1)+1);
areakml    = cell(1,length(kv));

% definition
DF_lat  = NaN(length(areacode),length(time)); DF_lon  = NaN(length(areacode),length(time));
DF_rsp  = NaN(length(areacode),length(time));

for i = 1:length(kv) % instead of test 2:2 Ameland
    num_tr  = m(i)-n(i)+1;
i
    for j = 1:num_tr
        altitude = nc_varget([OPT.ncpath2 OPT.Filename2],'altitude',[0 n(i)+j-2 0],[-1 1 -1]);
        
        for k = 1:length(time)
            if ~isempty(altitude(k,~isnan(altitude(k,:))))
                
            dumDF   = max(findCrossings((cross(~isnan(altitude(k,:))))',altitude(k,~isnan(altitude(k,:))),...
                [cross(1) cross(end)],[3.00 3.00]));
            if ~isempty(dumDF)
                DFlatdum = interp1(cross',lat(n(i)+j-1,:),dumDF);
                DFlondum = interp1(cross',lon(n(i)+j-1,:),dumDF);
                DF_lat(n(i)+j-1,k) = DFlatdum; DF_lon(n(i)+j-1,k) = DFlondum;
            else
                DFlatdum = NaN; DFlondum = NaN;
                DF_lat(n(i)+j-1,k) = DFlatdum; DF_lon(n(i)+j-1,k) = DFlondum;
            end
            
            end
        end
    end

%     DF per each kustvak
    DF_lat_kv  = DF_lat(n(i):m(i),:);   DF_lon_kv  = DF_lon(n(i):m(i),:);

end

% DF in RD coordinates
[DF_x, DF_y] = convertCoordinates(DF_lon,DF_lat,'CS1.code',4326,'CS2.code',28992);

% DF in cross-shore coordinates wrt RSP
for j = 1:length(id) 
    for k = 1:length(time)
        if (angle(j) >= 180 && angle(j) <= 360) % || (angle >= 90 && angle <= 180)
            if (DF_x(j,k) - rsp_x(j)) >= 0
                DF_rsp(j,k) = - sqrt((DF_x(j,k) - rsp_x(j))^2 + (DF_y(j,k) - rsp_y(j))^2);
            else
                DF_rsp(j,k) = sqrt((DF_x(j,k) - rsp_x(j))^2 + (DF_y(j,k) - rsp_y(j))^2);
            end
        else
            if (DF_x(j,k) - rsp_x(j)) >= 0
                DF_rsp(j,k) = sqrt((DF_x(j,k) - rsp_x(j))^2 + (DF_y(j,k) - rsp_y(j))^2);
            else
                DF_rsp(j,k) = - sqrt((DF_x(j,k) - rsp_x(j))^2 + (DF_y(j,k) - rsp_y(j))^2);
            end
        end
    end
end

% for j = 1:length(id)
%     for k = 1:length(time)
%         DF_rsp(j,k) = interp1(x(j,:), cross', DF_x(j,k));
%     end
% end


%% create empty file
nc_create_empty(OPT.outputfile)

%% Add overall meta info
%  http://cf-pcmdi.llnl.gov/documents/cf-conventions/1.5/cf-conventions.html#description-of-file-contents
%------------------

attval =  sprintf('created on %s', datestr(now));

nc_attput(OPT.outputfile, nc_global, 'title'         , 'dune foot position based on Jarkus data');
nc_attput(OPT.outputfile, nc_global, 'institution'   , 'Rijkswaterstaat'); % Deltares?
nc_attput(OPT.outputfile, nc_global, 'source'        , 'on shore measurements');
nc_attput(OPT.outputfile, nc_global, 'history'       , ['Data processed from Jarkus database, ', attval]);
nc_attput(OPT.outputfile, nc_global, 'references'    , 'original source: http://www.rijkswaterstaat.nl , Deltares storage: http://opendap.deltares.nl/thredds/dodsC/opendap/rijkswaterstaat/jarkus/ , converted with script convertDF2nc.m');
% nc_attput(OPT.outputfile, nc_global, 'email'         , '');
% nc_attput(OPT.outputfile, nc_global, 'comment'       , '');
% nc_attput(OPT.outputfile, nc_global, 'version'       , '0.1');
nc_attput(OPT.outputfile, nc_global, 'Conventions'   , 'CF-1.5');
% spatial reference using wkt. approach used by gdal.
nc_attput(OPT.outputfile, nc_global, 'spatial_ref', 'COMPD_CS["Amersfoort / RD New + NAP",PROJCS["Amersfoort / RD New",GEOGCS["Amersfoort",DATUM["Amersfoort",SPHEROID["Bessel 1841",6377397.155,299.1528128,AUTHORITY["EPSG","7004"]],TOWGS84[565.04,49.91,465.84,-0.40939438743923684,-0.35970519561431136,1.868491000350572,0.8409828680306614],AUTHORITY["EPSG","6289"]],PRIMEM["Greenwich",0.0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.017453292519943295],AXIS["Geodetic latitude",NORTH],AXIS["Geodetic longitude",EAST],AUTHORITY["EPSG","4289"]],PROJECTION["Oblique Stereographic",AUTHORITY["EPSG","9809"]],PARAMETER["central_meridian",5.387638888888891],PARAMETER["latitude_of_origin",52.15616055555556],PARAMETER["scale_factor",0.9999079],PARAMETER["false_easting",155000.0],PARAMETER["false_northing",463000.0],UNIT["m",1.0],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","28992"]],VERT_CS["Normaal Amsterdams Peil",VERT_DATUM["Normaal Amsterdams Peil",2005,AUTHORITY["EPSG","5109"]],UNIT["m",1.0],AXIS["Gravity-related height",UP],AUTHORITY["EPSG","5709"]],AUTHORITY["EPSG","7415"]]');
% nc_attput(OPT.outputfile, nc_global, 'CF:featureType', 'stationTimeSeries');
% nc_attput(OPT.outputfile, nc_global, 'terms_for_use' , 'These data can be used freely for research purposes provided that the following source is acknowledged: ARGUS??');
% nc_attput(OPT.outputfile, nc_global, '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', length(time)); % 0);
nc_adddim(OPT.outputfile, 'alongshore', length(along)); % one DF per transect
nc_adddim(OPT.outputfile, 'stringsize', 100);

%% 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'});

% id
i = i + 1;
nc(i).Name         = 'id';
nc(i).Nctype       = '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       = 'int';
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       = '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       = '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 within the 15 coastal areas as defined by rijkswaterstaat');

% angle
i = i + 1;
nc(i).Name         = 'angle';
nc(i).Nctype       = 'double';
nc(i).Dimension    = {'alongshore'};
nc(i).Attribute(1) = struct('Name', 'long_name' ,'Value', 'angle of transect');
nc(i).Attribute(2) = struct('Name', 'units'     ,'Value', 'degrees');
nc(i).Attribute(3) = struct('Name', 'comment'   ,'Value', 'positive clockwise 0 north');

% time
i = i + 1;
nc(i).Name         = 'time';
nc(i).Nctype       = '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 00:00 +1:00');
nc(i).Attribute(4) = struct('Name', 'comment'       ,'Value', 'measurement year');

% DF_lat
i = i + 1;
nc(i).Name         = 'lat';
nc(i).Nctype       = 'double';
nc(i).Dimension    = {'alongshore','time'};
nc(i).Attribute(1) = struct('Name', 'long_name' ,'Value', 'dune foot position');
nc(i).Attribute(2) = struct('Name', 'units'     ,'Value', 'degrees_north');
nc(i).Attribute(3) = struct('Name', 'comment'   ,'Value', 'dune foot is defined by the intersection of the profile with the line +3.00m nap');  

% DF_lon
i = i + 1;
nc(i).Name         = 'lon';
nc(i).Nctype       = 'double';
nc(i).Dimension    = {'alongshore','time'};
nc(i).Attribute(1) = struct('Name', 'long_name' ,'Value', 'dune foot position');
nc(i).Attribute(2) = struct('Name', 'units'     ,'Value', 'degrees_east');
nc(i).Attribute(3) = struct('Name', 'comment'   ,'Value', 'dune foot is defined by the intersection of the profile with the line +3.00m nap');  

% DF_x
i = i + 1;
nc(i).Name         = 'x';
nc(i).Nctype       = 'double';
nc(i).Dimension    = {'alongshore','time'};
nc(i).Attribute(1) = struct('Name', 'long_name' ,'Value', 'dune foot position');
nc(i).Attribute(2) = struct('Name', 'units'     ,'Value', 'm');
nc(i).Attribute(3) = struct('Name', 'axis'      ,'Value', 'X');
nc(i).Attribute(4) = struct('Name', 'comment'   ,'Value', 'dune foot is defined by the intersection of the profile with the line +3.00m nap');

% DF_y
i = i + 1;
nc(i).Name         = 'y';
nc(i).Nctype       = 'double';
nc(i).Dimension    = {'alongshore','time'};
nc(i).Attribute(1) = struct('Name', 'long_name' ,'Value', 'dune foot position');
nc(i).Attribute(2) = struct('Name', 'units'     ,'Value', 'm');
nc(i).Attribute(3) = struct('Name', 'axis'      ,'Value', 'Y');
nc(i).Attribute(4) = struct('Name', 'comment'   ,'Value', 'dune foot is defined by the intersection of the profile with the line +3.00m nap');

% indicator_name (DF_rsp)
i = i + 1;
nc(i).Name         = 'indicator_name';
nc(i).Nctype       = 'char';
nc(i).Dimension    = {'stringsize'};
nc(i).Attribute(1) = struct('Name', 'long_name' ,'Value', 'name of the indicator (dune foot position wrt rijkstrandpaal)');

% indicator_value (DF_rsp)
i = i + 1;
nc(i).Name         = 'indicator_value';
nc(i).Nctype       = 'double';
nc(i).Dimension    = {'alongshore','time'};
nc(i).Attribute(1) = struct('Name', 'long_name' ,'Value', 'dune foot position wrt rijkstrandpaal');
nc(i).Attribute(2) = struct('Name', 'units'     ,'Value', 'm');
nc(i).Attribute(3) = struct('Name', 'comment'   ,'Value', 'dune foot is defined by the intersection of the profile with the line +3.00m nap');

%% 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, 'angle'           , angle);
nc_varput(OPT.outputfile, 'time'            , time);
nc_varput(OPT.outputfile, 'lat'             , DF_lat);
nc_varput(OPT.outputfile, 'lon'             , DF_lon);
nc_varput(OPT.outputfile, 'x'               , DF_x);
nc_varput(OPT.outputfile, 'y'               , DF_y);
nc_varput(OPT.outputfile, 'indicator_name'  , dune_foot_name);
nc_varput(OPT.outputfile, 'indicator_value' , DF_rsp);

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

%% end
