%
% PURPOSE:  Cut a large .grd map up into equally-sized tiles.  Each tile
% is also a .grd.  Then use grdOctreeMaker to convert the .grd tiles to
% .bo's. 
%
% First, overlay Dave's EMARK_vents_MAUVoverEM124_Topo1m_UTM map with grid
% Rob M., 12 July 2018 (Axial)
% Rob M., 28 Mar 2023
%

%
% Path name to map:
% parentMapName = '/Volumes/AUVDevelopment/rob/Maps/20180715AxialTopoForRob/AxialSummit_MAUV_All_Topo1m_UTM.grd';
% parentMapName = '~/Maps/AxialTiles2019/AxialSummit_MAUV_and_Topo1mUTM_2019.grd';
  parentMapName = '~/Maps/AxialTiles2020/AxialSummit_MAUV_All_Topo1m_UTM.grd';
  parentMapName = '~/Maps/AtlanticTiles2023/EMARK_vents_MAUVoverEM124_Topo1m_UTM.grd';
%
% Define the tile boundaries.  The tiles are 1.5kmx1.5km, with center
% separations of 1.5km, which creates non-overlapped regions of
% 1x1km. So there are 13.5x21/1.5^2 126 tiles.
% 1x1km. So there are 2x2.5/1.5^2 = 2.2 tiles??
%
% Start with arrays of the lower left corners.
%
% For ~/Maps/AxialTiles2019/AxialSummit_MAUV_and_Topo1mUTM_2019.grd
  ec = [ 415500:1500: 428500]; % Easting corners
  nc = [5078500:1500:5099000];  % Northing corners

if(~exist('position'))
  position = get(0,'screensize');
  npposition(1:2) = position(1:2) + round([position(3)*.02 position(4)*.04]);
  npposition(3:4) = round([position(3)*.95 position(4)*.90]);
  position = npposition;
end

  lstr={};
  if exist( 'plgFigPtr' ) ~= 1
    plgFigPtr= figure;
    set( plgFigPtr, 'Name', 'Axial Grid', 'DeleteFcn', 'clear plgFigPtr' );
    set(gcf,'position',position)
    if(exist('fontsize')) 
      MbTrnFontsize=fontsize;
    else
      MbTrnFontsize=20;
    end
    MbTrnFontweight='normal';
    set(gca,'fontsize',MbTrnFontsize,'fontweigh',MbTrnFontweight)
    set(gcf,'defaultaxesfontsize',MbTrnFontsize,'defaulttextfontsize', ...
	    MbTrnFontsize)
  else
    figure(plgFigPtr);
  end

% parentMapName  
  if( exist('plg') == 0 )
    if( exist(parentMapName) == 2 )
      [plg.eastings plg.northings plg.z]=grdread2rob(parentMapName);
    else
      fprintf('Map %s doesn''t exist.  Exit.\n', parentMapName );
    end
  end
  
%
% Compute the tile boundaries using ec and nc defined above.  
  llc = [];
  for( i=1:length(nc) )
    llc = [llc;[ec' nc(i)*ones(length(ec),1)]];
  end
%
% Find the centers:
  nTiles = length(llc);
  tc =  llc + ones(nTiles,1)*[750 750];
  
%
% Find the boundaries:
  tb = [tc(:,1)-1000 tc(:,1)+1000 tc(:,2)-1000 tc(:,2)+1000];
  
  if( exist('tiles.csv')  ~= 2 && exist(parentMapName) == 2 )
    %
    % Create .csv file:
    tilesFid = fopen('tiles.csv','w');
    fprintf( tilesFid, 'TileName , Easting , Northing , %d\n', nTiles );
    %
    % Build grdcut argument strings and loop through:
    gcutmap = sprintf('/usr/local/bin/gmt grdcut  ');
    for( i=1:nTiles )
      tileName    = sprintf('AxialTile%d', i);
      grdTileName = strcat(tileName,'.grd');
      boTileName  = strcat(tileName,'.bo');
      if( exist(grdTileName) == 2 )
	%
	% If it exists, don't recompute it.
	fprintf('Tile %s already exists; don''t recreate.\n',grdTileName);
	continue;
      end
      gcutNameArg = sprintf(' -G%s ', grdTileName);
      gcutBounds  = sprintf(' -R%d/%d/%d/%d -N ',tb(i,1),tb(i,2),...
			    tb(i,3),tb(i,4));
      gcutCmd = [gcutmap parentMapName gcutNameArg gcutBounds];
      fprintf('%s', gcutCmd);
      fprintf(tilesFid, '%s , %d , %d\n', boTileName, tc(i,1), tc(i,2)); 
      tic;
      system(gcutCmd);
      fprintf('Created tile %s.\n', boTileName);
      toc;
    end
  else
    fprintf('tiles.csv exists, skip tile generation loop.\n');
  end
  
  Levels = [-105:1:-50];
  tic
  fprintf('Plotting Map\n');
  %[cc hh] = contourf(plg.eastings,plg.northings,plg.z,Levels);
  deci = 10;
  [cc hh] = contourf(downsample(plg.eastings,deci),downsample(plg.northings,deci),downsample(downsample(plg.z,deci)',deci)');
  %[cc hh] = contourf(plg.eastings,plg.northings,plg.z);
  clabel(cc,hh);
  axis equal
  %set(gca,'xminorgrid','on')
  set(gca,'ytick',[nc nc(end)+1500]);
  %set(gca,'xminorgrid','on');
  set(gca,'xtick',[ec ec(end)+1500]);
  toc
  grid on
  hold on

  if( exist( 'missionFileName') );
    dataWP=corridor(['~/missions/',missionFileName]);
  end
  
  linewidth=1;
  if( exist( 'dataWP' ) == 1 & min(size(dataWP)) > 1 )
    plot( dataWP(:,3), dataWP(:,2), 'r', 'linewidth', linewidth);
    lstr(length(lstr)+1,1) = {'Trackline'};
    plot( dataWP(:,3), dataWP(:,2), 'gd','linewidth', linewidth);
    lstr(length(lstr)+1,1) = {'Waypoints'};
    sizeDataWP = size(dataWP);
  end
  
  for( i=1:nTiles )
    locStr = sprintf('%d',i);
    tHandle = text(tc(i,1),tc(i,2),locStr);
    set(tHandle,'color','m');
    set(tHandle,'fontsize',12);
  end
