%% root mean square error

% Verify the calculation methods for the dunefoot position per
% area.

df_ver6 = load('D:\diama_ey\Desktop\Scripts Elly\new methodology\datasets\DF_combined_methods_thr1_0_001.mat');

% Get data from jarkus
url             = jarkus_url;
id              = nc_varget(url,'id');
areacode        = nc_varget(url,'areacode');
areaname        = nc_varget(url,'areaname');
lat             = nc_varget(url,'lat',[0,0],[-1,-1]);
lon             = nc_varget(url,'lon',[0,0],[-1,-1]);
xRSP            = nc_varget(url,'cross_shore');

% Retrieve data dunefoot

url_df = 'http://opendap.deltares.nl/thredds/dodsC/opendap/rijkswaterstaat/strandlijnen/strandlijnen.nc';

year        = round(1970 + nc_varget(url_df,'time')/365);
yrs         = 1965;
start_year  = find(year == yrs)-1;
yrs2        = year(end);
n_year      = yrs2-yrs+1;

id_DF   = nc_varget(url_df,'id');
DF      = nc_varget(url_df,'DF',[0 start_year],[-1 n_year]);
DF_x    = nc_varget(url_df,'DF_x',[0 start_year],[-1 n_year]);
DF_y    = nc_varget(url_df,'DF_y',[0 start_year],[-1 n_year]);
DF_lat  = nc_varget(url_df,'DF_lon',[0 start_year],[-1 n_year]);
DF_lon  = nc_varget(url_df,'DF_lat',[0 start_year],[-1 n_year]);

[idID_STa, idID_STb] = ismember(id,id_DF);

idID_STb(idID_STb == 0) = [];               %a = [9 9 8 8 7 7 7 6 6 6 5 5 4 4 2 1 1 1]
                                            %b = [1 1 1 3 3 3 3 3 4 4 4 4 4 9 9 9]
                                            %[lia1,locb1] = ismember(a,b)
                                            %returns
                                            %lia1 = [1 1 0 0 0 0 0 0 0 0 0 0 1 1 0 1 1 1]
                                            %locb1 = [14 14 0 0 0 0 0 0 0 0 0 0 9 9 0 1 1 1]

DF_filt(:,:)     = DF(idID_STb,:);
DF_filt_x(:,:)   = DF_x(idID_STb,:);
DF_filt_y(:,:)   = DF_y(idID_STb,:);
DF_filt_lat(:,:) = DF_lat(idID_STb,:);
DF_filt_lon(:,:) = DF_lon(idID_STb,:);

% Calculate X_rsp positions for the strandlijnen

[~, n] = size(DF);

DF_cross = zeros(size(DF_filt));

latdum = lat(idID_STa,:);
dd=[];
for i = 1:n;
    for j = 1:length(idID_STb)
        if ~isnan(DF_filt_lat(j,i))
            dummy    = abs(latdum(j,:) - DF_filt_lat(j,i));
            dd=[dd; min(dummy)];
            id_cross = find(dummy == min(dummy));
            DF_cross(j,i) = xRSP(id_cross);
        else
            DF_cross(j,i) = NaN;
        end
    end
end

% Jarkus dunefoot data
url_jks = 'http://opendap.deltares.nl/thredds/dodsC/opendap/rijkswaterstaat/DuneFoot/DF.nc';
id_DF2   = nc_varget(url_jks,'id');
[idID_STa, idID_STb] = ismember(id_DF2,id_DF);

for i = 1:n
    
    year            = round(1970 + nc_varget(url_df,'time')/365);
    yrs             = 1964 + i;
    start_year      = find(year == yrs)-1;
    yrs2            = year(end);
    n_year          = yrs2-yrs+1;
    
    df_uppermkl     = nc_varget(url_jks,'dune_foot_upperMKL_cross',[start_year,0],[1,-1]);
    df_threenap     = nc_varget(url_jks,'dune_foot_threeNAP_cross',[start_year,0],[1,-1]);
    
    DF_upperMKL(:,i)     = df_uppermkl(idID_STa);
	DF_threeNAP(:,i)     = df_threenap(idID_STa);
    
end

%% second derivative method dunefoot data

for i = 1:n
          
    xdf_m1(:,i)     = df_ver6.xDF(i,idID_STa);
    zdf_m1(:,i)     = df_ver6.zDF(i,idID_STa);
    
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% choose the dunefoot
dfcalc      = xdf_m1;
dfcalc_3    = DF_threeNAP;    %% xdf_m1 % DF_upperMKL; DF_threeNAP;  % cross shore direction 
%dfcalc_upper= DF_upperMKL;

% Look at data per pre-defined area
[idID_STa, idID_STb] = ismember(id,id_DF);
areacode2 = areacode(idID_STa);
areaname2 = areaname(idID_STa,:);

[kv, m, dummy]  = unique(areacode2);
n               = vertcat(m(2:end)-1,length(areacode2));
x = -1000:100:1000;
y = x;

%% For the entire Dutch coast
% new method
sq_error    = (DF_cross-dfcalc).^2;
se          = sq_error(~isnan(sq_error)); 
rmse_total  = sqrt(mean(se));

% old method
%dfcalc_3(isnan(dfcalc))= NaN; %m :: same results
sq_error_old    = (DF_cross-dfcalc_3).^2;
se_old          = sq_error_old(~isnan(sq_error_old)); 
rmse_total3     = sqrt(mean(se_old));

% % Upper MKL
% sq_error_upper  = (DF_cross-dfcalc_upper).^2;
% se_upper        = sq_error_upper(~isnan(sq_error_upper)); 
% rmse_total_upper= sqrt(mean(se_upper))

%% rmse calculations

rmse        = NaN(12,1);
rmse3       = NaN(12,1);
% rmse_upper  = NaN(12,1);
for i=1:12 % different areas

    DF_msrd  = DF_cross(m(i):n(i),:);  % matrix for one area only 
    DF_comp  = dfcalc(m(i):n(i),:);
    DF_comp3 = dfcalc_3(m(i):n(i),:);
%    DF_comp_upper =dfcalc_upper(m(i):n(i),:);
       
    sq_error    = (DF_msrd-DF_comp).^2;
    sq_error3   = (DF_msrd-DF_comp3).^2;
%    sq_error_upper= (DF_msrd-DF_comp_upper).^2;
    
    se  = sq_error(~isnan(sq_error)); 
    se3 = sq_error3(~isnan(sq_error3)); 
%    se_upper=sq_error_upper(~isnan(sq_error_upper)); 
    
    rmse(i,1)   = sqrt(mean(se));
    rmse3(i,1)  = sqrt(mean(se3));
%    rmse_upper(i,1) = sqrt(mean(se_upper))
end
