% script per maneggiare i dati EMODnet
% il file Mean_depth_Emodnet.dat viene da un dowload file csv di EMODNET
% che è stato modificato, togliendo i Nan (sost con 0) e togliendo le
% virgole (sost con gli spazi)
% -- il codice interpola i dati di batimetria 
% -- due esempi di griglia raffinata e più rada (coarse)
% -- coarse grid viene anche fatto uno smoothing e stampato in output

clear 

dataEmodnet = load("Mean_depth_Emodnet.dat");
np = length(dataEmodnet(:,1));

last   = np;
coordx = zeros(last,1); 
coordy = zeros(last,1); 
coordz = zeros(last,1); 

for i=1:last
    lat = dataEmodnet(i,1);
    lon = dataEmodnet(i,2);
    h   = dataEmodnet(i,3);
    [x, y, z] = geodetic2ecef(lat, lon, h);

    coordx(i,1) = x;
    coordy(i,1) = y;
    coordz(i,1) = z;
end

% raw data da EMODnet
figure(1)
scatter3(coordx,coordy,coordz,".")

% Remove outlier
idx = isoutlier(z);
x = x(~idx);
y = y(~idx);
z = z(~idx);

%% Griglia - fine
xq = linspace(min(coordx), max(coordx), 2000);
yq = linspace(min(coordy), max(coordy), 2000);
[Xq, Yq] = meshgrid(xq, yq);

% --> Interpolazione
F  = scatteredInterpolant(coordx, coordy, coordz, 'nearest', 'none');
Zq = F(Xq, Yq);

% --> Plot
figure(2)
surf(Xq, Yq, Zq)
shading interp



%% Griglia 2 - coarse
npunti = 250;
%xcq = linspace(min(coordx), max(coordx), npunti);
%ycq = linspace(min(coordy), max(coordy), npunti);
xcq = linspace(coordx(1), coordx(end), npunti);
ycq = linspace(coordy(1), coordy(end), npunti);
[Xcq, Ycq] = meshgrid(xcq, ycq);

% --> Interpolazione
F  = scatteredInterpolant(coordx, coordy, coordz, 'nearest', 'none');
Zcq = F(Xcq, Ycq);

% --> Smooth
Zcsmooth = movmean(Zcq, 5, 1);
Zcsmooth = movmean(Zcsmooth, 5, 2);

% --> Plot
figure(3)
surf(Xcq, Ycq, Zcsmooth)
shading interp

Zcsmooth = -Zcsmooth;

% ---> print out 
fid = fopen('Trieste_coast.bot','w');
for j = 1:size(Zcsmooth,1) 
    fprintf(fid,'%f ', Zcsmooth(j,1:end-1));  
    fprintf(fid,'%f\n', Zcsmooth(j,end));    
end
fclose(fid);

% check
load("Trieste_coast.bot")
figure(4)
surf(Trieste_coast)



% ===== FUNZIONE =====
function [x, y, z] = geodetic2ecef(lat, lon, h)
    lat = deg2rad(lat);
    lon = deg2rad(lon);

    a = 6378137.0;
    f = 1 / 298.257223563;
    e2 = f * (2 - f);

    N = a ./ sqrt(1 - e2 .* sin(lat).^2);

    x = (N + h) .* cos(lat) .* cos(lon);
    y = (N + h) .* cos(lat) .* sin(lon);
    z = h;
end
