
%
% Plots coast and bathymetry.
% Bathymetry is plotted with contoured lines
% 

function [hCoast, hBath, hRecif] = trace_coast_bathy(lonlat_minmax,colCoast,coast_bathy_region,ibath,levels)

if ~exist('coast_bathy_region','var')
    coast_bathy_region = 'med';
end

if ~exist('ibath','var')
    ibath = 0;
end

if ~exist('levels','var')
    levels = [-100 -500 -1000 -2000];
end


%% Bathymetry
if ibath == 1 && ~exist('xbath','var')
    if exist('coast_bathy_region','var')
        load(['bathyrelief_' coast_bathy_region '.mat']);
    else
        load bathyrelief_med.mat;
    end
    xbath = bathy(:,1); %#ok<NODEF>
    ybath = bathy(:,2);
    zbath = bathy(:,3);
    valid = (xbath >= lonlat_minmax(1) & xbath <= lonlat_minmax(2)) & ...
            (ybath >= lonlat_minmax(3) & ybath <= lonlat_minmax(4));
    xbath = xbath(valid);
    ybath = ybath(valid);
    zbath = zbath(valid);
    xbath_un = unique(xbath);
    ybath_un = unique(ybath);
    xbath = reshape(xbath,[length(xbath_un) length(ybath_un)]);
    ybath = reshape(ybath,[length(xbath_un) length(ybath_un)]);
    zbath = reshape(zbath,[length(xbath_un) length(ybath_un)]);
    clear bathy xbath_un ybath_un
    
    [cLines, hBath] = contour(xbath,ybath,zbath,levels,'Color',[0.5 0.5 0.5],'Linewidth',0.5);
    clabel(cLines,hBath,'LabelSpacing',500,'FontSize',9,'Color',[0.5 0.5 0.5]);
    hold on;
end


%% Coastline
if ~exist('xcot','var')
    if exist('coast_bathy_region','var')
        load(['coast_' coast_bathy_region '.mat']);
    else
        load coast_med.mat;
    end
    xcoast = coast(:,1); %#ok<NODEF>
    ycoast = coast(:,2);
end

hCoast = plot(xcoast,ycoast,colCoast,'Linewidth',1);


%% Recifs
if strcmp(coast_bathy_region,'centralpacific')
    load recif_centralpacific.mat;
    xrecif = recif(:,1); %#ok<NODEF>
    yrecif = recif(:,2);
    valid = (xrecif >= lonlat_minmax(1) & xrecif <= lonlat_minmax(2)) & ...
            (yrecif >= lonlat_minmax(3) & yrecif <= lonlat_minmax(4));
    xrecif = xrecif(valid);
    yrecif = yrecif(valid);
    
    hRecif = plot(xrecif,yrecif,[colCoast '.'],'MarkerSize',5);
end


%% Axes limits and aspectratio
set(gca,'Xlim',lonlat_minmax(1:2), ...
        'Ylim',lonlat_minmax(3:4),'Linewidth',1);
% set(gca,'Plotboxaspectratio',[1 1*(lonlat_minmax(4)-lonlat_minmax(3))/(lonlat_minmax(2)-lonlat_minmax(1))/ ...
%                                 cos( 0.5*(lonlat_minmax(3)+lonlat_minmax(4))*pi/180 ) 1]);
[dX, dY] = lonlat2km(lonlat_minmax(1),lonlat_minmax(3),lonlat_minmax(2),lonlat_minmax(4));
set(gca,'Plotboxaspectratio',[1 abs(dY/dX) 1]);
box on;
axis vis3d;

set(gca,'TickDir','out');
xlabel('Lon (deg E)');
ylabel('Lat (deg N)');
