close all; clear all; clc;

%% jarkus

disp('Retrieve JARKUS data...');

t = jarkus_transects('areacode',7,'id',7003500:7004500);

t_m = jarkus_merge(t);
t_mi = jarkus_interpolatenans(t_m, 'dim', 2);

Xj = t_mi.x;
Yj = t_mi.y;
Aj = t_mi.altitude;

%% vaklodingen

disp('Retrieve Vaklodingen data...');

Xv = []; Yv = []; Av = [];
for j = [2726] % 2928]
    for i = 119:121
        url = ['P:\mcdata\opendap\rijkswaterstaat\vaklodingen\vaklodingenKB' num2str(i) '_' num2str(j) '.nc'];

        info = nc_info(url);

        z = nc_varget(url, 'z', [info.Dimension(1).Length-1 1 1], [1 -1 -1]);
        x = nc_varget(url, 'x');
        y = nc_varget(url, 'y');

        x = x(2:end)-10;
        y = y(2:end)-10;

        [X Y] = meshgrid(x, y);

        Xv = [Xv X];
        Yv = [Yv Y];
        Av = [Av z];
    end
end

%% ahn 25

disp('Retrieve AHN 25m data...');

load('P:\1002266-xbeach\TheDayWeLostEgmond\bathymetry\AHN\cropped_egmond.mat');
As = Z;
[Xs Ys] = meshgrid(X, Y);

%% ahn 5

disp('Retrieve AHN 5m data...');

load('P:\1002266-xbeach\TheDayWeLostEgmond\bathymetry\GIS2\elevations.mat');
Ae = A / 100;

load('P:\1002266-xbeach\TheDayWeLostEgmond\bathymetry\GIS2\houses.mat');
Ah = A / 100;

load('P:\1002266-xbeach\TheDayWeLostEgmond\bathymetry\GIS2\roads.mat');
Ar = A / 100;

% determine grid
x = info.XLLCORNER+[0:info.NCOLS-1]*info.CELLSIZE;
y = info.YLLCORNER+[0:info.NROWS-1]*info.CELLSIZE;
[X Y] = meshgrid(x, y);

%% rotate grids

disp('Rotate grids...');

alpha = -8/180*pi;
R = [cos(alpha) sin(alpha); -sin(alpha) cos(alpha)];
Xr = info.XLLCORNER+R(1,1)*(X-info.XLLCORNER)+R(1,2)*(Y-info.YLLCORNER);
Yr = info.YLLCORNER+R(2,1)*(X-info.XLLCORNER)+R(2,2)*(Y-info.YLLCORNER);
Xvr = info.XLLCORNER+R(1,1)*(Xv-info.XLLCORNER)+R(1,2)*(Yv-info.YLLCORNER);
Yvr = info.YLLCORNER+R(2,1)*(Xv-info.XLLCORNER)+R(2,2)*(Yv-info.YLLCORNER);
Xjr = info.XLLCORNER+R(1,1)*(Xj-info.XLLCORNER)+R(1,2)*(Yj-info.YLLCORNER);
Yjr = info.YLLCORNER+R(2,1)*(Xj-info.XLLCORNER)+R(2,2)*(Yj-info.YLLCORNER);
Xsr = info.XLLCORNER+R(1,1)*(Xs-info.XLLCORNER)+R(1,2)*(Ys-info.YLLCORNER);
Ysr = info.YLLCORNER+R(2,1)*(Xs-info.XLLCORNER)+R(2,2)*(Ys-info.YLLCORNER);

%% modify grids

% remove measurements below water level
Ae(Ae < 3) = NaN;
%As(As < 0) = NaN;

% straighten offshore boundary
xx = [];
yy = [];
for i = 1:size(Av,1)
    j = min(find(Av(i,:) >= -15));
    
    xx(i) = Xvr(i,j);
    yy(i) = Yvr(i,j);
    
    k = min(find(Xvr(i,:) >= Xvr(i,j)-500));
    
    Av(i,1:k) = -20;
    Av(i,k:j) = interp1(Xvr(i,[k j]), Av(i,[k j]), Xvr(i,k:j));
end

%% plot

disp('Plot data...');

% plot data
figure;
pcolor(Xsr, Ysr, As); hold on;
pcolor(Xjr, Yjr, Aj); hold on;
pcolor(Xvr, Yvr, Av); hold on;
pcolor(Xr, Yr, Ae); hold on;
pcolor(Xr, Yr, Ar); hold on;
pcolor(Xr, Yr, Ah); hold on;
shading flat; axis equal; grid on;
title('Raw data');

% draw vertical comparison line
%line([1.029e5 1.029e5], [5.135e5 5.165e5], 'LineWidth', 3, 'Color', [1 0 0]);

%% create grid

disp('Create calculation grid...');

% cross-shore
y = 514800;
iv = find(abs(Yvr(:,1)-y)<10);

x = 93800;
jv = find(abs(Xvr(iv,:)-x)<10);

Tp = 18;
s = warning('off', 'all');

dxmin = 5; %5 original
X = profilegen(Av(iv,jv:end), Xvr(iv,jv:end), 20, Tp, dxmin);
warning(s(1).state, 'all');

CS = [X(2,:) X(2,end)+dxmin:dxmin:1000];

% longshore
ia = find(abs(Yr(:,1)-y)<1);

dy1 = 10;
dy2 = 20;
[f n T2] = grid_transition(dy1, dy2, 500);
T1 = cumsum([fliplr(diff(T2)) T2(1)]);

LS = Yr(ia,1)+[-2500:dy2:-1000 -1000+T1 -500+dy1:dy1:500 500+T2 1000+dy2:dy2:2500];

% mesh
[Xa Ya] = meshgrid(CS, LS);

figure;
subplot(2,2,1); plot(CS); title({'Calculation grid','cross-shore'});
subplot(2,2,2); plot(diff(CS)); title({'Calculation grid','cross-shore (diff)'});
subplot(2,2,3); plot(LS); title({'Calculation grid','long-shore'});
subplot(2,2,4); plot(diff(LS)); title({'Calculation grid','long-shore (diff)'});

%% grid data

disp('Interpolate data to calculation grid...');

figure;
title('Interpolated data on calculation grid');

Zs = griddata(Xsr, Ysr, As, Xa, Ya, 'linear'); disp('    AHN 25m data');
pcolor(Xa, Ya, Zs); hold on; shading flat; axis equal; grid on; drawnow;
Zj = griddata(Xjr, Yjr, Aj, Xa, Ya, 'linear'); disp('    JARKUS data');
pcolor(Xa, Ya, Zj); hold on; shading flat; axis equal; grid on; drawnow;
Zv = griddata(Xvr, Yvr, Av, Xa, Ya, 'linear'); disp('    Vaklodingen data');
pcolor(Xa, Ya, Zv); hold on; shading flat; axis equal; grid on; drawnow;
Ze = griddata(Xr, Yr, Ae, Xa, Ya, 'linear'); disp('    AHN 5m data (elevations)');
pcolor(Xa, Ya, Ze); hold on; shading flat; axis equal; grid on; drawnow;
Zr = griddata(Xr, Yr, Ar, Xa, Ya, 'linear'); disp('    AHN 5m data (roads)');
pcolor(Xa, Ya, Zr); hold on; shading flat; axis equal; grid on; drawnow;
Zh = griddata(Xr, Yr, Ah, Xa, Ya, 'linear'); disp('    AHN 5m data (houses)');
pcolor(Xa, Ya, Zh); hold on; shading flat; axis equal; grid on; drawnow;
%Ze = nan(size(Zv)); Zr = nan(size(Zv)); Zh = nan(size(Zv));

%% sources

disp('Determine data sources...');

S = nan(size(Xa));

ns = isnan(Zs);
nj = isnan(Zj);
nv = isnan(Zv);
ne = isnan(Ze);
nr = isnan(Zr);
nh = isnan(Zh);

S(~ns) = 10;
S(~nj) = 20;
S(~nv) = 30;
S(~ne | ~nr | ~nh) = 40;
S(~nr | ~nh) = 50;

disp('Plot data sources...');

figure;
colormap([1 0 0; 0 1 0; 0 0 1; 0 0 0]);
pcolor(Xa, Ya, S);
shading flat; axis equal; grid on;
colorbar;
title('Datasources');

%% finalize grid

% glue layers together
Za = zeros(size(Xa));
Znea = 20*ones(size(Xa));
nZs = isnan(Zs); Za(~nZs) = Zs(~nZs);
nZj = isnan(Zj); Za(~nZj) = Zj(~nZj);
nZv = isnan(Zv); Za(~nZv) = Zv(~nZv);
nZe = isnan(Ze); Za(~nZe) = Ze(~nZe);
%nZr = isnan(Zr); Za(~nZr) = Zr(~nZr); %Znea(~nZr) = 0;
nZh = isnan(Zh); Za(~nZh) = Zh(~nZh); Znea(~nZh) = 0;

ind1 = min(find(Xa(1,:) >= 98400));
ind2 = max(find(Xa(1,:) <= 103800));
Xaa = Xa(:,ind1:ind2);
Yaa = Ya(:,ind1:ind2);
Zaa = Za(:,ind1:ind2);
Zneaa = Znea(:,ind1:ind2);

% determine origin
Xori = min(min(Xaa));
Yori = min(min(Yaa));

X = Xaa-Xori;
Y = Yaa-Yori;
Z = Zaa;
Zne = Zneaa;

% extend grid in longshore direction
gs = 20; gn = 5;
X = [ones(gn,1)*X(1,:); X; ones(gn,1)*X(end,:)];
Y = [   ones(gn,1)*Y(1,:)+gs*[-gn:-1]'*ones(1,size(Y,2)); ...
        Y; ...
        ones(gn,1)*Y(end,:)+gs*[1:gn]'*ones(1,size(Y,2)) ...
    ];

Z = [ones(gn,1)*Z(1,:); Z; ones(gn,1)*Z(end,:)];
Zne = [ones(gn,1)*Zne(1,:); Zne; ones(gn,1)*Zne(end,:)];

% add polder
gs = 20; pl = 750; gn = pl/gs;
X = [X X(:,end)*ones(1,gn)+ones(size(X,1),1)*[1:gn]*gs];
Y = [Y Y(:,end)*ones(1,gn)];
Z = [Z max(-5, Z(:,end)*ones(1,gn)-ones(size(Z,1),1)*[1:gn])];
Zne = [Zne Zne(:,end)*zeros(1,gn)+20];

% add neuman proof borders
gn = 5;
Z1 = Z(1,:); Ze = Z(end,:); m1=-Inf; me=-Inf;
for i = 1:length(Z1)
    if Z1(i)>5 || i >= size(Z,2)-pl/gs;
        m1 = max(Z1(i),m1); Z(1:gn,i) = m1*ones(gn,1);
        me = max(Ze(i),me); Z(end-gn+1:end,i) = me*ones(gn,1);
    end
end

%% write data

disp('Write bathymetry to files...');

save('x.grd', 'X', '-ascii');
save('y.grd', 'Y', '-ascii');
save('bed.dep', 'Z', '-ascii');
save('ne_bed.dep', 'Zne', '-ascii');

disp('Plot written grids...');

figure;
pcolor(X, Y, Z);
shading flat; axis equal; grid on;
print bathy.png -dpng
title('Written bathymetry');

figure;
pcolor(X, Y, Zne);
shading flat; axis equal; grid on;
title('Written non-erodable layers');

%% extras

disp(['nx = ' num2str(size(X,2)-1)]);
disp(['ny = ' num2str(size(X,1)-1)]);
disp(['xori = ' num2str(Xori)]);
disp(['yori = ' num2str(Yori)]);

disp('Done');