% nbp1201_read_plot_jgofs.m
%
%@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
%@@@@@@@@@@@@@@@@@   INPUT  @@@@@@@@@@@@@@@@@@@@
%@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
% JGOFS format data files
% http://science.whoi.edu/users/mcgillic/nbp1201/process/1201JGOF.tar
%@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
% PLOT data locations and values in colormap,
% one day data from days [361:365 1:36].

clear;
t0 = clock;

% Path, names.
path_part_1='/data2/antarctic_data2/nbp1201/underway/process_ship/';
path_part_2='JGOF/';
path_name=[path_part_1 path_part_2];
fname_part=['jg']; % e.g. jg001.dat
path_name_out=[path_name 'process/'];
fname_ext='dat';
npar=26;
%******************
% C O N S T A N T S
%******************
% Font size.
fnt_axis=10;
fnt_lbl=10;
fnt_title=10;
fnt_text=8;
fnt_name='helvetica';
marker_size=2;
marker_size_time=9;
fig_pos=[360.00        460.00        599.00        461.00];
crName='NBP1201';
%@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
%@@@@@@@@@@   Data file columns:  @@@@@@@@@@@
%@@@@@@@@@@@ Date/Time/Location  @@@@@@@@@@@@
%@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
colDay_M=1;
colMonth_M=2;
colYear_M=3;
colHour_M=4;
colMin_M=5;
colSec_M=6;
colLat_M=7; % NGL Latitude tt.tttt (negative is South)
colLon_M=8; % NGL Longitude ggg.gggg (negative is West)
%@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
%@@@@@@@@@@@@@  Day/Date list   @@@@@@@@@@@@@
%@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
day_list_string=strvcat('361: Dec 27,2011',...
    '362: Dec 28, 2011',...
    '363: Dec 29, 2011',...
    '364: Dec 30, 2011',...
    '365: Dec 31, 2011',...
    '001: Jan 1, 2012',...
    '002: Jan 2, 2012',...
    '003: Jan 3, 2012',...
    '004: Jan 4, 2012',...
    '005: Jan 5, 2012',...
    '006: Jan 6, 2012',...
    '007: Jan 7, 2012',...
    '008: Jan 8, 2012',...
    '009: Jan 9, 2012',...
    '010: Jan 10, 2012',...
    '011: Jan 11, 2012',...
    '012: Jan 12, 2012',...
    '013: Jan 13, 2012',...
    '014: Jan 14, 2012',...
    '015: Jan 15, 2012',...
    '016: Jan 16, 2012',...
    '017: Jan 17, 2012',...
    '018: Jan 18, 2012',...
    '019: Jan 19, 2012',...
    '020: Jan 20, 2012',...
    '021: Jan 21, 2012',...
    '022: Jan 22, 2012',...
    '023: Jan 23, 2012',...
    '024: Jan 24, 2012',...
    '025: Jan 25, 2012',...
    '026: Jan 26, 2012',...
    '027: Jan 27, 2012',...
    '028: Jan 28, 2012',...
    '029: Jan 29, 2012',...
    '030: Jan 30, 2012',...
    '031: Jan 31, 2012',...
    '032: Feb 1, 2012',...
    '033: Feb 2, 2012',...
    '034: Feb 3, 2012',...
    '035: Feb 4, 2012',...
    '036: Feb 5, 2012');
%@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
%@@@@@@@@@@@@@  Quantites list   @@@@@@@@@@@@
%@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
par_list=strvcat('09 Speed over ground, Knots',...
    '10 GPS HDOP',...
    '11 Gyro Heading, Degrees (azimuth)',...
    '12 Course made good, Degrees (azimuth)',...
    '13 Mast PAR, MkEinsteins/meter^2 sec',...
    '14 Sea surface temperature, degC',...
    '15 Sea surface conductivity, siemens/meter',...
    '16 Sea surface salinity, PSU',...
    '17 Sea depth (uncorrected, calc. sw sound vel. 1500 m/s), meters',...
    '18 True wind speed (max speed windbird), meters/sec',...
    '19 True wind direction (max speed windbird), degrees (azimuth)',...
    '20 Ambient air temperature, degC',...
    '21 Relative humidity, %',...
    '22 Barometric pressure, mBars',...
    '23 Sea surface fluorometry, Mkg/l (mg/m^3)',...
    '24 Transmissometry, %',...
    '25 PSP, W/m^2',...
    '26 PIR, W/m^2');
%@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
%@@@@@@@@@@@@@@@  MENUS   @@@@@@@@@@@@@@@@@@@
%@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
%------------------------------------------------
%------ MENU (1): Day of the year: Date  --------
%------------------------------------------------
day_name=f_menu_of_strings('The day of year: Date',day_list_string);
if isempty(day_name), return; end
%-----------------------------------------------
day=str2num(day_name(1:3));
day_date=deblank(day_name(6:end));
date_vec=datevec(day_date);
year=date_vec(1,1);    
month=date_vec(1,2);    
date_day=date_vec(1,3);    
if day<10
    day_str=['00' num2str(day)];
elseif day>=10 & day<100
    day_str=['0' num2str(day)];
else
    day_str=['' num2str(day)];
end
%------------------------------------------------
%-------------------- MENU (2): QUANTITY --------
%------------------------------------------------
colNum_name=f_menu_of_strings('QUANTITY Column: Name',par_list);
if isempty(colNum_name), return; end
%-----------------------------------------------
colNum=str2num(colNum_name(1:2));
q_name=deblank(colNum_name(4:end));
%*******************************************************************
%******************    Open file for reading   *********************
%*******************************************************************
fname=[fname_part day_str];
pfid = fopen([path_name fname '.' fname_ext]);
if pfid<0
    disp(['File not found: ' fname '.' fname_ext ]);
    return;
end
count_line=0;
data=[];
while 1
    count_line=count_line+1;
    line = fgetl(pfid);
    if ~ischar(line), break, end
    %if ~isstr(line), break, end
    line_cor1 = strrep(line,':',char(9));
    line_cor2 = strrep(line_cor1,'/',char(9));
    line_cor2 = strrep(line_cor2,['NA' char(9)],['NAN' char(9)]);
    % Read line under format control.
    line_data = (sscanf(line_cor2,'%f'))';
    % CHECK vectors length. If 'good' - MAKE MATRIX, if 'bad' display a message:
    if length(line_data)==npar
        data = [data;line_data];
    else
        disp(['Character (not number) is found in a line ' num2str(count_line) ':']);
        disp(line);
        disp('This part of line is read sucessfully:');
        disp(line_data);
        bad_line=count_line
        return;
    end
    clear line line_data line_cor1 line_cor2;
end  % For while 1
% Close file.
fclose(pfid);
d_uw=data;
clear data;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%  Make West Lon positive %%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Make W Lon positive - all E lon
r_w=find(d_uw(:,colLon_M)<0);
d_uw(r_w,colLon_M)=d_uw(r_w,colLon_M)+360;
clear r_w;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%  E X C L U D E  NaN data  %%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
r_nan=find(isnan(d_uw(:,colNum)));
d_uw(r_nan,:)=[];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%  DATA %%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if isempty(d_uw)
    disp(['Day #' num2str(day) ': ' channel_id  ' ' q_name ' ' subname_str ': no data'])
    return;
else % For if isempty(d_uw)
    % Decimal hour.
    h_dec_N=d_uw(:,colHour_M) + ( d_uw(:,colMin_M) + d_uw(:,colSec_M)/60 )/60;
    % Quantity and location.
    q_uw=d_uw(:,colNum);
    lat_uw=d_uw(:,colLat_M);
    lon_uw=d_uw(:,colLon_M);
    % Data min, max.
    data_min_max=[min(q_uw) max(q_uw)];
    % Year, Month.
    year=d_uw(1,colYear_M)+2000;
    month=d_uw(1,colMonth_M);
    %#############################################
    %############   Map limits  ##################
    %#############################################
    lonmin_d = min(lon_uw);
    lonmax_d = max(lon_uw);
    latmin_d = min(lat_uw);
    latmax_d = max(lat_uw);
    klat=0.05;
    klon=0.05;
    if day==361 | day==362 | day==3 | day==4 | day==6
        klat=0.15;
    elseif day==19  | day==25  | day==35
        klat=1.3;
    elseif day==33
        klon=1.3;
    end
    d_lon=range([lonmin_d lonmax_d])*klon;
    d_lat=range([latmin_d latmax_d])*klat;
    lonmin=lonmin_d-d_lon;
    lonmax=lonmax_d+d_lon;
    latmin=latmin_d-d_lat;
    latmax=latmax_d+d_lat;
    lon_int=1;
    lat_int=0.5;
    % --------  Calculate relation relyx = axisy/axisx (for PlotBoxAspectRatio)
    weg = 111*cos(pi*(latmin+(latmax-latmin)/2)/180);% Weight of one degree longitude
    xc = (lonmax-lonmin)*weg;                        % Longitude, km
    yc = (latmax-latmin)*111;                        % Latitude, km
    relyx = yc/xc;
    %#############################################
    %######## X,Y ticks, tickLabels ##############
    %#############################################
    % ---- X-ticks.
    x_tick=floor(lonmin):lon_int:lonmax;
    x_tick_label='';
    for ii=1:length(x_tick)
        if x_tick(ii)<180
            x_tick_str=[num2str(x_tick(ii)) 'E'];
        else
            x_tick_str=[num2str(360-x_tick(ii)) 'W'];
        end

        if ii==length(x_tick)
            delim='';
        else
            delim='|';
        end
        x_tick_label=[x_tick_label x_tick_str delim];
    end
    clear ii;
    % ---- Y-ticks.
    y_tick=floor(latmin):lat_int:latmax;
    y_tick_label='';
    for ii=1:length(y_tick)
        if ii==length(y_tick)
            delim='';
        else
            delim='|';
        end
        y_tick_label=[y_tick_label num2str(-y_tick(ii)) 'S' delim];
    end
    clear ii;
    %#############################################
    %############  COLORMAP data  ################
    %#############################################
    % Load Colormap.
    load palet_13_ant;
    cmap=palet_13_ant;
    n_bands=size(cmap,1);
    n_dec_chars=2;
    cval=1:size(cmap,1)+1;
    % Quantity min, max.
    q_min_d = data_min_max(1,1);
    q_max_d = data_min_max(1,2);
    %***********************************************************
    %****************  Set figure properties  ******************
    %***********************************************************
    figure;
    set(gcf,'position',fig_pos,'color',[1 1 1],...
        'paperPositionMode','auto');
    %***********************************************************
    %*****************  PLOT TIME Markers    *******************
    %***********************************************************
    d_hour=[0 4 8 12 16 20 23];
    %date_format='yyyy-mm-dd HH:MM:SS';
    %date_format='mm/dd/yyyy HH:MM:SS';
    %date_format='mm/dd';
    %date_format='yyyy/mm/dd HH:MM';
    date_format='HH:MM';
    hold on;
    for nn=1:length(d_hour)
        dnum_n=datenum([year month date_day d_hour(nn) 0 0]);
        r_h=find(d_uw(:,colHour_M)==d_hour(nn) & d_uw(:,colMin_M)==0);
        if ~isempty(r_h)
            % use the point with min HOUR time.
            [mi,rmi]=min(h_dec_N(r_h,1));
            lat=lat_uw(r_h(rmi),1);
            lon=lon_uw(r_h(rmi),1);
            text(lon,lat-range([latmin latmax])*0.06,['  ' datestr(dnum_n,date_format)],...
                'fontsize',fnt_text,'fontname',fnt_name,'horizontalAlignment','left','color',[1 0 0]);
            plot(lon,lat,'k+','markersize',marker_size_time,'markerEdgeColor',[1 0 0]);
        end
    end
    hold off;
    if q_min_d==q_max_d
        %***********************************************************
        %*******************       MIN=MAX      ********************
        %***************   PLOT locations only  ********************
        %***********************************************************
        hold on;
        plot(lon_uw,lat_uw,'o','color',[0 0 1],'markersize',marker_size);
        hold off;
    else % For if q_min_d==q_max_d
        % Calculate colormap vector.
        vv=f_colormap_vector(q_min_d,q_max_d,n_bands,n_dec_chars);
        if isempty(vv)
            disp(strvcat('Data range is too small to calculate the colormap vector',...
                'Try to increase the number of decimal points (function f_colormap_vector input variable n_dec_chars'));
           return;
        end
        colormap(cmap);
        %@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
        %@@@@@@  PLOT markers in COLORMAP color @@@@@@@@@@@@
        %@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
        hold on;
        D=[];
        for jj=1:size(q_uw,1)
            q_j=q_uw(jj,1);
            lat=lat_uw(jj,1);
            lon=lon_uw(jj,1);
            countSearch=0;
            for ii=1:length(vv)-1
                countSearch=countSearch+1;
                if q_j<=vv(1)
                    color_j=cmap(1,:);
                elseif q_j>=vv(end)
                    color_j=cmap(end,:);
                elseif q_j>=vv(ii) & q_j<vv(ii+1)
                    color_j=cmap(ii,:);
                else
                    color_j=[];
                end
                if ~isempty(color_j), break; end
            end
            if ~isempty(color_j)
                % Check.
                D=[D;q_j ii cmap(ii,:)];
                % Plot marker in colormap color.
                plot(lon,lat,'o','markersize',marker_size,...
                    'markeredgecolor',color_j,'markerfacecolor',color_j,...
                    'linewidth',0.3);
            else % For if ~isempty(color_j)
                disp(['Color not found for the point: #' num2str(jj) ', value=' num2str(q_j)]);
                return;
            end % For  if ~isempty(color_j)
        end
        hold off;
        %***********************************************************
        %***********************   COLORBAR   **********************
        %***********************************************************
        cb_ticklabel=vv;
        caxis([min(cval) max(cval)]);
        h1 = colorbar;
        aPos = get(h1,'Position');
        set(h1,'Position',[aPos(1,1)+aPos(1,3)*1.25 aPos(1,2)+aPos(1,4)*0.05 aPos(1,3)*0.5 aPos(1,4)*0.9]);
        set(h1,'ytick',cval,'yticklabel',cb_ticklabel,'yLim',[min(cval) max(cval)],...
            'Fontsize',fnt_axis,'ticklength',get(h1,'ticklength')*2.6,...
            'xcolor',[0 0 0],'ycolor',[0 0 0],'fontname',fnt_name);
    end % For if q_min_d==q_max_d
    %***********************************************************
    %*****************  Set Axes Properties  *******************
    %***********************************************************
    axPos = get(gca,'Position');
    kk=1.0;
    set(gca,'position',[axPos(1)+axPos(3)*(-0.04) axPos(2)+axPos(4)*0.03 axPos(3)*kk axPos(4)*kk]);
    axis([lonmin lonmax latmin latmax]);
    set(gca,'PlotBoxAspectRatio',[1 relyx*1 1],...
        'tickdir','out',...
        'ticklength',(get(gca,'ticklength')*0.6),...
        'xtick',x_tick,'xticklabel',x_tick_label,...
        'ytick',y_tick,'yticklabel',y_tick_label,...
        'fontname',fnt_name,'Fontsize',fnt_axis,...
        'box','on');
    xlabel('Longitude','fontname',fnt_name,...
        'Fontsize',fnt_lbl);
    ylabel('Latitude','fontname',fnt_name,...
        'Fontsize',fnt_lbl);
    title([crName ' JGOFS ' day_date ' (day #' num2str(day) '):  ' q_name],...
        'fontname',fnt_name,...
        'Fontsize',fnt_title,'color',[0 0 0]);
    %***********************************************************
    %********************   PLOT Min/Max   *********************
    %***********************************************************
    str_mi_ma(1)={['Min = ' num2str(data_min_max(1,1))]};
    str_mi_ma(2)={['Max = ' num2str(data_min_max(1,2))]};
    text(lonmin+range([lonmin lonmax])*0.1,latmin+range([latmin latmax])*0.85,str_mi_ma,...
        'fontsize',fnt_text,'fontname',fnt_name,'horizontalAlignment','left');
    %***********************************************************
    %*************  Find measurements interval *****************
    %***********************************************************
    % --- choose one hour which have max data points to compute measurements interval.
    h_n=[];
    for h_i=[0:23]
        rr=find(d_uw(:,colHour_M)==h_i);
        if ~isempty(rr), h_n=[h_n;h_i length(rr)]; end
    end
    if min(h_n(:,2))==max(h_n(:,2)) & min(h_n(:,2))==1 & sum(h_n(:,2))==24
        time_int_str=[ ', Time interval = 1hour'];
    elseif min(h_n(:,2))==max(h_n(:,2)) & sum(h_n(:,2))<24
        time_int_str='';
    else
        [ma,ma_r]=max(h_n(:,2));
        hour_use_for_interval=h_n(ma_r,1);
        rrr=find(d_uw(:,colHour_M)==hour_use_for_interval);
        % --- time interval in seconds.
        time_dif=diff(h_dec_N(rrr,1));
        % --- exclude max to cosider missing data.
        [mad,madr]=max(time_dif);
        time_dif(madr)=[];
        % --- average time interval.
        time_int_sec=mean(time_dif)*3600;
        if time_int_sec>=60
            % --- time interval string - in minutes.
            time_int_str=[ ', Time interval = ' num2str(time_int_sec/60,'%2.2f') ' minute(s)']; 
        else
            % --- time interval string - in seconds.
            time_int_str=[ ', Time interval = ' num2str(time_int_sec,'%2.2f') ' second(s)'];
        end
    end % For if min(h_n(:,2))==max(h_n(:,2)) & sum(h_n(:,2))==24
    %***********************************************************
    %*************   PLOT Text - file info   *******************
    %***********************************************************
    % --- control text position depending on plot box aspect ratio.
    kx_text=0.05;
    if relyx>=2.4
        ky_text=0.14;
        kx_text=0.7;
    elseif relyx>=1.69
        ky_text=0.14;
        kx_text=0.4;
    elseif relyx>=0.87
        ky_text=0.14;
        kx_text=0.1;
    elseif relyx>=0.82
        ky_text=0.14;
    elseif relyx>=0.67
        ky_text=0.16;
    elseif relyx<=0.195
        ky_text=0.6;
    elseif relyx<=0.32
        ky_text=0.4;
    else
        ky_text=0.3;
    end
    % --- plot text.
    str_t(1)={['File name: ' fname '.' fname_ext ' (column #'  num2str(colNum) ')' time_int_str]};
    text(lonmin-range([lonmin lonmax])*kx_text,latmin-range([latmin latmax])*ky_text,str_t,...
        'fontsize',fnt_text,'fontname',fnt_name,...
        'horizontalAlignment','left','fontAngle','it','color',[0 0 1]);
end % For if isempty(d_uw)
estimated_time_minutes=etime(clock,t0)/60
