% netcdfStrandbreedte.m
%
% Generates strandbreedte.nc which contains the following coastal indicator: 
% strandbreedteGHW   (de strandbreedte bij Gemiddeld Hoog Water)
% strandbreedteGLW   (de strandbreedte bij Gemiddeld Laag Water)
%                              
%
% strandbreedteGHW and strandbreedteGLW are obtained by HKV.
%     The information are available from Deltares repository @
%     https://repos.deltares.nl/repos/OpenEarthRawData/trunk/rijkswaterstaat/kustlijnkaartenboek/raw/mkl/
%
% giorgio.santinelli@deltares.nl
%
% adapted by Wim van Balen and Vincent Vuik
% vuik@hkv.nl
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%% settings

% jarkus database for transect info
OPT.ncpath3      = 'http://opendap.deltares.nl/thredds/dodsC/opendap/rijkswaterstaat/jarkus/profiles/';
OPT.Filename3    = 'transect.nc';

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

%% Read data form jarkus transect.nc
id         = nc_varget([OPT.ncpath3 OPT.Filename3],'id');
MHW        = nc_varget([OPT.ncpath3 OPT.Filename3],'mean_high_water');
MLW        = nc_varget([OPT.ncpath3 OPT.Filename3],'mean_low_water');
lat        = nc_varget([OPT.ncpath3 OPT.Filename3],'lat');
lon        = nc_varget([OPT.ncpath3 OPT.Filename3],'lon');
cross      = nc_varget([OPT.ncpath3 OPT.Filename3],'cross_shore');
along      = nc_varget([OPT.ncpath3 OPT.Filename3],'alongshore');
areaname   = nc_varget([OPT.ncpath3 OPT.Filename3],'areaname');
areacode   = nc_varget([OPT.ncpath3 OPT.Filename3],'areacode');
time       = nc_varget([OPT.ncpath3 OPT.Filename3],'time');

jarkusyear = str2num(datestr(nc_cf_time([OPT.ncpath3 OPT.Filename3]),'yyyy'));

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

%% create the matrix for each parameter, with a structure similar to the jarkus file

% define the strandbreedte
bestand                = ['strandbreedtes_droog.xls'];
raaien                 = xlsread(bestand,'totaal','A2:A1424');
jaren                  = xlsread(bestand,'totaal','B1:AU1');
output                 = xlsread(bestand,'totaal','B2:AU1424');
for i=1:length(raaien);
    lokid(i)           = find(id==raaien(i)); 
end
strandbreedte          = NaN.*ones(46,2178);
strandbreedte(:,lokid) = output';
strandbreedteGHW       = strandbreedte;
clear('output','lokid');

% define the strandbreedte
bestand                = ['strandbreedtes_totaal.xls'];
raaien                 = xlsread(bestand,'totaal','A2:A1418');
jaren                  = xlsread(bestand,'totaal','B1:AU1');
output                 = xlsread(bestand,'totaal','B2:AU1418');
for i=1:length(raaien);
    lokid(i)           = find(id==raaien(i)); 
end
strandbreedte          = NaN.*ones(46,2178);
strandbreedte(:,lokid) = output';
strandbreedteGLW       = strandbreedte;
clear('output','lokid');



%% 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'         , 'beach width at MHW and MLW');
nc_attput(OPT.outputfile, nc_global, 'institution'   , 'HKV Lijn in Water');
nc_attput(OPT.outputfile, nc_global, 'email'         , 'balen@hkv.nl');
nc_attput(OPT.outputfile, nc_global, 'version'       , '1.0');

%% create dimensions
% add dimensions
nc_adddim(OPT.outputfile, 'time'      , length(time) ); 
nc_adddim(OPT.outputfile, 'alongshore', length(along)); 
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);
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');

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

% strandbreedteGHW
i = i + 1;
nc(i).Name         = 'beach_width_at_MHW';
nc(i).Nctype       = 'double';
nc(i).Dimension    = {'time','alongshore'};
nc(i).Attribute(1) = struct('Name', 'long_name' ,'Value', 'Beach width at MHW');
nc(i).Attribute(2) = struct('Name', 'units'     ,'Value', 'meters');
nc(i).Attribute(3) = struct('Name', 'comment'   ,'Value', 'Distance between dune boundary and sea boundary at MHW');

% strandbreedteGLW
i = i + 1;
nc(i).Name         = 'beach_width_at_MLW';
nc(i).Nctype       = 'double';
nc(i).Dimension    = {'time','alongshore'};
nc(i).Attribute(1) = struct('Name', 'long_name' ,'Value', 'Beach width at MLW');
nc(i).Attribute(2) = struct('Name', 'units'     ,'Value', 'meters');
nc(i).Attribute(3) = struct('Name', 'comment'   ,'Value', 'Distance between dune boundary and sea boundary at MLW');  

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

nc_varput(OPT.outputfile, 'beach_width_at_MHW'      , strandbreedteGHW);
nc_varput(OPT.outputfile, 'beach_width_at_MLW'      , strandbreedteGLW);

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

%% end