% sat_plot_bathy.m
%========================================================================
% Bathymetry
% https://download.gebco.net/
%========================================================================

clear;

path_name_bathy='~/satellites/pacific/bathymetry/';
fname_bathy='gebco_2020.nc';

%ncdisp('gebco_2020.nc')
%lon=ncread('gebco_2020.nc','lon');
%lat=ncread('gebco_2020.nc','lat');
%h=ncread('gebco_2020.nc','elevation');

p=['h=ncread(''' path_name_bathy fname_bathy ''',''elevation'');'];
eval(p);
p=['lon=ncread(''' path_name_bathy fname_bathy ''',''lon'');'];
eval(p);
p=['lat=ncread(''' path_name_bathy fname_bathy ''',''lat'');'];
eval(p);
% whos
%  h                    14400x9600            1105920000  double
%  lat                   9600x1                    76800  double
%  lon                  14400x1                   115200  double
h=double(h);

h_10=h(1:10:end,1:10:end);
lat_10=lat(1:10:end);
lon_10=lon(1:10:end);

% Font size.
fnt_ax=10;
fnt_lbl=10;
fnt_title=10;
font_name='times';
% South Pacific area
latmin=-67;
latmax=-45;
lonmin=-175;
lonmax=-125;
%latmin=-62-3;
%latmax=-30;
%lonmin=-155-10;
%lonmax=-145+10;
fig_pos=[136.00        265.00        779.00        679.00];
fig_color=[1 1 1];
%*****************
% Calculate relation relyx = axisy/axisx (for PlotBoxAspectRatio)
%*****************
weg = 111*cos(pi*(latmin+(latmax-latmin)/2)/180);% Weight of one degree longitude
xc = (lonmax-lonmin)*weg;                        % Longitude, km
yc = (latmax-latmin)*111;                        % Latitude, km
relyx = yc/xc;
%%%set(gca,'PlotBoxAspectRatio',[1 relyx*1 1]);
figure;
set(gcf,'position',fig_pos,'color',fig_color,...
    'paperPositionMode','auto');
colormap jet;
%*****************
% Plot image.
%*****************
pcolor(lon_10,lat_10,h_10');
shading flat;
caxis([-6000 0]);
colorbar('hor');
title('Bathymetry (m)');
xlabel('Longitude','fontname',font_name,'Fontsize',fnt_lbl);
ylabel('Latitude','fontname',font_name,'Fontsize',fnt_lbl);
set(gca,'Fontsize',fnt_ax,'box','on',...
    'tickdir','out',...
    'ticklength',get(gca,'ticklength')*0.7,...
    'fontname',font_name,'Fontsize',fnt_ax);
axis([lonmin lonmax latmin latmax]);    % Range
set(gca,'PlotBoxAspectRatio',[1 relyx*1 1]);
%*****************
% Set axes position.
%*****************
ah=get(gca,'position');
set(gca,'position',[ah(1)-ah(3)*0.05 ah(2)-ah(4)*0.09 ah(3)*1.12 ah(4)*1.12]);
%print -r120 bathymetry_pacific.png -dpng

