% ti_2006_alex_horsec.m
% Tioga-171, Tioga-177, ... Live counts data: Plot horizontal map 
% Olga Kosnyreva, AOPE/WHOI, 5/5/2006
%
% I N P U T:
% (1) Live count data.
% Filenames: live_count_ti171_06.txt, live_count_ti177_06.txt, ...
% (2) CTD Stations list data.
% Filenames: ti171_ctd_list.mat (or .txt), ti177_ctd_list.mat (or .txt), ...

clear;
flops(0);
%
% Save png files.
%saveOption='yes';
saveOption='no';
%
% Vector of contour levels.
vv=[0 10 20 30 50 100 200 300 400 500 700 1000 2000];% like 'old'
% ---- MENU: Tioga cruise
crNum=f_menu_of_numbers_title('Tioga cruise',[171 177]);
if isempty(crNum), return; end
if crNum==171
   ctdList=1:29;
   cr_name='R/V Tioga TI171';
   cr_date='April (11 - 12), 2006';
   cr_id='ti171';
elseif crNum==177
   ctdList=1:29;
   cr_name='R/V Tioga TI177';
   cr_date='April (25 - 26), 2006';
   cr_id='ti177';
end
%*******************************************************************
%**************         C O N S T A N T S        ******************* 
%*******************************************************************
% Draw coast line.
drawCoast='yes';
%drawCoast='no';
if strcmp(drawCoast,'yes')
   %  load 'structure fields' for MARMAP1 mesh:
   %  mp - cartesian, mm - geographycal coordinate 
   %marm_ll_g2s_moody
   load land_data;
   mm=land_data;
   %load land_data_tioga_2006;
   %mm=land_data_tioga_2006;
   mm.x=-mm.x;
   mm.lnd(:,1)=-mm.lnd(:,1);
end
%
% Land, isobaths - color.
bnd_color=[0.71 0.66 0.43];
isobaths_color=[0.38 0.41 1];
color_lbl=[0 0 0.43];
%
% Option: colored or black and white map.
%optColor='inBlackWhite';
optColor='inColor';
if strcmp(optColor,'inBlackWhite')
   ocean_color=[1 1 1];
else
   ocean_color=[0.8 0.96 1];    % Ocean
end
%
% Print file type.
printFile='ps'; % and PNG
%printFile='jpg';
if strcmp(printFile,'jpg')
   fnt_lbl_isobath=8;
   fnt_ax=9;
   fnt_lbl=9;
   fnt_title=10;
   fnt_title_sub=fnt_title;
   fnt_text=9;
   fnt_alex=8;
   marker_size=4;
   text_title_color=[0.3 0.3 0.26];
   text_min_max_color=[0.3 0.3 0.26];
   text_cb_color=text_min_max_color;
   marker_edge_color=[0.5 0.5 0.5];
   marker_face_color=[1 1 0.9];
   plot_id='o';
   marker_edge_color=[0 0 0];
   marker_face_color=[0 0 0];
   fnt_name_ax='lucida';
   fnt_name_lbl='lucida';
   fnt_name_title='lucida';
   line_width=1.2;
   k_text=0.16;
   k_text_y=0.04;
   land_color=[0.8 0.8 0.66];
   d_yc=0.8; 
   ky_tcb=-0.21;
   d_lat=0.033;
else
   fnt_lbl_isobath=8;
   fnt_ax=7;
   fnt_lbl=9;
   fnt_title=10;
   fnt_alex=8;
   fnt_title_sub=fnt_title;
   fnt_text=8;
   marker_size=2;
   text_title_color=[0 0 0];
   %text_min_max_color=[0.4 0.4 0.26];
   text_min_max_color=[0 0 0];
   text_cb_color=text_min_max_color;
   plot_id='o';
   marker_edge_color=[0 0 0];
   marker_face_color=[0 0 0];
   fnt_name_ax='helvetica';
   fnt_name_lbl='helvetica';
   fnt_name_title='helvetica';
   line_width=0.5;
   k_text=0.18;
   k_text_y=0.05;
   land_color=[0.855 0.783 0.663];
   d_yc=0.8;
   ky_tcb=-0.218;
   d_lat=0.029;
end
fig_pos = [134.00        271.00        393.00        510.00];
% COLORMAP DATA:
load palet_12_section_v3;
cmap=palet_12_section_v3;
nband=size(cmap,1);
% FOR COLORMAP GRAY:
% --- Gray colormap: 
%colormap_name='gray';
% --- Any other colormap: 
colormap_name='';
% Alex file columns:
colSt_a=1;
colAlex=8;
% CTD list file: lat, lon columns#.
colSt=1;
colLat=8;
colLon=9;
% Mapping C O N S T A N T S:
makeContourLabel='no';
makeShadingFlat='yes';
makeContours='no';
% Mapping area C O N S T A N T S
%mapArea='full';
%mapArea='area1';
%mapArea='area2';
%mapArea='forJPG';
%mapArea='spec1';
mapArea='spec1_Mass';
%mapArea='spec1_four_sources';
if strcmp(mapArea,'full')
   latmin = 41;
   latmax = 45.5;
   lonmax = 71;
   lonmin = 65;
elseif strcmp(mapArea,'area1')
   latmin = 43.1;
   latmax = 44.1;
   lonmin = 68.3;
   lonmax = 70;
elseif strcmp(mapArea,'area2')
   latmin = 43.1;
   latmax = 44.2;
   lonmin = 67.9;
   lonmax = 70.25;
elseif strcmp(mapArea,'forJPG')
   latmin = 42;
   latmax = 45.2;
   lonmax = 71;
   lonmin = 65.8-0.2;
elseif strcmp(mapArea,'spec1')
   latmin = 41.45;
   latmax = 45.2;
   lonmax = 71+0.1;
   lonmin = 65.8-0.2;
elseif strcmp(mapArea,'spec1_Mass')
   latmin = 41.45;
   latmax = 42.8;
   lonmax = 71+0.1;
   lonmin = 69.7;
elseif strcmp(mapArea,'spec1_four_sources')
   latmin = 41.45-0.2;
   latmax = 42.8;
   lonmax = 71+0.1;
   lonmin = 69.7-0.5;
else
   latmin = 43.1;
   latmax = 44.2;
   lonmin = 67.9;
   lonmax = 70;
end
% Patch for coloring right lower corner appearing as Land. 
patch_x=[66.5 lonmin  lonmin 66.5];
patch_y=[latmin 42.5 latmin latmin];
% Meshgrid limits.
meshLonMin=lonmin;
meshLonMax=lonmax;
meshLatMin=latmin;
meshLatMax=latmax;
% Meshgrid interval for mapping [deg].
%meshInt=0.025;
meshInt=0.01;
% 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;
%**************************************************************************
%****************************     LOAD ALEX DATA     **********************
%**************************************************************************
p=['load live_count_ti' num2str(crNum) '_06.txt;'];
eval(p);
p=['d_a=live_count_ti' num2str(crNum) '_06;'];
eval(p);
%**************************************************************************
%***********************   LOAD CTD DATA LIST   ***************************
%**************************************************************************
p=['load ti' num2str(crNum) '_ctd_list;'];
eval(p);
p=['c_list=ti' num2str(crNum) '_ctd_list;'];
eval(p);
%**************************************************************************
%*******************************    QUANTITY NAMES     ********************
%**************************************************************************
jN='alex_live';
j_title='Surface Live Counts';
j_label='Alexandrium (cells/L)';
%*************************
% Set figure properties;
%*************************
figure;
set(gcf,'color',[1 1 1],'position',fig_pos,...
   'paperPositionMode','auto');
colormap(cmap);
%**************************************************************************
%*************************     LOOP# 1:    STATION#    ********************
%**************************************************************************
DATAC=[];
countNotUsed=0;
count_bad=0;
str_disp='';
for stNum=ctdList
   % Current station Alexandrium data. 
   rst_a=find(d_a(:,colSt_a)==stNum);
   if ~isempty(rst_a)
      parSt=d_a(rst_a,colAlex);
      if length (parSt)>1;
          disp(['Station#' num2str(stNum) ' Alex data: total ' num2str(length(parSt)) ' points']);
          return;
          % Find mean of SHA values.
          %disp(['Station#' num2str(stNum) ' SHA data (total ' num2str(length(parSt)) ' points) is averaged']);
          %parSt=mean(parSt);
      end
      if isnan(parSt)
         % Make a string matrix  - message about NaNs, remove NaN data.
         count_bad=count_bad+1;
         str_disp=strvcat(str_disp,['Station No ' num2str(stNum) ', depth ' num2str(depth) 'm: data not available']);
      else
         % Current station coordinates. 
         rst=find(c_list(:,colSt)==stNum);
         latSt = c_list(rst,colLat); 
         lonSt = c_list(rst,colLon);
         % Make a matrix with the columns: station#, Alex, Lat, Lon data.
         DATAC = [DATAC;stNum parSt latSt lonSt];
      end
   end % For if ~isempty(rst_a)
   clear parSt latSt lonSt;
end % For for stNum=ctdList
%**************************************************************************
%****************************      PLOT MAPS         **********************
%**************************************************************************
%********************
% Lon, lat, Alex  - all stations.
%********************
if size(DATAC,1)<=2, return; end
   LON = DATAC(:,4);
   LAT = DATAC(:,3);
   PAR = DATAC(:,2);
   ST = DATAC(:,1);
%********************
% Meshgrid and griddata. 
%********************
   [LONi,LATi] = meshgrid(meshLonMin:meshInt:meshLonMax,meshLatMin:meshInt:meshLatMax);
   PARi = griddata(LON,LAT,PAR,LONi,LATi,'linear');
%********************
% Qantity Min, Max. 
%********************
   parMin = min(PAR(:));
   [parMax,line_num] = max(PAR(:));
   parMax_st=ST(line_num);
%********************
% Qantity Min, Max after Griddata.
%**************************
   parMin_i = min(PARi(find(~isnan(PARi(:)))));
   parMax_i = max(PARi(find(~isnan(PARi(:)))));
   if parMin_i==parMax_i, return; end    
%*************************
% Calculate scaled Z matrix.
%*************************
   [Z_new,cval] = f_oc412_spec_scale(PARi,vv); 
%*************************
% Plot filled contours.
%*************************
   contourf(LONi,LATi,Z_new,cval);
   caxis([min(cval) max(cval)]);
   set(gca,'xdir','reverse');
%*************************
% Set shading interpolation.
%*************************
   if strcmp(makeShadingFlat,'yes'), shading flat; end
%*************************
% Set axes.     
%*************************
   axPos=get(gca,'position');  
   d_x=-0.03;
   d_y=-0.11;
   k=1.1;
   set(gca,'position',[axPos(1)+axPos(3)*(d_x) axPos(2)+axPos(4)*d_y axPos(3)*k axPos(4)*k]);
%************************* 
% Colorbar.
%*************************
   h1 = colorbar('horiz');
   % Put scaled values to  colorbar ticks and primary data values to colorbar ticks labels.
   set(h1,'xtick',cval,'xticklabel',vv,'xLim',[min(cval) max(cval)],...
      'Fontsize',fnt_ax,'ticklength',get(h1,'ticklength')*1.9,...
      'xcolor',text_cb_color,'ycolor',text_cb_color,'fontname',fnt_name_ax);
   % Put every other ticklabel if the number of bands >10.
   if nband>20
      tlb=f_ticklabel(vv,2);
      set(h1,'xtickLabel',tlb); 
   else
      set(h1,'xticklabel',vv);
   end
   aPos = get(h1,'Position');
   d_xc=-0.09;
   kc=0.85;
   d_xc=(1-kc)/2;
   set(h1,'Position',[aPos(1,1)+aPos(1,3)*d_xc aPos(1,2)+aPos(1,4)*d_yc aPos(1,3)*kc aPos(1,4)*0.3]);
   % Text - quantity name.
   text(lonmax-range([lonmax lonmin])*0.5,latmin+range([latmax latmin])*ky_tcb,j_label,...
      'FontSize',fnt_lbl,....
      'horizontalAlignment','center','color',text_cb_color,'fontname',fnt_name_lbl);
%******************************
% Set axes labels and properties.
%******************************
   axis([lonmin lonmax latmin latmax]);
   set(gca,'FontSize',fnt_ax,...
      'box','on','PlotBoxAspectRatio',[1 relyx*1 1],...
      'ticklength',get(gca,'ticklength')*0.6,'tickdir','out','color',ocean_color,'fontname',fnt_name_ax);
   xlabel('Longitude (degW)','FontSize',fnt_lbl,'fontname',fnt_name_lbl);
   ylabel('Latitude (degN)','FontSize',fnt_lbl,'fontname',fnt_name_lbl);
%************************
% Subplot title.
%************************
   str_title(1)={[j_title]};
   str_title(2)={[cr_name ', ' cr_date]};
   title(str_title,...
     'FontSize',fnt_title,'verticalAlignment','bottom',...
     'color',text_title_color,'fontname',fnt_name_title);
%*********************************************************************************
%**************************       PLOT  LAND          ****************************
%*********************************************************************************
   if strcmp(drawCoast,'yes')
      % Draw land fill.
      hold on;
      lndfill(mm,land_color);
      hold off;
   end % For if strcmp(drawCoast,'yes')
%*************************
% Plot station locations.
%*************************
   hold on;
   plot(LON,LAT,plot_id,'markersize',marker_size,...
      'markerEdgeColor',marker_edge_color,'markerFaceColor',marker_face_color,...
      'linewidth',line_width);
   hold off;
%***********************************************************
%*********************    PLOT  Patch     ******************
%***********************************************************
   % Plot the Patch for coloring right lower corner appearing as Land. 
   if strcmp(mapArea,'spec1_Mass')
   else
      hold on;
      patch(patch_x,patch_y,ocean_color,'faceColor',ocean_color,'edgeColor','none');
      hold off;
   end
%******************************
% Text: min, max.
%****************************** 
   str_txt(1)={['Min=' num2str(parMin)]};
   str_txt(2)={['Max=' num2str(parMax)]};
   %str_txt(2)={['Max=' num2str(parMax) ', station #' num2str(parMax_st)]};
   text(lonmin+range([lonmin lonmax])*k_text,latmin+range([latmin latmax])*k_text_y,str_txt,...
      'FontSize',fnt_text,'horizontalAlignment','left','color',text_min_max_color);
%*****************
% Plot text - alex numbers.
%*****************
for stNum=ctdList
   rst_loc=find(c_list(:,colSt)==stNum);
   lat = c_list(rst_loc,colLat);
   lon = c_list(rst_loc,colLon);
   rst_a=find(d_a(:,colSt_a)==stNum);
   alex = d_a(rst_a,colAlex);
   if crNum==177 & stNum==16
      text(lon,lat+d_lat,[num2str(alex)],...
         'horizontalAlignment','center','verticalAlignment','middle',...
         'fontsize',fnt_alex+1,'color','r');
   else
      text(lon,lat+d_lat,[num2str(alex)],...
         'horizontalAlignment','center','verticalAlignment','middle',...
        'fontsize',fnt_alex);
   end
end 
%**************************
% Plot axes box 
% (if it gets invisible).
%**************************
   hold on;
   plot([lonmax lonmax],[latmin latmax],'-k',...
      [lonmin lonmin],[latmin latmax],'-k',...
      [lonmax lonmin],[latmin latmin],'-k',...
      [lonmax lonmin],[latmax latmax],'-k');
   hold off;
%*******************   
% Get print file - png
%*******************
if strcmp(saveOption,'yes')
   pause(1);
   p=['print -r100 ' cr_id '_alex_liveCounts.png -dpng'];
   eval(p);
end % For if strcmp(saveOption,'yes')
%*******************   
% Display the number of stations where Alex not counted.
%*******************   
%if count_bad>0,...
%   disp(strvcat(blanks(1),[' Alexandrium ' j_title ...
%       ' , total ' num2str(count_bad) ' stations where Alexandrium not counted:'],str_disp,blanks(1))); ...
%end
%***********************************************************************************
%****************************      GRAY COLORMAP         ***************************
%***********************************************************************************
if strcmp(colormap_name,'gray')
   colormap gray;
   gmap=colormap;       % Max values are white, Min are black
   gmap=flipud(gmap);   % FLIPUD: Max values are black now, Min are white
   colormap(gmap);
end
