% Use Argovis API to create a map for the number of Argo profiles in 1
% degree bins during a period of interest. Also plot some metadata.

clear all
close all


set(0, 'DefaultFigureVisible', 'on')

% select region and years of interest
dx      = 1;
dy        = 1;
x_bin     = -179.5:dx:179.5;
y_bin     = -89.5:dy:89.5;
basin_bin = 0:56;
dp        = 10;
press_bin = 5:10:2500;

years     = 2004:2019;
mm_ALL    = 1:12;
% to include months, but not an entire year
end_month = 8;
end_year  = 2019;

ArgovisNaN = -999;

% define bins edges
x_edges     = x_bin(1)-dx/2:dx:x_bin(end)+dx/2;
y_edges     = y_bin(1)-dy/2:dy:y_bin(end)+dy/2;
basin_edges = basin_bin(1)-1/2:1:basin_bin(end)+1/2;
press_edges = press_bin(1)-dp/2:dp:press_bin(end)+dp/2;

% set options for webread
options = weboptions('Timeout',40,'UserAgent', 'http://www.whoishostingthis.com/tools/user-agent/');

% other info to get
other_info                = {'PLATFORM_TYPE' 'POSITIONING_SYSTEM' 'DATA_MODE' 'dac'};

% allocate the output var
NUM                         = nan(length(x_bin),length(y_bin),length(mm_ALL)*length(years));
NUM_by_LAT_by_max_PRES_temp = nan(length(y_bin),length(press_bin),length(mm_ALL)*length(years));

% all the vars of interest
all_info = [other_info {'lon' 'lat' 'date' 'position_qc' 'pres_max_for_TEMP' ...
    'pres_min_for_TEMP' 'pres_max_for_PSAL' 'pres_min_for_PSAL'}];
n =0;
disp(datestr(now))

for yy=1:length(years)
    for imm=1:length(mm_ALL)
        if ~(years(yy)==end_year && mm_ALL(imm)==end_month)
            clear mm
            mm  = mm_ALL(imm);
            
            disp([num2str(years(yy)) ', Month ' num2str(mm)])
            % build url to query Argovis in the month/year of interest
            url = ['https://argovis.colorado.edu/selection/profiles/' ...
                num2str(mm) '/' num2str(years(yy)) '/'];
            
            % use url to query Argovis and retrieve the information in a matlab
            % structure
            clear data x y t_dnum bfr* pmax_temp basin
            
            % get data from Argovis
            tic;
            data   = webread(url,options);
            toc;
            
            % store the info of interest
            n = n+1;
            if ~isempty(data)
                %%% NOTE that since not all the profiles have the same fields,
                %%% data may be a cell of structures. In this case, 
                %%% we convert the cell to a merged
                %%% structure that includes only the fields of interest (which
                %%% should be there for every profile for this to work)
                %
                try
                    x = [data(:).lon];
                catch
                    whos data
                    clear data_select
                    for lll=1:length(data)
                        for lll_vars=1:length(all_info)
                            eval(['data_select(lll).' all_info{lll_vars} ' = ' ...
                                'data{lll}.' all_info{lll_vars} ';'])
                        end
                    end
                    clear data
                    data = data_select';
                end
                % store the location
                x = [data(:).lon];
                y = [data(:).lat];
                position_qc = [data(:).position_qc];
                
                % set to nan positions that don't have QC == 1
                x(position_qc~=1) = nan;
                y(position_qc~=1) = nan;
                %basin = [data(:).basin];
                pmax_temp = [data(:).pres_max_for_TEMP];
                
                % store the time
                bfr_t  = reshape([data(:).date],24,[])';
                t_dnum = datenum(str2num(bfr_t(:,1:4)),str2num(bfr_t(:,6:7)),...
                    str2num(bfr_t(:,9:10)),str2num(bfr_t(:,12:13)),...
                    str2num(bfr_t(:,15:16)),str2num(bfr_t(:,18:19)));
                
                % check that data are in the right month/year
                if min(t_dnum)<datenum(years(yy),mm,1) || ...
                        max(t_dnum)>datenum(years(yy),mm,eomday(years(yy),mm),23,59,59)
                    disp('Unexpected dates from the url, will stop')
                    disp(datestr(min(t_dnum)))
                    disp(datestr(max(t_dnum)))
                    return
                end
                
                % other info
                for j=1:length(other_info)
                    eval(['bfr_other = {data(:).' other_info{j} '};'])
                    bfr_unique = unique(bfr_other(~cellfun('isempty', bfr_other)));
                    if n==1
                        % pre-allocate
                        eval([other_info{j} '_NUM   = zeros(1000,length(mm_ALL)*length(years));'])
                        eval([other_info{j} '_types = bfr_unique;'])
                    end
                    eval(['bfr_types = ' other_info{j} '_types;'])
                    % count
                    for jj=1:length(bfr_unique)
                        clear I
                        I = find(strcmp(bfr_types,bfr_unique{jj}));
                        if ~isempty(I)
                            eval([other_info{j} '_NUM(I,n) = ' ...
                                other_info{j} '_NUM(I,n) + sum(strcmp(bfr_other,bfr_unique{jj}));'])
                        else
                            % if a new type appears, include it in the list
                            eval([other_info{j} '_NUM(length(' other_info{j} ...
                                '_types)+1,n) = sum(strcmp(bfr_other,bfr_unique{jj}));'])
                            eval([other_info{j} '_types{end+1} = bfr_unique{jj};'])
                            
                        end
                    end
                end
                
                % number of data by 1 degree bin in the month/year of interest
                NUM(:,:,n)                          = histcounts2(x,y,x_edges,y_edges);
                
                % check for nans
                clear msk
                msk = y==ArgovisNaN | pmax_temp==ArgovisNaN;
                y(msk)         = nan;
                pmax_temp(msk) = nan;
                
                NUM_by_LAT_by_max_PRES_temp(:,:,n)  = histcounts2(y,pmax_temp,y_edges,press_edges);
                
            end
        else
            break
        end
    end
    if (years(yy)==end_year && mm_ALL(imm)==end_month)
        break
    end
    disp(datestr(now))
end
disp(datestr(now))

% reduce size of NUM for other info
for j=1:length(other_info)
    eval([other_info{j} '_NUM = ' other_info{j} '_NUM(1:length(' other_info{j} ...
        '_types),:);'])
end

save(['test' datestr(now,'ddmmmyy') '.mat'],...
    '*NUM*','*_types','mm_ALL','years','x_bin','y_bin')

%% %%%% PLOTS
close all
fig_pos  = [0.1        0.1       1420        700];

fntsz = 28;

% tag to save fig
clear month_tag
if sum(mm_ALL==1:12)==length(mm_ALL)
    month_tag = 'all';
end

d_tag_file_out = [num2str(min(years)) '_' num2str(max(years)) '_' month_tag];

%%%%%%%% plot the max pressure available for temperature in different years
figure('color','w','position',fig_pos.*[1 1 1 1]);
years_select = [2004 2005 2007 2009 2011 2013 2015 2017];
[mm_bfr,yrs_bfr] = ndgrid(mm_ALL, years);
yrs_bfr = yrs_bfr(:);
for l = 1:length(years_select)
    subplot(2,4,l)
    set(gca, 'ColorOrder', flipud(parula(length(press_bin))));
    set(gca,'NextPlot','replacechildren')
    
    pcolor(y_bin,press_bin,sum(NUM_by_LAT_by_max_PRES_temp(:,:,yrs_bfr==years_select(l)),3)')
    shading flat
    colorbar('horizontal');colormap([[.85 .85 .85]; flipud(copper(length(press_bin)))]);
    
    axis ij 
    axis tight
    
    if l==1 || l==5
        ylabel('Pressure, dbar')
    end
    
    set(gca,'fontsize',22,'linewidth',2,'fontweight','bold')
    pos_bfr = get(gca,'position');
    
    text(5,2300,['year ' num2str(years_select(l))],'fontsize',22)
    
    %set_ticklabel_to_latitude(gca,'xtick')
    
end

set(gcf,'PaperPositionMode','auto');
print('-dpng',['Nprof_in_maxPressForTemp_and_lat_bins.png'],'-r150')

%% %%%%%% plot map for the number of profiles in each bin

clear d d_x d_y
d       = nansum(NUM,3);
d_x     = x_bin;
d_y     = y_bin;

% Rescale data 1-64
d(d==0) = nan;
d       = log10(d);
mn      = min(d(:));
rng     = max(d(:))-mn;
d       = 1+63*(d-mn)/rng; % Self scale data

L       = [0.01 0.02 0.05 0.1 0.2 0.5 1 2 5 10 20 50 100 200 500 1000 2000 5000];%
L_tag   = {'0.01' '0.02' '0.05' '0.1' '0.2' '0.5' '1' '2' '5' '10' '20' '50' '100' '200' '500' '1000' '2000' '5000'};%];
% Choose appropriate or somehow auto generate colorbar labels
l       = 1+63*(log10(L)-mn)/rng; % Tick mark positions

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% make a quick plot
figure('color','w','position',fig_pos.*[1 1 1 1]);
[Y,X] = meshgrid(d_y,d_x);
scatter(X(:),Y(:),10,d(:),'filled')
colormap('parula')
set(gca,'linewidth',2,'fontsize',24)
axis tight
xlabel('Longitude')
ylabel('Latitude')
hC = colorbar;
set(hC,'Ytick',l,'YTicklabel',L_tag);
set(hC,'position',[.925 .38 .02 .3],'fontsize',24)

set(gcf,'PaperPositionMode','auto');
print('-dpng',['Nprof_in' num2str(dx) 'x' num2str(dy) 'bins_' ...
    d_tag_file_out '_quick.png'],'-r150')

%% %%%%%%%%%% now plot the number of profiles by *
clear t M YR
[YR,M]   = meshgrid(years,mm_ALL);
t        = datenum(YR(:),M(:),15);
xtickall = [t(1:24:end)' datenum(year(t(end)+30),month(t(end)+30),1)];

for icase=1:2
    for j=1:length(other_info)
        close all
        figure('color','w','position',fig_pos.*[1 1 1 1]);
        clear bfr* fntsz_leg
        eval(['bfr_other       = ' other_info{j} '_NUM'';'])
        eval(['bfr_other_types = ' other_info{j} '_types;'])
        
        %%%%%%%%%%%%%%%%%%%%%%%% combine some of the types %%%%%%%%%%%%%%%%
        if icase==2
            clear bfr_other_new bfr_other_types_new
            switch other_info{j}
                case 'PLATFORM_TYPE'
                    for ityp=1:length(bfr_other_types)
                        if contains(upper(bfr_other_types{ityp}),'APEX') && ...
                                ~contains(upper(bfr_other_types{ityp}),'D')
                            bfr_other_types{ityp} = 'APEX';
                        elseif contains(upper(bfr_other_types{ityp}),'APEX') && ...
                                contains(upper(bfr_other_types{ityp}),'D')
                            bfr_other_types{ityp} = 'APEX-D';
                            
                        elseif contains(upper(bfr_other_types{ityp}),'ARVOR') && ...
                                ~contains(upper(bfr_other_types{ityp}),'D')
                            bfr_other_types{ityp} = 'ARVOR';
                        elseif contains(upper(bfr_other_types{ityp}),'ARVOR') && ...
                                contains(upper(bfr_other_types{ityp}),'D')
                            bfr_other_types{ityp} = 'ARVOR-D';
                            
                        elseif contains(upper(bfr_other_types{ityp}),'HM2000')
                            bfr_other_types{ityp} = 'HM2000';
                            
                        elseif contains(upper(bfr_other_types{ityp}),'NAVIS')
                            bfr_other_types{ityp} = 'NAVIS-A';
                            
                        elseif contains(upper(bfr_other_types{ityp}),'NEMO')
                            bfr_other_types{ityp} = 'NEMO';
                            
                        elseif contains(upper(bfr_other_types{ityp}),'NINJA') && ...
                                ~contains(upper(bfr_other_types{ityp}),'D')
                            bfr_other_types{ityp} = 'NINJA';
                        elseif contains(upper(bfr_other_types{ityp}),'NINJA') && ...
                                contains(upper(bfr_other_types{ityp}),'D')
                            bfr_other_types{ityp} = 'NINJA-D';
                            
                        elseif contains(upper(bfr_other_types{ityp}),'NOVA')
                            bfr_other_types{ityp} = 'NOVA';
                            
                        elseif contains(upper(bfr_other_types{ityp}),'PALACE')
                            bfr_other_types{ityp} = 'PALACE';
                            
                        elseif contains(upper(bfr_other_types{ityp}),'PROVOR')
                            bfr_other_types{ityp} = 'PROVOR';
                            
                        elseif contains(upper(bfr_other_types{ityp}),'SOLO') && ...
                                ~contains(upper(bfr_other_types{ityp}),'D') && ...
                                ~contains(upper(bfr_other_types{ityp}),'II') && ...
                                ~contains(upper(bfr_other_types{ityp}),'S2A') && ...
                                ~contains(upper(bfr_other_types{ityp}),'S2X') && ...
                                ~contains(upper(bfr_other_types{ityp}),'ALTO')
                            bfr_other_types{ityp} = 'SOLO';
                        elseif contains(upper(bfr_other_types{ityp}),'SOLO') && ...
                                ~contains(upper(bfr_other_types{ityp}),'D') && ...
                                (contains(upper(bfr_other_types{ityp}),'II') || ...
                                contains(upper(bfr_other_types{ityp}),'S2A') || ...
                                contains(upper(bfr_other_types{ityp}),'S2X') || ...
                                contains(upper(bfr_other_types{ityp}),'ALTO'))
                            bfr_other_types{ityp} = 'SOLO-II';
                        elseif contains(upper(bfr_other_types{ityp}),'SOLO') && ...
                                contains(upper(bfr_other_types{ityp}),'D') && ...
                                ~contains(upper(bfr_other_types{ityp}),'II') && ...
                                ~contains(upper(bfr_other_types{ityp}),'S2A') && ...
                                ~contains(upper(bfr_other_types{ityp}),'S2X') && ...
                                ~contains(upper(bfr_other_types{ityp}),'ALTO')
                            bfr_other_types{ityp} = 'SOLO-D';
                        else
                            bfr_other_types{ityp} = 'OTHER';
                        end
                    end
                    bfr_other_types_new = {'APEX' 'APEX-D' ...
                        'ARVOR' 'ARVOR-D' 'PROVOR' ...
                        'SOLO' 'SOLO-II' 'SOLO-D' 'NAVIS-A' ...
                        'HM2000' 'NEMO' 'NINJA' 'NINJA-D' ...
                        'NOVA' 'PALACE'  'OTHER'};
                case 'POSITIONING_SYSTEM'
                    for ityp=1:length(bfr_other_types)
                        if ~contains(upper(bfr_other_types{ityp}),'ARGOS') && ...
                                ~contains(upper(bfr_other_types{ityp}),'GPS') && ...
                                ~contains(upper(bfr_other_types{ityp}),'IRIDIUM')
                            bfr_other_types{ityp} = 'OTHER';
                        end
                    end
                    bfr_other_types_new = {'ARGOS' 'GPS' 'IRIDIUM' 'OTHER'};
                    
                case 'DATA_MODE'
                    % this code is to sort the types
                    %                     isort = [2 1 3];
                    %                     if length(unique(isort))~=length(bfr_other_types) || ...
                    %                             max(isort)~=length(bfr_other_types)
                    %                         stop
                    %                     end
                    %                     bfr_other_types_new = {bfr_other_types{isort}};
                    bfr_other_types_new = bfr_other_types;
                case 'dac'
                    % this code is to sort the types
                    %                     isort = [1 3 5 7 4 2 10 6 8 9 11];
                    %                     if length(unique(isort))~=length(bfr_other_types) || ...
                    %                             max(isort)~=length(bfr_other_types)
                    %                         stop
                    %                     end
                    %                     bfr_other_types_new = {bfr_other_types{isort}};
                    bfr_other_types_new = bfr_other_types;
            end
            if exist('bfr_other_types_new','var')
                bfr_other_new       = ...
                    nan(size(bfr_other,1),length(bfr_other_types_new));
                for ityp=1:length(bfr_other_types_new)
                    bfr_other_new(:,ityp) = ...
                        sum(bfr_other(:,...
                        strcmp(bfr_other_types,bfr_other_types_new{ityp})),...
                        2);
                end
                clear bfr_other bfr_other_types
                bfr_other       = bfr_other_new;
                bfr_other_types = bfr_other_types_new;
                clear bfr_other_new bfr_other_types_new
            end
        end
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
        if length(bfr_other_types)>20
            fntsz_leg = 8;
            
        else
            fntsz_leg = 20;
        end
        bfr_other(bfr_other==0) = nan;
        
        if length(bfr_other_types)>5
            
            cbr = parula(length(bfr_other_types));

        else
            
            cbr = flipud([[178 102 255]; [102 102 255]; [102 178 255]; [200 200 200]; [130 130 130]]./255);%[255 102 102]
        end
        
        set(0,'DefaultAxesColorOrder',cbr)
        set(gca,'NextPlot','replacechildren')
    
        if size(bfr_other,2)<100
            bar(t,bfr_other,'stacked','edgecolor','none','BarWidth',1.125)
            
            set(gca,'linewidth',2,'fontsize',24,'xtick',xtickall(1:end-1),'xlim',...
                [xtickall(1) xtickall(end)])
            
            datetick('x','yyyy','keepticks','keeplimits')
            
            hl = legend(bfr_other_types,'location','eastoutside');
            set(hl,'interpreter','none','fontsize',fntsz_leg)
            
            set(gca,'position',[0.1197    0.1100    0.7133    0.8150],...
                'linewidth',2,'fontweight','bold')
            
            set(gcf,'PaperPositionMode','auto');
            print('-dpng',['Nprof_' d_tag_file_out '_by_' other_info{j} ...
                '_n' num2str(icase,'%02d') '.png'],'-r150')
        end
        
    end
end