% Map of West Antarctic Peninsula
% Mapping help from Schuyler Nardelli & Jack Conroy, March 2020
% Updated by Jessie Turner June 2021 through October 2023

clear; close all hidden;
% Load data, saved from script "calculate_bloomIndices_WAP.m"
datafolder = '/Users/jsturner/satellite_data/matlab_WAP/daily_CMEMS/';
load([datafolder 'StartDates_CMEMS_WAP_v4.mat']);
% Load map information
mapfolder = '/Users/jsturner/satellite_data/matlab_WAP/WAP_map_SchuylerJack/';
load([mapfolder 'bathy_and_coast.mat'],'coastlon','coastlat','shlflon','shlflat','lon','lat','depth'); % load Jack+Schuyler's coastline:
load([mapfolder 'WAP_regions_2022-07-15.mat']);

%% Subplot 1: Map with Bathymetry

fig = figure(1);
clf
set(gcf,'position',[158   496   829   420]);

ax(1) = subplot(1,3,1);

% Can change map limits by lat and long in the next two lines
m_proj('lambert','long',[-80 -55],'lat',[-70.05 -59.95]); hold on
plt = m_pcolor(lon,lat,-depth); shading flat;
cmap_purple = [236,231,242
    208,209,230
    166,189,219
    116,169,207
    54,144,192
    5,112,176]./255;
colormap(ax(1),cmap_purple);
caxis([0 3000]);
hold on
% thicker line
m_grid('tickdir','in','xtick',-80:10:-55,'ytick',-70:2:-60,'fontsize',14,'linewidth',2,'fontname','Arial','backcolor','none');

% Add coast
coastcolor = [.8 .8 .8];
for i=1:length(coastlon)
    p_coast = m_patch(coastlon{i},coastlat{i},coastcolor,'LineWidth',1);
end; hold on;

% Add ice shelf
shelfcolor = [1 1 1];
for i=1:length(shlflon)
    p_shelf = m_patch(shlflon{i},shlflat{i},shelfcolor,'LineWidth',1);
end

cb = colorbar;
cb.Location  = 'southoutside';
cb.LineWidth = 1.5;
cb.FontSize = 16;
cb.Label.String = 'Depth (m)';
cb.Label.FontSize = 16;
cb.Position = [0.1    0.225   0.225    0.03];

% Plot ecoregions
plt1 = m_plot(lon_sos,lat_sos,'-','linewidth',3);
plt2 = m_plot(lon_mis,lat_mis,'-','linewidth',3);
plt3 = m_plot(lon_nos,lat_nos,'-','linewidth',3);
plt4 = m_plot(lon_miz,lat_miz,'-','linewidth',3);
plt5 =  m_plot(lon_spf,lat_spf,'-','linewidth',3);

% Number ecoregions
str1 = {'1','2','3','4','5'};
str_x = [-77.69,-77.69,-61.4,-70.23,-74.9];
str_y = [-61.6,-66,-63.23,-66.35,-68.56];
m_text(str_x,str_y,str1,'fontsize',18, 'fontweight','bold');

pos = get(gca,'position');
set(gca,'pos',[pos(1)-0.06 pos(2)-0.02 pos(3)+0.08 pos(4)+0.08]);

%% Inset map to show where on Antarctica
subplot(1,3,3)
m_proj('stereographic','lat',-90,'lon',0,'radius',35);
m_coast('patch',coastcolor);
m_grid('xtick',12,'tickdir','in','ytick',[-90:10:-60],'fontsize',...
    10,'linewidth',2,'fontname','Arial','backcolor','none',...
    'XAxisLocation','top');
hold on
bndry_lon=[-80 -80 -75 -70  -65 -60  -55 -55 -73  -80];
bndry_lat=[-70 -60 -60 -60 -60 -60  -60 -70 -70  -70];
m_line(bndry_lon,bndry_lat,'linewi',2,'color','k');     % Area outline ...
    m_hatch(bndry_lon,bndry_lat,'single',30,5,'color','k'); % ...with hatching added.
pos = get(gca,'position');
set(gca,'position',[pos(1)+0.08 pos(2)+0.02 pos(3)-0.02 pos(4)-0.02]);

%% Subplot 2: Map with bloom timing
% Average over all years for all bloom indices
% Long-term means of all indices
bi_avg_thres = nanmean(bi_thres,3);

% Plot each bloom index long-term average start date (map)
ax(2) = subplot(1,3,2);

% Set up ticks for day of year
dates_ticks = [244 274 305 335 366 397 425];
datestr_bi = list_doy_strings(dates_ticks);

m_proj('lambert','long',[-80 -55],'lat',[-70.05 -59.95]); hold on
m_pcolor(lon_chl,lat_chl,bi_avg_thres); shading flat;
colormap(ax(2),turbo(30));
hold on
m_grid('tickdir','in','xtick',-80:10:-55,'ytick',-70:2:-60,'fontsize',14,'linewidth',2,'fontname','Arial','backcolor','none');
caxis([dates_ticks(1) dates_ticks(end)]);
cb2 = colorbar;
cb2.Location = 'southoutside';
cb2.Label.String = 'Start Date';
cb2.LineWidth = 1.5;
cb2.Ticks = dates_ticks;
cb2.TickLabels = {'S','O','N','D','J','F','M'};
cb2.FontSize = 16;
cb2.Position = [0.462    0.225    0.225    0.03];

% Add coast
coastcolor = [.8 .8 .8];
for i=1:length(coastlon)
    p_coast = m_patch(coastlon{i},coastlat{i},coastcolor,'LineWidth',1);
end; hold on;
% Add ice shelf
shelfcolor = [1 1 1];
for i=1:length(shlflon)
    p_shelf = m_patch(shlflon{i},shlflat{i},shelfcolor,'LineWidth',1);
end

% Plot regions %%% NEW as of July 15, 2022
plt1 = m_plot(lon_sos,lat_sos,'-','linewidth',3);
plt2 = m_plot(lon_mis,lat_mis,'-','linewidth',3);
plt3 = m_plot(lon_nos,lat_nos,'-','linewidth',3);
plt4 = m_plot(lon_miz,lat_miz,'-','linewidth',3);
plt5 =  m_plot(lon_spf,lat_spf,'-','linewidth',3);

% str1 = {'1','2','3','4','5'};
% str_x = [-77.69,-77.69,-61.4,-70.23,-74.9];
% str_y = [-61.6,-66,-63.23,-66.35,-68.56];
% m_text(str_x,str_y,str1,'fontsize',18, 'fontweight','bold','color','w');
pos = get(gca,'position');
set(gca,'pos',[pos(1)+0.02 pos(2)-0.02 pos(3)+0.08 pos(4)+0.08]);

%%
sublabel(ax,40,10,'fontsize',18);

%% Save - Make sure to use the correct date.
print(gcf,'Fig_1_Map_WAP.png','-dpng','-r300');
print(gcf,'Fig_1_Map_WAP.eps','-depsc');

