% Use Argovis API to get metadata globally
clear all
% select region and years of interest
dx      = 1;
dy      = 1;
x_bin   = -179.5:dx:179.5;
y_bin   = -89.5:dy:89.5;

years   = 2004%:2018;
mm_ALL  = 1:12;

% 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;

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

% other info to get
other_info                = {'PLATFORM_TYPE' 'platform_number' 'PI_NAME' ...
    'POSITIONING_SYSTEM' 'POSITIONING_SYSTEM' ...
    'DATA_MODE' 'x_id'};

% still need to deal with following
% lat, lon, position_qc, 'STATION_PARAMETERS_inMongoDB'  'pres_max_for_TEMP' 'pres_min_for_TEMP' 'pres_max_for_PSAL' ...
%   'pres_min_for_PSAL' 'cycle_number'
% allocate the output var
NUM     = nan(length(x_bin),length(y_bin),length(mm_ALL)*length(years));

n =0;
disp(datestr(now))
for yy=1%:length(years)
    for imm=1%:length(mm_ALL)
        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_name ='https://argovis.colorado.edu'; %http://206.189.72.50:3000
        url = [URL_name '/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*
        
        % get data from Argovis
        tic;
        data   = webread(url,options);
        toc;
        
        % store the info of interest
        n = n+1;
        if ~isempty(data)
            % store the location
            % store the location
            for jk = 1:length(data)
                x(jk)  = [data{jk}.lon];
                y(jk)  = [data{jk}.lat];
                
                
                % store the time
                bfr_t  = reshape([data{jk}.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(t_dnum(l)))
                    stop
                end
                
                % other info
                for j=1:length(other_info)
                    eval(['bfr_other = {data{jk}.' other_info{j} '};'])
                    bfr_unique = unique(bfr_other(~cellfun('isempty', bfr_other)));
                    if n==1
                        % pre-allocate
                        eval([other_info{j} '_NUM   = zeros(length(bfr_unique),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
            end
            
            % number of data by 1 degree bin in the month/year of interest
            NUM(:,:,n) = histcounts2(x,y,x_edges,y_edges);
        end
        
    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('info.mat','*NUM*','*_types','mm_ALL','years','x_bin','y_bin')

% plot NUM using x_bin and y_bin as coord