% routine to get MHW and MLW as mat files.
%
% Intersection for MHW and MLW is defined as the most seaward point of the
% profile that has a descending profile, and save mat files`
%
% $Id$
% $Date$
% $Author$
% $Revision$
% $HeadURL$
% $Keywords$ 

%% settings
% jarkus database - dataset 1965 - ...
OPT.ncpath2      = 'C:\Users\fuentesm\Marilu\Deltares\Projects\Coastviewer\jarkus\scripts\';
OPT.Filename2    = 'transect_r20240405.nc';
OPT.refdatenum   = datenum(1970,1,1);

% this will be temporarly adapted to use the old tideinfo values (as the once in jarkus2017 file)

%% jarkus database - dataset 1965 - ...
%% Read data
id       = ncread([OPT.ncpath2 OPT.Filename2],'id');
MHWm     = ncread([OPT.ncpath2 OPT.Filename2],'mean_high_water');
MLWm     = ncread([OPT.ncpath2 OPT.Filename2],'mean_low_water');
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');
time     = ncread([OPT.ncpath2 OPT.Filename2],'time');

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

timenum = OPT.refdatenum+ncread([OPT.ncpath2 OPT.Filename2],'time');

%% Evaluation of MHW

% counters
[kv, m, dummy] = unique(areacode);
% n          = vertcat(1, m(1:end-1)+1);
n          = vertcat(m(2:end)-1,length(id));

% definition
% MHWm & MLWm are the same as in the source file https://svn.oss.deltares.nl/repos/openearthrawdata/trunk/rijkswaterstaat/jarkus/raw/info/TIDEINFO.txt
MHW_cross = NaN(length(areacode),length(time));
MHW_lat   = NaN(length(areacode),length(time)); 
MHW_lon   = NaN(length(areacode),length(time));
MHW_rsp   = NaN(length(areacode),length(time));

%% MHW
for i = 1:length(kv) 
    %num_tr  = m(i)-n(i)+1;
    num_tr  = n(i)-m(i)+1;

    for j = 1:num_tr % loop for transect
        % altitude = nc_varget([OPT.ncpath2 OPT.Filename2],'altitude',[0 m(i)+j-2 0],[-1 1 -1]);
        altitude = ncread([OPT.ncpath2 OPT.Filename2],'altitude',[1 m(i)+j-1 1],[Inf 1 Inf]);
        altitude = squeeze(altitude)';
        disp(['MHW: transect ' num2str(j) ' of area ' num2str(kv(i)) ...
            ' (' num2str(m(i)+j-1) '/' num2str(length(id)) ')']);
        
        for k = 1:length(time) % loop for year
            if ~isempty(altitude(k,~isnan(altitude(k,:)))) % altitude measurement exists that year
                % disp(num2str(datestr(timenum(k))))
                 % intersection, crossdir=2 is descending; *1_new_out is
                 % the new array with intersection points
                [x_int z_int x1_new_out z1_new_out dum1_new_out dum2_new_out crossdir] = ...
                    findCrossings((cross(~isnan(altitude(k,:))))',altitude(k,~isnan(altitude(k,:))),...
                    [cross(1) cross(end)],[MHWm(m(i)+j-1) MHWm(m(i)+j-1)]);

                if any(crossdir == 2 | crossdir == 1)
                    MHWxintersection = max(x_int(crossdir == 2 | crossdir == 1));
                else
                    MHWxintersection = NaN;
                end
                
                MHW_cross(m(i)+j-1,k) = MHWxintersection;
                
                if ~isnan(MHWxintersection) % is it is non-NaN, calculate WGS84 coords.
                    MHWlatdum = interp1(cross',lat(m(i)+j-1,:),MHWxintersection);
                    MHWlondum = interp1(cross',lon(m(i)+j-1,:),MHWxintersection);
                    
                    MHW_lat(m(i)+j-1,k) = MHWlatdum; 
                    MHW_lon(m(i)+j-1,k) = MHWlondum;
                else
                    MHWlatdum = NaN; 
                    MHWlondum = NaN;
                    
                    MHW_lat(m(i)+j-1,k) = MHWlatdum; 
                    MHW_lon(m(i)+j-1,k) = MHWlondum;
                end
                
                clear MHWxintersection
            end
        end
    end

end

% MHW in RD coordinates
[MHW_x, MHW_y] = convertCoordinates(MHW_lon,MHW_lat,'CS1.code',4326,'CS2.code',28992);

% remove before saving new one
if exist('..\raw\MHW.mat', 'file')==2
  delete('..\raw\MHW.mat');
end

% save
save('..\raw\MHW.mat','MHW_cross', 'MHW_rsp', 'MHW_lat', 'MHW_lon', 'MHW_x', 'MHW_y');

%% Evaluation of MLW

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

% definition
MLW_cross = NaN(length(areacode),length(time));
MLW_lat   = NaN(length(areacode),length(time)); 
MLW_lon   = NaN(length(areacode),length(time));
MLW_rsp   = NaN(length(areacode),length(time));

%% MLW
for i = 1:length(kv)
    %num_tr  = m(i)-n(i)+1;
    num_tr  = n(i)-m(i)+1;

    for j = 1:num_tr % loop for transect
        % altitude = nc_varget([OPT.ncpath2 OPT.Filename2],'altitude',[0 m(i)+j-2 0],[-1 1 -1]);
        altitude = ncread([OPT.ncpath2 OPT.Filename2],'altitude',[1 m(i)+j-1 1],[Inf 1 Inf]);
        altitude = squeeze(altitude)';
        disp(['MLW: transect ' num2str(j) ' of area ' num2str(kv(i)) ...
            ' (' num2str(m(i)+j-1) '/' num2str(length(id)) ')']);
        
        for k = 1:length(time) % loop for year
            if ~isempty(altitude(k,~isnan(altitude(k,:)))) % altitude measurement exists that year
                
                 % intersection, crossdir=2 is descending; *1_new_out is
                 % the new array with intersection points
                [x_int z_int x1_new_out z1_new_out dum1_new_out dum2_new_out crossdir] = ...
                    findCrossings((cross(~isnan(altitude(k,:))))',altitude(k,~isnan(altitude(k,:))),...
                    [cross(1) cross(end)],[MLWm(m(i)+j-1) MLWm(m(i)+j-1)]);

                if any(crossdir == 2 | crossdir == 1)
                    MLWxintersection = max(x_int(crossdir == 2 | crossdir == 1));
                else
                    MLWxintersection = NaN;
                end
                
                MLW_cross(m(i)+j-1,k) = MLWxintersection;
                
                if ~isnan(MLWxintersection) % is it is non-NaN, calculate WGS84 coords.
                    MLWlatdum = interp1(cross',lat(m(i)+j-1,:),MLWxintersection);
                    MLWlondum = interp1(cross',lon(m(i)+j-1,:),MLWxintersection);
                    
                    MLW_lat(m(i)+j-1,k) = MLWlatdum; 
                    MLW_lon(m(i)+j-1,k) = MLWlondum;
                else
                    MLWlatdum = NaN; 
                    MLWlondum = NaN;
                    
                    MLW_lat(m(i)+j-1,k) = MLWlatdum; 
                    MLW_lon(m(i)+j-1,k) = MLWlondum;
                end
                
                clear MLWxintersection
            end
        end
    end

end

% MLW in RD coordinates
[MLW_x, MLW_y] = convertCoordinates(MLW_lon,MLW_lat,'CS1.code',4326,'CS2.code',28992);

% remove before saving new one
if exist('..\raw\MLW.mat', 'file')==2
  delete('..\raw\MLW.mat');
end


% save
save('..\raw\MLW.mat','MLW_cross', 'MLW_rsp', 'MLW_lat', 'MLW_lon', 'MLW_x', 'MLW_y');


