% BKL_TKL_TND2nc
%
% In order to be used in MorphAn, the file BKL_TKL_TND.nc in 
% https://repos.deltares.nl/repos/OpenEarthRawData/trunk/rijkswaterstaat/coastal_indicators_morphAn\processed
% has been changed.
% 
% Generates BKL_TKL_TND.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@deltares.nl
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%% settings
% BKL, TKL, TND
OPT.ncpath1      = 'http://opendap.deltares.nl/thredds/dodsC/opendap/rijkswaterstaat/BKL_TKL_MKL/';
OPT.Filename1    = 'BKL_TKL_TND.nc';

% jarkus database for transect info
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\BKL_TKL_TND.nc';
OPT.disp         = true;
OPT.dump         = true;

%% BKL, TKL, TND from Deltares repository
% Read data form jarkus transect.nc
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');

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

time     = nc_varget([OPT.ncpath1 OPT.Filename1],'time'); % fairly one year more!!

basal_coastline    =  (nc_varget([OPT.ncpath1,OPT.Filename1],'basal_coastline'))';
testing_coastline  =  (nc_varget([OPT.ncpath1,OPT.Filename1],'testing_coastline'))';
trend              =  (nc_varget([OPT.ncpath1,OPT.Filename1],'trend'))';

%% 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'         , 'basal coastline, testing coastline, trend, based on Jarkus data');
nc_attput(OPT.outputfile, nc_global, 'institution'   , 'Rijkswaterstaat'); % Deltares?
nc_attput(OPT.outputfile, nc_global, 'source'        , 'measurements and extrapolation');
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/ , converted with script BKL_TKL_TND2nc.m');
nc_attput(OPT.outputfile, nc_global, 'email'         , 'giorgio.santinelli@deltares.nl');
nc_attput(OPT.outputfile, nc_global, 'comment'       , '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.');
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: Rijkswaterstaat');
% 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 parameter per transect
nc_adddim(OPT.outputfile, 'stringsize', 100);

%% add variables
clear nc
i = 0;

%% add epsg

% smaller variables first
[lon_dummy, lat_dummy, 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');

% indicator_value (basal_coastline)
i = i + 1;
nc(i).Name         = 'basal_coastline';
nc(i).Nctype       = 'double';
nc(i).Dimension    = {'alongshore','time'};
nc(i).Attribute(1) = struct('Name', 'indicator_name','Value', 'basal_coastline');
nc(i).Attribute(2) = struct('Name', 'long_name'     ,'Value', 'BKL position');
nc(i).Attribute(3) = struct('Name', 'units'         ,'Value', 'm');
nc(i).Attribute(4) = struct('Name', 'comment'       ,'Value', 'BKL position defined by rijkswaterstaat');
nc(i).Attribute(5) = struct('Name', '_FillValue'    ,'Value', NaN);

% indicator_value (testing_coastline)
i = i + 1;
nc(i).Name         = 'testing_coastline';
nc(i).Nctype       = 'double';
nc(i).Dimension    = {'alongshore','time'};
nc(i).Attribute(1) = struct('Name', 'indicator_name','Value', 'testing_coastline');
nc(i).Attribute(2) = struct('Name', 'long_name'     ,'Value', 'TKL position');
nc(i).Attribute(3) = struct('Name', 'units'         ,'Value', 'm');
nc(i).Attribute(4) = struct('Name', 'comment'       ,'Value', 'TKL position defined by rijkswaterstaat');
nc(i).Attribute(5) = struct('Name', '_FillValue'    ,'Value', NaN);

% indicator_value (trend)
i = i + 1;
nc(i).Name         = 'trend';
nc(i).Nctype       = 'double';
nc(i).Dimension    = {'alongshore','time'};
nc(i).Attribute(1) = struct('Name', 'indicator_name','Value', 'trend');
nc(i).Attribute(2) = struct('Name', 'long_name'     ,'Value', 'erosion trend');
nc(i).Attribute(3) = struct('Name', 'units'         ,'Value', 'm / year');
nc(i).Attribute(4) = struct('Name', 'comment'       ,'Value', 'erosion trend defined by rijkswaterstaat');
nc(i).Attribute(5) = struct('Name', '_FillValue'    ,'Value', NaN);

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

% with respect to the Rsp
nc_varput(OPT.outputfile, 'basal_coastline'         , basal_coastline);
nc_varput(OPT.outputfile, 'testing_coastline'       , testing_coastline);
nc_varput(OPT.outputfile, 'trend'                   , trend);

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

%% end

