%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% script provided by @Laurent Oziel to visualize binned data from 
% the dataset publushed in SEANOE under the doi:
% associated with the manuscript published in Remote Sensing of Environment
% under the doi: 
% The user must have: 
% (1) the 'color_mark.m' function in the working path,
% (2) m_map (https://www.eoas.ubc.ca/~rich/map.html), 
% (3) cmocean (https://www.eoas.ubc.ca/~rich/map.html), and 
% (4) ncload
% (https://github.com/nctoolbox/nctoolbox/blob/master/cdm/utilities/misc/ncload.m)
% (5) cptcmap
% (https://fr.mathworks.com/matlabcentral/fileexchange/28943-color-palette-tables-cpt-for-matlab)
% 
% 
% Feel free to adapt to your needs, provided without any warranty
% contact: laurent.oziel@awi.de
% 01/04/2022
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all
close all
clc

czcs_annual = ncload('CZCS_1979_1984_ANNUAL_binned.nc');
seawifs_annual = ncload('SEAWIFS_1998_2010_ANNUAL_binned.nc');
modis_annual = ncload('MODIS_2007_2016_ANNUAL_binned.nc');

delta1 = log10(seawifs_annual.Chla) - log10(czcs_annual.Chla);
delta2 = log10(modis_annual.Chla) - log10(seawifs_annual.Chla);

pts= 2;

figure;
ind=find(isnan(delta1)==0);
m_proj('stereographic','lat',90,'long',-90,'radius',40);
[XX,YY]=m_ll2xy(czcs_annual.Longitudes(ind),czcs_annual.Latitudes(ind),'clip','on');
h = color_mark(XX,YY,delta1(ind),pts);
m_coast('patch',[.6 .6 .6],'edgecolor','none');
m_grid('linest','none','linewi',2,'tickdir','out','fontsize',16,'YTickLabel',[],'XTickLabel',[],'backcolor',[.3 .3 .3]);
colorbar
caxis([-.7 .7])
cmocean('balance');
title('SEAWIFS - CZCS')

figure;
ind=find(isnan(delta1)==0);
m_proj('stereographic','lat',90,'long',-90,'radius',40);
[XX,YY]=m_ll2xy(czcs_annual.Longitudes(ind),czcs_annual.Latitudes(ind),'clip','on');
h = color_mark(XX,YY,delta2(ind),pts);
m_coast('patch',[.6 .6 .6],'edgecolor','none');
m_grid('linest','none','linewi',2,'tickdir','out','fontsize',16,'YTickLabel',[],'XTickLabel',[],'backcolor',[.3 .3 .3]);
colorbar
caxis([-.7 .7])
cmocean('balance');
title('MODIS - SeaWiFS')


czcs_summer = ncload('CZCS_1979_1984_SUMMER_binned.nc');
seawifs_summer = ncload('SEAWIFS_1998_2010_SUMMER_binned.nc');
modis_summer = ncload('MODIS_2007_2016_SUMMER_binned.nc');

figure;
x = log10(modis_summer.Chla);
ind=find(isnan(x)==0);
m_proj('stereographic','lat',90,'long',-90,'radius',40);
[XX,YY]=m_ll2xy(czcs_annual.Longitudes(ind),czcs_annual.Latitudes(ind),'clip','on');
h = color_mark(XX,YY,x(ind),pts);
m_coast('patch',[.6 .6 .6],'edgecolor','none');
m_grid('linest','none','linewi',2,'tickdir','out','fontsize',16,'YTickLabel',[],'XTickLabel',[],'backcolor',[.3 .3 .3]);
colorbar
caxis([-1.1 1]);
cmap=cptcmap('seminf-haxby');
cmap=interp_cmap(cmap,201);
colormap(cmap);
title('MODIS summer')