clc; clear all;close all;
load polynyas_data% Data clasified inside the polynyas plus the ocntours for the polynyas

% Creates an structure of strings to call data from different years
years=(2004:2019)';
labyear=strcat('y_',num2str(years));
%
%% To plot the polynyas contours
figure('units','normalized','outerposition',[0 0 0.8 1]);
% pcolorps(Ybat,Xbat,-batim);hold on;shading flat;
for i=1:length(years) % loop into the 16 years
    for j=1:length(poly_data.(labyear(i,:)).name)% loop into the 18 polynyas
        % The next line plots the polynya polygons that are used to help define the polynyas contours
        plotps(poly_data.(labyear(1,:)).polygon{j}(:,2),poly_data.(labyear(i,:)).polygon{j}(:,1),'g');hold on
        % The next loop plots the polynya contours for each year
        
        pol=poly_contour.(labyear(i,:)){j};% "pol" gives the polynya contour for each month. We'll plot here the ones for e.g. August
        if any(~isnan(pol{8}(:,1)))
            plotps(pol{8}(:,1),pol{8}(:,2),'r');colormap(flipud(bone(10)));caxis([0 2000]);hold on
        end
    end
end
xlim([1e6 3e6])
ylim([-2.5e6 2.2e6])
%% To get the data inside polynyas
for i=1:length(years) % loop into the 16 years
    
    for j=1:length(poly_data.(labyear(i,:)).name)% loop into the 18 polynyas
        lon=poly_data.(labyear(i,:)).lon{j};
        if any (~isnan(lon))% when there is data
            lat=poly_data.(labyear(i,:)).lat{j};% latitude (°North)
            time=poly_data.(labyear(i,:)).time{j};% longitude (° East)
            mld=poly_data.(labyear(i,:)).mld{j};% mixed layer depth(m)
            sal=poly_data.(labyear(i,:)).sal{j};% Practical Salinity (in EOS80)
            temp=poly_data.(labyear(i,:)).temp{j};% Potential temperature (in EOS80)
            N2=poly_data.(labyear(i,:)).N2{j};% Stratification (Brunt-Vaïsälä frecquency)
            p=poly_data.(labyear(i,:)).pres{j};% Pressure
            tag=poly_data.(labyear(i,:)).tag{j};% Tag number
            ro=poly_data.(labyear(i,:)).rho{j}; % Potential density anomaly (potential density-1000 kg m-3)
           
            %%% to plot any variable for a given polynya and year
            % example: Cape Darnley polynya 
            % polynyas names are contained inside the estructure for each
            % year: e.g. poly_data.(labyear(1)).name
            polynya_name='Cape Darnley';
            if strcmp(polynya_name,poly_data.(labyear(i,:)).name{j})
                auxt=repmat(time,2000,1);
                
                figure                    
                s1=subplot(2,1,1);
                scatter(auxt(:),p(:),10,sal(:),'filled');cc=colorbar;ylabel(cc,'Salinity');colormap(s1,parula);hold on
                plot(time,mld,'ok','markersize',2,'MarkerFaceColor','k')
                datetick('x','dd/mm/yy');ylabel('Depth')
                set(gca,'ydir','reverse');xlim([min(auxt(:),[],'omitnan') max(auxt(:),[],'omitnan')])
                
                s2=subplot(2,1,2)  ;
                scatter(auxt(:),p(:),10,temp(:),'filled');cc=colorbar;ylabel(cc,'\theta (\circC)'); colormap(s2,autumn);hold on
                plot(time,mld,'ok','markersize',2,'MarkerFaceColor','k')
                set(gca,'ydir','reverse')
                xlabel('Time');ylabel('Depth')
                datetick('x','dd/mm/yy')
                set(gca,'ydir','reverse');xlim([min(auxt(:),[],'omitnan') max(auxt(:),[],'omitnan')])                
                
            end
            polynya_year=[poly_data.(labyear(i,:)).name{j},'_',num2str(years(i))]% prints the polynya name and year
        end
    end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% How to read the whole dataset (only clasiffied by years but not inside polynyas)
close all
load full_seals_data

figure('units','normalized','outerposition',[0 0 0.4 1]);
cmap2=cmocean('-haline',9);
for i=1:length(sealdata.lon)% each year
    if ~isempty(sealdata.lon{i})
        sal=sealdata.sal{i};
        temp=sealdata.temp{i};
        pressure=sealdata.p{i};
        lat=sealdata.lat{i};
        lon=sealdata.lon{i};
        t=sealdata.time{i};% julian days
        
        %%% Here we plot surface salinity as an example
        scatterps(lat,lon,5,mean(sal(1:10,:),'omitnan'),'filled');hold on
        graticuleps(-80:10:-60,0:30:150)
        textps(-72,60,'60\circE','rotation',30,...% rotate it 30 degrees clockwise
            'color',0.5*[1 1 1],'verticalalignment','bottom') % make the pinning point below the text (rather than centered on it)
        textps(-72,90,'90\circE','rotation',0,...% rotate it 30 degrees clockwise
            'color',0.5*[1 1 1],'verticalalignment','bottom') % make the pinning point below the text (rather than centered on it)
        textps(-72,120,'120\circE','rotation',-30,...% rotate it 30 degrees clockwise
            'color',0.5*[1 1 1],'verticalalignment','bottom') %
        set(gca,'tickdir','out','fontsize',16,'xtick',[],'ytick',[])
        xlim([1e6 3.5e6])
        ylim([-2.5e6 3e6])
        c=colorbar;ylabel(c,'Salinity')
        colormap(cmap2)
        caxis([33 34.8])
    end
end
