% MW2nc
%
% during years from 1843 to 1964 strandlijnen.nc is used. For next years,
% the intersection of MHW and MLW with Jarkus profile transect.nc is used.
% MHW and MLW height values are considered as static for the dataset 1965-now.
% These values are stored in the Jarkus dataset, calculated in 2018. For more info:
% https://svn.oss.deltares.nl/repos/openearthrawdata/trunk/rijkswaterstaat/jarkus/raw/info/1230043-000-ZKS-0006-m-Specialistisch%20advies.pdf
% 
% definition of Jarkus intersection and test are respectively in
% MW_Jarkus.nc and compareMW.nc
%
% Generates Mean Low Water and Mean High Water position in NetCDF file format in
% http://opendap.deltares.nl/thredds/dodsC/opendap/rijkswaterstaat/MHW_MLW/MHW_MLW.nc
%
% MHW_MLW.nc contains:
% mean_high_water position defined as the mean high water intersection 
% with JarKus profile altitude. Most seaward point with descending profile is taken.
% 
% mean_low_water position defined as the mean low water intersection 
% with JarKus profile altitude. Most seaward point with descending profile is taken.
% 
% NB. 
% a) Mean High Water and Mean Low Water positions are taken from STRANDLIJNEN dataset (1843 - 1998), already in .nc format @
% https://repos.deltares.nl/repos/OpenEarthRawData/trunk/rijkswaterstaat/strandlijnen/processed/strandlijnen.nc
% b) Mean High Water and Mean Low Water positions are 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) Mean High Water and Mean Low Water in 1843 - 1964 are kept as they are in STRANDLIJNEN dataset
% e) Definition of intersection with JarKus profile and script for tests are respectively in MW_Jarkus.mat and compareMW.mat
%
%
% $Id$
% $Date$
% $Author$
% $Revision$
% $HeadURL$
% $Keywords$ 
%
%%

%% read files
% general settings
OPT.fillvalue    = NaN;
OPT.refdatenum   = datenum(1970,1,1);
OPT.outputfile   = '..\processed\MHW_MLW.nc';

% jarkus database for info about RSP, latlon, RDNew
OPT.ncpath2      = 'http://opendap.deltares.nl/thredds/dodsC/opendap/rijkswaterstaat/jarkus/profiles/';
OPT.Filename2    = 'transect_r20240405.nc';

% read jarkus information
id       = ncread([OPT.ncpath2 OPT.Filename2],'id');
timej    = ncread([OPT.ncpath2 OPT.Filename2],'time');
timebj   = ncread([OPT.ncpath2 OPT.Filename2],'time_bounds');
rsp_lat  = ncread([OPT.ncpath2 OPT.Filename2],'rsp_lat');
rsp_lon  = ncread([OPT.ncpath2 OPT.Filename2],'rsp_lon');
lat      = (ncread([OPT.ncpath2 OPT.Filename2],'lat'))';
lon      = (ncread([OPT.ncpath2 OPT.Filename2],'lon'))';
cross    = ncread([OPT.ncpath2 OPT.Filename2],'cross_shore');
along    = ncread([OPT.ncpath2 OPT.Filename2],'alongshore');
areaname = (ncread([OPT.ncpath2 OPT.Filename2],'areaname'))';
areacode = ncread([OPT.ncpath2 OPT.Filename2],'areacode');
rsp_x    = ncread([OPT.ncpath2 OPT.Filename2],'rsp_x');
rsp_y    = ncread([OPT.ncpath2 OPT.Filename2],'rsp_y');
x        = (ncread([OPT.ncpath2 OPT.Filename2],'x'))';
y        = (ncread([OPT.ncpath2 OPT.Filename2],'y'))';

% changed back to old tideinfo, for the moment 
mean_high_water = ncread([OPT.ncpath2 OPT.Filename2],'mean_high_water');
mean_low_water  = ncread([OPT.ncpath2 OPT.Filename2],'mean_low_water');

time_year = str2num(datestr(timej+OPT.refdatenum,'yyyy')); %#ok<ST2NM>

transectid = id - areacode*1e6;

% standlijnen.nc for info meanhigwater, meanlowwater, dunefoot 1843-1998
OPT.ncpath3   = 'http://opendap.deltares.nl/thredds/dodsC/opendap/rijkswaterstaat/strandlijnen/';
%OPT.ncpath3   = 'c:\Users\santinel\Documents\openearthrawdata\rijkswaterstaat\strandlijnen\';
OPT.Filename3 = 'strandlijnen.nc';

% read from 1843-1998 strandlijnen.nc file
id_STRAND      = ncread([OPT.ncpath3 OPT.Filename3],'id');
time_STRAND    = ncread([OPT.ncpath3 OPT.Filename3],'time'); % year
MHW_STRAND     = (ncread([OPT.ncpath3 OPT.Filename3],'MHW'))';
MHW_STRANDx    = (ncread([OPT.ncpath3 OPT.Filename3],'MHW_x'))';
MHW_STRANDy    = (ncread([OPT.ncpath3 OPT.Filename3],'MHW_y'))';
% LAT and LON are reversed in strandlijnen.nc
MHW_STRANDlat  = (ncread([OPT.ncpath3 OPT.Filename3],'MHW_lon'))';
MHW_STRANDlon  = (ncread([OPT.ncpath3 OPT.Filename3],'MHW_lat'))';

MLW_STRAND     = (ncread([OPT.ncpath3 OPT.Filename3],'MLW'))';
MLW_STRANDx    = (ncread([OPT.ncpath3 OPT.Filename3],'MLW_x'))';
MLW_STRANDy    = (ncread([OPT.ncpath3 OPT.Filename3],'MLW_y'))';
% LAT and LON are reversed in strandlijnen.nc
MLW_STRANDlat  = (ncread([OPT.ncpath3 OPT.Filename3],'MLW_lon'))';
MLW_STRANDlon  = (ncread([OPT.ncpath3 OPT.Filename3],'MLW_lat'))';

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

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

[idIDa idIDb]       = ismember(id,id_STRAND); % 324 common ones of 334 from strandlijnen
[timeIDa timeIDb]   = ismember(time_year, time_STRAND_year); % 34 common years

% leave out the part after 1964 (1964 included)
MHW_STRAND(:,timeIDb(1):end) = [];    MLW_STRAND(:,timeIDb(1):end) = [];
MHW_STRANDx(:,timeIDb(1):end) = [];   MLW_STRANDx(:,timeIDb(1):end) = [];
MHW_STRANDy(:,timeIDb(1):end) = [];   MLW_STRANDy(:,timeIDb(1):end) = [];
MHW_STRANDlat(:,timeIDb(1):end) = []; MLW_STRANDlat(:,timeIDb(1):end) = [];
MHW_STRANDlon(:,timeIDb(1):end) = []; MLW_STRANDlon(:,timeIDb(1):end) = [];

% Preallocation and definition
MHW_tot     = NaN(length(id),length(time_totaldatenum)); MLW_tot     = NaN(length(id),length(time_totaldatenum));
MHW_x_tot   = NaN(length(id),length(time_totaldatenum)); MLW_x_tot   = NaN(length(id),length(time_totaldatenum));
MHW_y_tot   = NaN(length(id),length(time_totaldatenum)); MLW_y_tot   = NaN(length(id),length(time_totaldatenum));
MHW_lat_tot = NaN(length(id),length(time_totaldatenum)); MLW_lat_tot = NaN(length(id),length(time_totaldatenum));
MHW_lon_tot = NaN(length(id),length(time_totaldatenum)); MLW_lon_tot = NaN(length(id),length(time_totaldatenum));

% fill the variables with old values from strandlijnen.nc
for ii = 1:length(id)
        if idIDa(ii)
            MHW_tot(ii,1:(timeIDb(1)-1))     = MHW_STRAND(idIDb(ii),:);
            MHW_x_tot(ii,1:(timeIDb(1)-1))   = MHW_STRANDx(idIDb(ii),:);
            MHW_y_tot(ii,1:(timeIDb(1)-1))   = MHW_STRANDy(idIDb(ii),:);
            MHW_lat_tot(ii,1:(timeIDb(1)-1)) = MHW_STRANDlat(idIDb(ii),:);
            MHW_lon_tot(ii,1:(timeIDb(1)-1)) = MHW_STRANDlon(idIDb(ii),:);
            
            MLW_tot(ii,1:(timeIDb(1)-1))     = MLW_STRAND(idIDb(ii),:);
            MLW_x_tot(ii,1:(timeIDb(1)-1))   = MLW_STRANDx(idIDb(ii),:);
            MLW_y_tot(ii,1:(timeIDb(1)-1))   = MLW_STRANDy(idIDb(ii),:);
            MLW_lat_tot(ii,1:(timeIDb(1)-1)) = MLW_STRANDlat(idIDb(ii),:);
            MLW_lon_tot(ii,1:(timeIDb(1)-1)) = MLW_STRANDlon(idIDb(ii),:);
        end
end

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


% read from .mat file created with MW_Jarkus.m
load('..\raw\MHW.mat');
load('..\raw\MLW.mat');

% fill the variables with new values (intersection from transect.nc profiles)
MHW_tot(:,timeIDb(1):length(time_totaldatenum))     = MHW_cross;
MHW_x_tot(:,timeIDb(1):length(time_totaldatenum))   = MHW_x;
MHW_y_tot(:,timeIDb(1):length(time_totaldatenum))   = MHW_y;
MHW_lat_tot(:,timeIDb(1):length(time_totaldatenum)) = MHW_lat;
MHW_lon_tot(:,timeIDb(1):length(time_totaldatenum)) = MHW_lon;

MLW_tot(:,timeIDb(1):length(time_totaldatenum))     = MLW_cross;
MLW_x_tot(:,timeIDb(1):length(time_totaldatenum))   = MLW_x;
MLW_y_tot(:,timeIDb(1):length(time_totaldatenum))   = MLW_y;
MLW_lat_tot(:,timeIDb(1):length(time_totaldatenum)) = MLW_lat;
MLW_lon_tot(:,timeIDb(1):length(time_totaldatenum)) = MLW_lon;

%% assembly variables from strandlijnen.nc and .mat files all together
s.mean_high_water       = mean_high_water;
s.mean_high_water_cross = MHW_tot;
s.mean_high_water_x     = MHW_x_tot;
s.mean_high_water_y     = MHW_y_tot;
s.mean_high_water_lat   = MHW_lat_tot;
s.mean_high_water_lon   = MHW_lon_tot;
s.mean_low_water        = mean_low_water;
s.mean_low_water_cross = MLW_tot;
s.mean_low_water_x     = MLW_x_tot;
s.mean_low_water_y     = MLW_y_tot;
s.mean_low_water_lat   = MLW_lat_tot;
s.mean_low_water_lon   = MLW_lon_tot;
% include relevant information from Jarkus
s.id   = id;   
s.cross = cross; 
s.along = along;
s.time        = time_totaldatenum - datenum('1970','yyyy');
s.time_bounds = ([time_totaldatenum_beg - datenum('1970','yyyy'), time_totaldatenum_end - datenum('1970','yyyy')])';
s.areaname = areaname'; 
s.areacode = areacode;
% add for convenience
s.jarkusyear = time_total; 

mwl2netcdf(OPT.outputfile, s,'processing_level', 'final');
