plot_map.m

function plot_map(hTsgGUI, hPlotAxes)
% Function to plot the earth map and ship trackline
%
% Input
% -----
% hTsgGUI ............ Handel to the main user interface
% hPlotAxes .......... Handels to the graphic axes
%
% Output
% ------
%
% Library : 'M_MAP'
%
% TODOS
%
% note that a decimal degree notation is used, so that a longitude of
% 120?30'W is specified as -120.5
%
% PROBLEM NOT YET RESOLVED
%    Does M_MAP use 0-360? or -180? +180? covention ?
%

global DEBUGGING

% Retrieve named application data
% -------------------------------
tsg = getappdata( hTsgGUI, 'tsg_data');

% save pointer and change it to watch
% -----------------------------------
pointer = get(hTsgGUI, 'pointer');
set(hTsgGUI, 'pointer', 'watch');
set(findobj('tag','MAP_FIGURE'), 'pointer', 'watch');

% Retrieve map border
% --------------------
border = tsg.preference.map.border;

% Retrieve map resolution
% set 'checked' property to 'on' for current map resolution if user
% change it from menu option/preferences
% ------------------------------------------------------------
hdl_map_resolution = flipud(findobj( '-regexp', 'tag', 'TAG_UIMENU_MAP_RESOLUTION_'));
set(hdl_map_resolution, 'Checked', 'off');  % set all menu off
set(hdl_map_resolution(tsg.preference.map.resolution), 'checked', 'on');

% Get the Geographic limit of the TSG time series
% -----------------------------------------------
dateLim = get(hPlotAxes(1), 'Xlim');
ind = find( tsg.DAYD >= dateLim(1) & tsg.DAYD <= dateLim(2));

% m_grid need the use of double instead single
% --------------------------------------------
latx = double(tsg.LATX);
lonx = double(tsg.LONX);

if ~isempty( ind )
  
  latMin = min(latx(ind) );
  if latMin < -90
    latMin = -90;
  end
  latMax = max( latx(ind) );
  if latMax > 90
    latMax = 90;
  end
  lonMin = min( lonx(ind) );
  lonMax = max( lonx(ind) );
  
  % Oversize window due to the large frame
  %--------------------------------------
  latRange = (latMax-latMin);
  latMin = max(floor(latMin), -90);
  latMax = min(ceil(latMax), 90);
  
  lonRange = (lonMax-lonMin);
  lonMin = floor(lonMin);
  lonMax = ceil(lonMax);
  
  lonRange = (lonMax-lonMin);
  if lonRange>=360
    lonMin = -183.6;     %to account for fancy frame
    lonMax = 183.6;
    tsg.lonplus = 180;
    tsg.lonmod = 360;
  else
    tsg.lonplus = 0;
    tsg.lonmod = 0;
  end
  
  % Positionning the right axes (set map current axe)
  % -------------------------------------------------
  axes(hPlotAxes(4));
  
  % Use of Mercator projection
  % --------------------------
  m_proj('Mercator','lat',[latMin-border latMax+border],'long',[lonMin-border lonMax+border]);
  
  % plot climatology mean contour if available
  % ------------------------------------------
  hdl_pushtool_climato = findobj( '-regexp', 'tag', 'PUSHTOOL_CLIM');
  
  % plot climato on map only when the main pushtool button climato is
  % selected
  % -----------------------------------------------------------------
  if strcmp(get(hdl_pushtool_climato, 'state'), 'on')
    
    % display debug infos
    if DEBUGGING
      fprintf(1, 'Plot climato on map:\n');
      tsg.levitus.data
    end
    
    % fill climato uicontrol
    set(findobj('Tag', 'CLIMATO_MIN_VALUE'), 'string', ...
      num2str(tsg.preference.map.climato.(tsg.plot.sample).min));
    set(findobj('Tag', 'CLIMATO_MAX_VALUE'), 'string', ...
      num2str(tsg.preference.map.climato.(tsg.plot.sample).max));
    set(findobj('Tag', 'CLIMATO_STEP_VALUE'), 'string', ...
      num2str(tsg.preference.map.climato.(tsg.plot.sample).step));
    
    climato_contour_scale = (tsg.preference.map.climato.(tsg.plot.sample).min : ...
      tsg.preference.map.climato.(tsg.plot.sample).step : ...
      tsg.preference.map.climato.(tsg.plot.sample).max);
    
    % get climato array for the proper parameter
    climato = tsg.levitus.data.(['WOA_MEAN_' tsg.plot.sample]);
    % get last selected climatology structure
    s = get(findobj('Tag', 'TAG_UIMENU_CLIMATO_MAIN'), 'Userdata');
    level = tsg.preference.levitus_depth_value;
    % extract the matrix
    climato = squeeze(climato(s.time, level,:,:));
    % use local var for quick access
    LATX = tsg.levitus.data.WOA_LATX;
    LONX = tsg.levitus.data.WOA_LONX;
    
    % swap negative value (west) to 180-360
    if ~isempty(tsg.indcross)
      LONX(LONX < 0) = LONX(LONX < 0) + 360;
      [val, indval] = sort(LONX);
      LONX          = LONX(indval);
      climato       = climato(:,indval);
    end
    
    % find  indices for the selected area
    indLat = find(LATX > latMin-border & LATX <latMax+border);
    indLon = find(LONX > lonMin-border & LONX <lonMax+border);
    lon = LONX(indLon);
    lat = LATX(indLat);
    
    % plot 2D climatogogy on map with pcolor
    switch tsg.preference.map.climatology
      case 'pcolor'
        % mask climato range uipanel on map
        set(findobj('tag', 'TAG_MAP_CLIMATO_UIPANEL'),'Visible', 'off');
        % resize map panel for axes
        set(findobj('tag', 'TAG_MAP_CLIMATO_AXES'),'Position',[0, 0, 1, 1])
        % We shift half a grid point
        dx = mean(diff(lon));
        dy = mean(diff(lat));
        % plot climato with pcolor
        m_pcolor(lon-dx/2,lat-dy/2, climato(indLat, indLon));
        shading flat;
        %colormap(m_colmap('jet','step',10));
        colorbar;
        
      case 'contourf'
        % resize map panel for axes
        set(findobj('tag', 'TAG_MAP_CLIMATO_AXES'),'Position',[0, 0, .8, 1])
        % set climato range uipanel on map visible
        set(findobj('tag', 'TAG_MAP_CLIMATO_UIPANEL'),'Visible', 'on');
        % plot climato with contourf
        m_contourf(lon, lat, climato(indLat, indLon), climato_contour_scale);
        colormap(m_colmap('jet','step',10));
        colorbar;
        
      case 'none'
        colorbar('off');
        % mask climato range uipanel on map
        set(findobj('tag', 'TAG_MAP_CLIMATO_UIPANEL'),'Visible', 'off');
        % resize map panel for axes
        set(findobj('tag', 'TAG_MAP_CLIMATO_AXES'),'Position',[0, 0, 1, 1])
    end
  else
    colorbar('off');
  end
  
  % select map type and resolution
  % ------------------------------
  switch tsg.preference.map.resolution
    
    case 1
      % Low-resolution coast lines,
      % --------------------------
      if tsg.preference.map.patch_value == 2
        m_coast('patch',[.7 .7 .7], 'TAG', 'TAG_PLOT4_LINE_COAST');
      else
        s = load('world_c.mat');
      end
      
    case 2
      % Medium-resolution coast lines
      % -----------------------------
      if tsg.preference.map.patch_value  == 2
        m_gshhs_l('patch',[.7 .7 .7], 'TAG', 'TAG_PLOT4_LINE_COAST');
      else
        s = load('world_l.mat');
      end
      
    case 3
      % Intermediate-resolution coast lines
      % ----------------------------
      if tsg.preference.map.patch_value  == 2
        m_gshhs_i('patch',[.7 .7 .7], 'TAG', 'TAG_PLOT4_LINE_COAST');
      else
        s = load('world_i.mat');
      end
      
    case 4
      % High-resolution coast lines
      % ----------------------------
      if tsg.preference.map.patch_value  == 2
        m_gshhs_h('patch',[.7 .7 .7], 'TAG', 'TAG_PLOT4_LINE_COAST');
      else
        s = load('world_h.mat');
      end
      
    otherwise
      tsg.preference.map.resolution = 1;
      % Low-resolution coast lines
      % --------------------------
      if tsg.preference.map.patch_value  == 2
        m_coast('patch',[.7 .7 .7], 'TAG', 'TAG_PLOT4_LINE_COAST');
      else
        s = load('world_c.mat');
      end
      
  end
  
  % Make a grid on the map with fancy box
  % -------------------------------------
  m_grid('box', 'fancy', 'tickdir', 'in', 'TAG', 'TAG_PLOT4_LINE_GRID', ...
    'Fontsize', tsg.fontSize);
  
  % use 'world_x.mat' as tsg.preference.map.patch_string == 'off'
  % --------------------------------------------------------------
  if tsg.preference.map.patch_value == 1
    x = s.ncst(:,1);
    y = s.ncst(:,2);
    % swap negative value (west) to 180-360
    if ~isempty(tsg.indcross)
      x(x < 0) = x(x < 0) + 360;
      % remove line horizontal lines
      x(x == 0) = nan;
    end
    
    % plot coast
    m_line(x,y, 'color','k');
  end
  
  % Plot ship track with QC color
  % -----------------------------
  
  % Get list of keys from hashtable tsg.qc.hash, defined inside
  % tsg_initialisation.m
  % -----------------------------------------------------------
  qc_list = keys(tsg.qc.hash);
  
  % Plot Sample/TSG differences on axe 2
  % iterate (loop) on each key store inside hastable
  % Iterate from the end of the cell array : fliplr
  % ------------------------------------------------
  for key = fliplr( qc_list )
    
    % get  some values in hashtable
    % -----------------------------
    %qcState = tsg.qc.hash.(key).state;
    qcCode  = tsg.qc.hash.(key).code;
    qcColor = tsg.qc.hash.(key).color;
    
    % plot tsg salinity sample with right code/color
    % ----------------------------------------------
    ind = find( tsg.SSPS_QC == qcCode & isnan(tsg.SSPS) == 0);
    if ~isempty( ind )
      
      % Create line on the map
      % ----------------------
      m_line( mod(lonx(ind) + tsg.lonplus, tsg.lonmod) - tsg.lonplus,...
        latx(ind), 'LineStyle', 'none', 'marker','*',...
        'markersize', tsg.preference.map.track.size, 'color', qcColor);
    end
  end
  
  % Write title
  % -----------
  if strcmp(tsg.preference.map.climatology, 'none')  % if climato selected on 2D map
    title(tsg.file.name,  'fontsize',tsg.fontSize,...
      'fontweight', 'bold', 'interpreter','none');
  else
    title([tsg.file.name, ' / ', ...
      tsg.levitus.version, ' ', tsg.levitus.type], 'fontsize',tsg.fontSize,...
      'fontweight', 'bold', 'interpreter','none');
  end
  
  % Set tag for each line
  % ---------------------
  hLines = get(hPlotAxes(4), 'Children');
  for i=1:length(hLines)
    set(hLines(i), 'Tag', ['TAG_PLOT4_LINE_MAP' int2str(i)]);
  end
  
  % Save tsg structure
  % ------------------
  setappdata( hTsgGUI, 'tsg_data', tsg);
  
  % restore pointer
  % ---------------
  set(hTsgGUI, 'pointer', pointer);
  set(findobj('tag','MAP_FIGURE'), 'pointer', 'arrow');
  
end